Geometry and kinetics determine the microstructure in arrested coalescence of Pickering emulsion droplets

TL;DR

This study investigates how the shape and rate of particle coverage influence the final microstructure of arrested coalescence in Pickering emulsions, revealing the importance of extrinsic geometry and kinetics.

Contribution

It demonstrates that the microstructure depends on the ratio of area change rate to diffusion constant, highlighting the role of extrinsic geometry and arrest kinetics.

Findings

Microstructure depends on the ratio of area change rate to diffusion constant.

Extrinsic geometry influences the arrest process more than static packings.

Negative Gaussian curvature regions can form around the neck during coalescence.

Abstract

An important strategy to stabilize emulsions is to arrest coalescence of the constituent droplets with an opposing rheological force. Colloidal particles adsorbed on the surface of emulsion droplets in a Pickering emulsion become increasingly crowded during successive coalescence events because the combined surface area of coalescing droplets is less than that of the constituent droplets. Beyond a critical density, the particles form a rigid shell around the droplet and inhibit both relaxation of the droplet shape and further coalescence. The resulting droplets have a nonuniform distribution of curvature and, depending on the initial coverage, may incorporate a region with negative Gaussian curvature around the neck that bridges the two droplets. Here, we resolve the relative influence of the curvature and the kinetic process of arrest on the microstructure of the final state. Identifying the dimensionless ratio of the rate of area change \dot{A} to the diffusion constant D as a measure of the importance of kinetics in this system, we show that this depends on the extrinsic geometry of the surface as opposed to the static packings that depend solely on intrinsic geometry.

Figures

Equations

$$\frac{a^2{\lambda }^2(x^2+y^2)+a^2{z}^2}{(x^2+y^2+z^2{)}^2}=1,$$

$$\Delta p=2\gamma H$$

$$n_d=\frac{1}{2}\left(\frac{{\sum }_i|q_i|}{12}-1\right),$$

$$\frac{\dot{A}}{D}=\frac{dA}{dt}\frac{{\tau }_d}{2r^2}=\frac{1}{2r^2}\frac{{\tau }_d}{{\tau }_r}\frac{dA}{dT}$$

$$\theta ={\tau }_d/{\tau }_r$$

$$\dot{A}=\int \nabla \cdot \mathbf{N}\left(\frac{d\mathbf{X}}{d\lambda }\cdot \mathbf{N}\right)\frac{d\lambda }{dt}dA,$$

$$\dot{A}=\int 2H\left(\frac{d\mathbf{X}}{d\lambda }\cdot \mathbf{N}\right)\frac{d\lambda }{dt}dA,$$

$${\mathbf{x}}_i^{{}^{\prime }}(t+\Delta t_p)={\mathbf{x}}_i(t)+{\eta }_i\sqrt{2D\Delta t_p},$$

$$2r=\sqrt{2D{\tau }_d}.$$

$$V_{overlap}=\begin{cases}\begin{array}[]{cc}r^{2}-rx&x<r\\ 0&x\geq r\end{array},\end{cases}$$

Full text

I Introduction

Emulsions are mixtures of immiscible fluids where one fluid is dispersed as droplets inside the other, continuous fluid. Such mixtures phase separate or coarsen by coalescence, fusion of droplets, and Ostwald ripening, exchange of fluid between droplets by diffusion through the continuous phase. To stabilize an emulsion, it is therefore necessary to inhibit these mechanisms. One strategy to do so is to disperse micron sized or nanoscale particles in the mixture, forming a Pickering emulsion(Pawar et al., 2011; Studart et al., 2009; Dockx et al., 2018). The particles adsorb on the surface of the droplets, reducing the interfacial contact between the two fluids. When two such droplets come into contact, and the particle coverage fraction is less than some critical value $\Phi_{c}$, the droplets fully coalesce. Because the combined droplet has a lower surface area than the precursor droplets, however, the coverage fraction necessarily increases. Conversely, if droplets come into contact with a particle coverage that exceeds the critical coverage fraction, the particles become crowded and prevent the droplets from fully relaxing. The end result is typically a non-spherical droplet with a solid interfacial shell that prevents further coalescence of additional droplets and also tends to inhibit Ostwald ripening. This phenomenon is an example of arrested coalescence, which can also be achieved by other offsetting rheological resistance such as internal viscoelastic fluids(Pawar et al., 2012).

Prior studies(Pawar et al., 2011) of arrested coalescence have shown that the point of arrest can be predicted by a simple geometric model incorporating $\Phi_{0}$, the particle coverage fraction of the precursor droplets prior to coalescence, and the relative size of the two droplets. Once the expected coverage fraction for the combined droplet would exceed $\Phi_{c}=\pi/\sqrt{12}\approx 0.9$, the value of random close packing in 2D, the coalescence will be arrested; if $\Phi_{0}$ is increased, coalescence is arrested at an earlier point.

In addition to the benefit of increased stability, arrested coalescence also constitutes a straightforward method for sculpting non-spherical Pickering emulsion droplets. While these studies have successfully shown how to control the resulting shape of the emulsion droplets, little attention has been paid to the microstructure of the final state, or the role of the arrest process itself. The purpose of this paper is, therefore, to determine how the particle microstructure is affected by the nonuniform curvature, and to disentangle the relative influence of geometric influences and the kinetic process of relaxation prior to arrest.

Valuable insight into the role played by geometry comes from the rich literature on packings on curved surfaces(Bowick and Giomi, 2009), for which spherical Pickering emulsion droplets—dubbed *colloidosomes—*have proven an ideal model system(Bausch et al., 2003). Unlike 2D Euclidean space where crystalline hexagonal packing is optimal, on curved surfaces defects—particles with a contact number $c_{i}$ other than $6$—are generally required to accommodate the curvature. Defining the defect charge $q_{i}=c_{i}-6$, the total charge of the curved surface must satisfy $\sum_{i}q_{i}=6\chi$, where $\chi$ is the Euler characteristic of the surface.

Continuum elastic theory(Seung and Nelson, 1988) can be used to analyze the defect distribution on curved surfaces(Pérez-Garrido et al., 1997; Bowick et al., 2000, 2002, 2006; Vitelli et al., 2006; Giomi and Bowick, 2007, 2008; Irvine et al., 2010). The free energy is dominated by the competition between the defect core energy term, which favors fewer defects, and the compensating elastic term that arranges the defect charge distribution to approximately match the Gaussian curvature. For packings on spherical surfaces where $\chi=2$, if the core energy is large, i.e. the number of particles is small, there are only the twelve 5-fold disclinations required by the Euler Theorem which form an icosahedron. Conversely, if the core energy is small, and hence the number of particles is large, excess dislocations appear with different signs to maintain the total defect charge required by topology. The excess dislocations tend to proliferate into grain boundaries, referred to as scars, to help screen the Gaussian curvature(Bowick et al., 2000). Experiments have confirmed that the extra disclinations start to emerge when the number of particles exceeds some limit(Bausch et al., 2003). Nonuniform curvature leads to localization of the defects(Irvine et al., 2010; Brojan et al., 2015; Jiménez et al., 2016; Burke and Atherton, 2016).

Insights from the packing literature are most applicable when the particle structure is in quasistatic equilibrium with the host shape, i.e. when relaxation of the droplets formed in arrested coalescence occurs much more slowly than processes such diffusion that relax the structure. We are, however, unaware of any prior attempt to understand the role of kinetics in arrested coalescence. To that end, we simulate as described in section II the structures produced by arrested coalescence using previous theoretical and experimental studies of liquid drop coalescence(Studart et al., 2009; Pawar et al., 2011; Eggers et al., 1999; Menchaca-Rocha et al., 2001; Wu et al., 2004). We then separate the influence of geometry and kinetics on the microstructure by comparing static packings produced with a fixed shape in section III and kinetic packings produced as the shape relaxes towards the final spherical ground state in section IV. Properties of the surface that predict these effects are described together with conclusions in section VI.

II Model

An analytical ansatz to describe the shape evolution of a pair of spheres as they relax following coalescence was proposed by Garabedian et al.(Garabedian and Helble, 2001). After initial contact, the surface is described by the level set,

[TABLE]

where $a$ is the half length of long axis and aspect ratio $\lambda\in\left[0,1\right]$ controls the extent of coalescence as shown in Fig. 1A. The value $\lambda=0$ corresponds to the two droplets just touching each other and $\lambda=1$ represents the final state as one spherical droplet. The center is located at the origin and the $z$ axis is the axis of rotational symmetry. The value $a$ is chosen for each $\lambda$ to hold the total volume of the surface constant.

The defect distribution in static packings is known to be controlled by the distribution of Gaussian curvature, which acts like a nonuniform background charge distribution in addition to the discrete defect charges for the elastic term in free energy.(Giomi and Bowick, 2008). We therefore display the Gaussian curvature along the axis of rotation $z$ for several different values of $\lambda$ with fixed volume as the sum of two spheres with diameter 1 in Fig. 1B. During the initial stages of relaxation, i.e. $\lambda\lesssim 0.7$, the neck of the doublet creates a region where the Gaussian curvature $K$ is negative. We show in the inset of Fig. 1B the transition of curvature at the neck from negative to positive. For $\lambda=0.3$, the curvature at the neck is large and negative while at the ends of the droplet the curvature approaches a constant value because here the surface is almost spherical. As $\lambda$ increases, the profile $K(z)$ becomes smoothed over time until it approaches a uniform constant value for the spherical final state.

Also shown in Fig. 1C is the Mean curvature $H$, which similarly exhibits a negative region at the neck for early stages of the relaxation. This quantity is related to the capillary pressure difference across the surface $\Delta p$ through the Young-Laplace equation,

[TABLE]

where $\gamma$ is the surface tension. As a measure, therefore, of the local generalized force acting to evolve the shape towards the equilibrium state, we shall show in section IV that the mean curvature plays a key role in the kinetically dominated regime.

To simulate particles embedded on the surface we employ two different algorithms fully described in the Methods section below. Static packings are produced on shapes of fixed $\lambda$ using a Monte Carlo inflation algorithm inspired by the Lubachevsky–Stillinger algorithm(Lubachevsky and Stillinger, 1990); we supplement this algorithm to ensure rigidity of the final packings. A second algorithm, with particles of fixed radius $r$ that diffuse with their centers of mass constrained to the evolving ansatz surface, is used to resolve the influence of kinetics. While the particles move by diffusion, the shape $\lambda(t)$ is slowly evolved at constant volume such that the radius of the neck scales $\propto t^{1/2}$ with time. This form was proposed for the inertial regime where Reynolds number is large by Eggers et al.(Eggers et al., 1999) and confirmed experimentally(Menchaca-Rocha et al., 2001; Wu et al., 2004) to hold in the early stages of coalescence. We use this power law for the whole process for simplicity.

III Geometry

Using the protocol described above with fixed shape, we generated 100 rigid packings of 1000 particles for several values of the coalescence parameter $\lambda\in[0.3,1]$. Representative results are shown in Fig. 2A, where the particles are colored by coordination number computed from the Delaunay triangulation. Visually, the packings appear largely crystalline with the expected chains of defects distributed over the whole surface. Plots of the density-density pair correlation function $g(s)$(Donev et al., 2005) as a function of arclength $s$ in units of the particle radius $r$ for the particles around the neck ($z=0$) and pole ($z=a$) are shown in Fig. 2B as dashed and solid lines respectively and appear similar. In agreement with visual inspection, these plots show long range order and split second peaks indicative of local crystalline order.

To visualize the scar network, we compute the neighbor graph from the Delaunay triangulation and delete all vertices with six-fold coordination(Mascioli et al., 2017). The disjoint components of this defect subgraph correspond to individual scars; typical scar networks are shown in Fig. 2A alongside the corresponding packing. The majority of the scars display a linear morphology with some branching. No strong variation in scar morphology is observed for different $\lambda$, however, the size distribution of the scars does vary. In 2A, the largest scar gets longer as $\lambda\to 1$. The scar length, here defined as the mean number of vertices in the disjoint subgraphs, is plotted in Fig. 2C. For low $\lambda$ packings, shorter scars are generated, and as shown in the inset of Fig. 2C, a higher fraction of them are isolated singleton defects. This is because for low $\lambda$ the radius of the two spheres is smaller and hence the effect of curvature is stronger. As the negative curvature becomes less distinct, the scar length increases until $\lambda\thickapprox 0.7$, the point at which the region of negative Gaussian curvature disappears. Beyond this value, the length of the scars doesn’t increase further, showing that positive Gaussian curvature generates similar scars.

Consistent with prior work, we define $n_{d}$, the* *excess dislocations beyond the twelve required by topology as(Burke et al., 2015),

[TABLE]

where the sum is over dislocations and $q_{i}$ represents the charge of the $i$th dislocation, and display $n_{d}$ as a function of $\lambda$ in Fig. 2D. For low $\lambda$, corresponding to early arrest, the large amount of negative curvature induces more dislocations. As the neck region becomes flatter, the defect number decreases and reaches a minimum at $\lambda=0.6$. For $\lambda>0.6$, more dislocations emerge to accommodate the positive curvature until the final spherical stage. The packing fraction $\phi$, defined as the fraction of the surface covered by particles, is also shown in Fig. 2D and follows an inverted trend because the emergence of dislocations reduces the packing fraction. Since the shape away from the neck is spherical, this variation in the $n_{d}$ appears mainly due to the curvature distribution around the neck.

To test this, spatially resolved plots of the bond orientational order parameter $\psi_{6}=\left\langle\exp(i6\theta)\right\rangle$(Nelson and Halperin, 1979) for several $\lambda$ are displayed in Fig. 2E. In these plots, the surface is divided into 24 equal-area axially symmetric regions and the symmetry of the shape is exploited by collapsing corresponding regions for positive and negative $z$. The orientational order is generally quite high, greater than $0.8$ everywhere except close to the neck for $\lambda=0.3$, and as expected is spatially uniform for the spherical case $\lambda=1$. The distribution of $\psi_{6}$ undergoes a transition around the neck area: for $\lambda=0.3$, the region close to the neck has low $\psi_{6}$, with a local maximum at around $z/a=0.2$. In contrast, for $\lambda=0.6$ the neck region has enhanced $\psi_{6}$, and a less ordered region away from the neck. As $\lambda$ increases further, the amplitude of this variation is reduced, with lower order around the neck.

In the inset of Fig. 2E, the mean value of $\psi_{6}$ at the innermost region is shown as a function of $\lambda$. As $\lambda$ increases, the neck becomes less pronounced, and the order increases; after $\lambda\sim 0.6$, when the neck becomes flat, the appearance of the positive curvature reduces the order again.

Consistent with this, the number density of the defects, shown in Fig. 2F, closely resembles the inverted form of $\psi_{6}$ since defects reduce the bond orientational order dramatically. More defects are induced by the large negative curvature around the neck while the flat region has fewer defects and the number of defects increases again as the neck gains positive curvature. As $\lambda\to 1.0$, a uniform distribution of defects occurs due to the uniform curvature.

Finally, in Fig. 2G and its inset, the defect charge density is shown in solid lines, together with a predicted distribution from the integrated Gaussian curvature in dashed lines. The defect charge distribution is well explained by the curvature distribution of the surface and also the constraint that the total charge must be $12$ from the Euler characteristic of the surface. For $\lambda=0.3$, strongly negative Gaussian curvature at the neck as apparent in Fig. 1B tends to induce negative defects. Hence, additional positive charges must be generated at the ends to satisfy the topological constraint. Conversely, for $\lambda=0.6$, fewer negative defects are induced and the neck region has almost zero mean curvature leading to an overall enhancement of the order. Defects are needed outside the neck area to satisfy the charge constraint and also to match the positive curvature.

This section has elucidated the role played by the nonuniform geometry on the final states observed in arrested coalescence, notably the key role that the neck plays in determining the distribution of defects. For shapes corresponding to early arrest, when the neck possesses strong negative curvature, additional defects are induced, while at shapes corresponding to intermediate arrest, the flat neck region promotes locally enhanced order. Scars at early arrest are shorter, because the radius of curvature at the caps is smaller; the scar length grows until the region of negative curvature vanishes. Consistent with previous work, the defect charge distribution is well predicted by the distribution of Gaussian curvature and the overall charge constraint imposed by topology.

IV Kinetics

We now examine the role played by kinetics on the arrested state. A natural measure of the importance of kinetics is the rate of total area change $\dot{A}$ at the arrest point, which may be compared to the diffusion coefficient $D$ because they have the same dimensions. We therefore construct a new dimensionless parameter,

[TABLE]

and nondimensionalize time using $\tau_{d}$ the corresponding diffusion timescale defined through $2r=\sqrt{2D\tau_{d}}$, $T$ is the dimensionless time $t/\tau_{r}$ and $\tau_{r}$ the time of full relaxation. The ratio

[TABLE]

therefore also emerges as a measure of the relative influence of diffusion and relaxation and is in practice the independent variable that is varied to construct our ensemble of simulations. We shall also make use of the quantity $\dot{A}/(ND)$ where $N$ is the number of particles to quantify the average rate of area change per particle, together with a local version $\frac{\dot{a}}{nD}$, where $a$ is the area of some region of interest on the surface and $n$ is the number of particles in that region. The quantity $dA/dT$ is almost linearly related to $\lambda$ for all but small values of $\lambda<0.1$ where the ansatz is a poor approximation to the experiment(Pawar et al., 2011), and so we do not investigate arrest in this regime.

Using the protocol described in section II, a set of arrested states was generated with $N=800$ particles and suitable particle radius to promote a point of arrest $\lambda_{a}$ for two different scenarios: an early arrest case with $\lambda_{a}\approx 0.3$, where the curvature of the neck is extremely negative and late arrest at $\lambda_{a}\approx 0.6$, where the neck is flat. The ratio of diffusion time scale to the total relaxation time $\theta=\tau_{d}/\tau_{r}$ is varied from $2^{-8}$ to $2^{4}$ with 50 samples for each value.

In Fig. 3A, the bond orientational order parameter $\psi_{6}$ is shown as a function of $|\dot{A}|/ND$, which is, as was shown in Eq. 4, linearly related to $\theta$. This choice of parameter accounts for the fact that the rate of area change is naturally larger for early arrest. Faster relaxation leads, as might be expected, to less ordered structures for both scenarios, and overall the order is higher for late arrest. This latter is consistent with the results of Fig. 2E that show more ordered packings for the $\lambda=0.6$ surface and is therefore geometric in origin.

Before examining microstructural differences between the two scenarios, we pause to discuss the effect of the shape evolution on the particle motion. A single free particle at rest on the initial surface moves only subject to constraint forces which lie along the direction locally normal to the surface. We neglect the inertia of the particles which is assumed to be damped by the surrounding fluid. Representative trajectories are visualized in Fig. 3C for several starting positions around the doublet. In some locations, for example the caps, the trajectories are convergent indicating that particles here tend to be compressed as the doublet relaxes; elsewhere the trajectories are weakly divergent indicating that here particles would be pushed apart.

We can use a similar approach to define a parametrization-independent local rate of change in area $\dot{a}$ by dividing the initial surface into annuli of equal area and evolving the boundaries along the normal direction. The instantaneous $\frac{da}{dT}$ as a function of position is plotted in Fig. 3B for the early and late arrest scenarios, where $T$ is the dimensionless time $t/\tau_{r}$. An important difference is that for early arrest $\lambda_{a}\approx 0.3$, a very strong region of compression exists close to the neck; this is missing for $\lambda_{a}\approx 0.6$ where particles at the neck are being pushed apart by the constraint forces. Compressive regions exist in both scenarios around the caps.

Identical particles initially uniformly distributed on the doublet therefore evolve differently depending on their position. Those near the caps are pushed together and should crowd first, suggesting that the crystalline region should develop first at the caps and grow, reminiscent of the growth fronts observed in (Waitukaitis and Jaeger, 2012; Han et al., 2016). If the arrest is early, crowding may also occur around the neck. Particles in the intermediate region between the neck and the cap have more opportunity to relax than particles in the compressive regions.

The spatially resolved bond orientational order parameter is shown for both scenarios in Fig. 4A and B for $\lambda\approx 0.6$ and $\lambda\approx 0.3$ respectively. Different traces correspond to different values of $|\dot{A}|/ND$ and, as before, the left and right portions of the doublet are combined into one plot. In both cases, the static order parameter distribution, shown by the red line, is recovered for sufficiently slow relaxation. For late arrest, this corresponds to an enhancement of $\psi_{6}$ in the flatter central area and reduced order closer to the cap. As $|\dot{A}|/ND$ increases, these features are softened and as $|\dot{A}|/ND\gtrsim 1$, the orientational order converges on a uniform constant value of $\sim 0.83$ which is the minimum order achieved by a static packing with similar shape.

In contrast, for early arrest, the distribution of $\psi_{6}$ remains similar as a function of $|\dot{A}|/ND$: there is a strong reduction in the order parameter around the neck and a uniform distribution near the cap. Increasing $|\dot{A}|/ND$ reduces the order globally, shifting the curves down by as much as $0.1$ but does not change their overall form.

To understand the role of the arrest point more carefully, an ensemble of simulations with different arrest points $\lambda$ was run in the kinetically dominated regime $|\dot{A}|/ND\approx 10$ and the resulting distributions of $\psi_{6}$ are shown in Fig. 4C. Comparing these with equivalent plots for static packings, Fig. 2E, we see that for arrest earlier than $\lambda\sim 0.5$, the static and dynamic order distributions appear similar, although the overall values of $\psi_{6}$ are attenuated. For arrest after $\lambda\sim 0.5$, the distributions no longer resemble the static distributions, instead having approximately uniform $\psi_{6}\approx 0.83$. Rapid relaxation therefore appears more readily able to wash out variations caused by the nonuniform curvature for late arrest.

The transition in behavior appears to coincide with the point at which the compressive region disappears from the neck, as can be seen in Fig. 4D which displays the instantaneous local dynamical parameter $\dot{a}/nD$ at the neck at the point of arrest. For early arrest, therefore, the signature of Gaussian curvature remains imprinted on the microstructure regardless of the overall rate of relaxation.

Time resolved analysis of the orientational order parameter allows us to understand the transition further. The particle number density is plotted at several time points $t/t_{f}$ in Fig. 5A for the late arrest scenario at $\lambda_{a}\approx 0.6$, showing growth of a dense region from the cap and neck as predicted from the geometric argument above. Corresponding plots of the the $\psi_{6}$ distribution are displayed in Fig. 5B, and some simulation snapshots are visualized in Fig. 5E. Here $t_{f}$ is defined to be the time at which arrest occurs. Corresponding plots and visualizations for early arrest $\lambda_{a}\approx 0.3$ are displayed in Fig. 5C, D and F respectively. Neighbors of particles are not determined from Delaunay triangulation since the particles are generally not densely packed during the coalescence. Instead, we define neighbors as particles within the center-to-center distance of $1.5$ times diameter. In practice, this definition changes the numerical values of $\psi_{6}$ very little and doesn’t affect the trend of evolution.

For late arrest particles are initially uniformly ordered, except close to the neck. As relaxation proceeds, compression at the cap produces a denser region that converges to about $\psi_{6}\sim 0.83$. This value is similar to the ordering observed in the static packing at the cap. At the neck, the particles are initially more disordered because of the strong concentration of negative Gaussian curvature, and the initial compressive region that exists here begins to widen the disordered region. At around $\lambda\sim 0.4$, however, the widening stops, and as $\lambda$ increases further the ordered region growing in from the caps completely overcomes the disordered neck region.

For early arrest, the ordering at the cap changes very little through the course of the simulation, remaining similar to the static value. The initially disordered region at the neck begins to widen, as it did for the late arrest scenario, and reaches a maximum width at around $\lambda=0.25$. After this, the neck region begins to uniformize but does not completely do so, freezing in the disordered region that is also seen in the static case.

The influence of kinetics on the microstructure therefore depends critically on the arrest point: the history of compression and expansion of the surface may freeze in disordered regions caused by the underlying curvature, even where relaxation might be expected to wash out the effect of geometry. For early arrest, the strong initial compression leads to a disordered region that is frozen in at all relaxation speeds tested, while for late arrest fast relaxation fully blurs out the nonuniform distribution of order promoted by geometry.

V Discussion

The static packings produced in section III reconfirmed the well-known role of the Gaussian curvature in determining the distribution of order. The role of kinetics was shown in section IV to be surprisingly complicated, with a disordered region present in the neck for early arrest, while for late arrest kinetics is able to fully wash out variations in the distribution of order caused by geometry. The final microstructure therefore depends critically upon the history of compressive and and expansive regions as the relaxation proceeds. These regions are identified by examining local contributions to the the rate of change of the area $\dot{A}$.

We can compute $\dot{A}$ from more fundamental quantities as follows. Suppose we have a one-parameter family of surfaces $\mathbf{X}(\lambda)$ that describes the shape evolution of the surface. The rate of change of area can then be written,

[TABLE]

where $\mathbf{N}$ is the local outward surface normal and the integral is over the surface. The first factor is the divergence of the normal and can be rewritten as the mean curvature,

[TABLE]

while the second factor in (6) and (7) is the normal component of the velocity as the surface evolves according to $\lambda(t)$. We identify the integrand of (7) as $\dot{a}$. We note that (7) reveals an elegant distinction: the microstructure of the static packings is determined by the Gaussian curvature, an intrinsic quantity, while the role of kinetics is determined by extrinsic quantities—those that depend on how the surface is embedded—such as the mean curvature.

The discriminant $\text{sign}(\dot{a}/D)$ that identifies regions of compressed and expansion furthermore includes both properties of the surface at arrest and quantities that depend on its evolution. Importantly, the discriminant arises both from the sign of the mean curvature and whether the motion is locally inward or outward. At the caps, $H$ is positive but the motion is inward; conversely at the neck $H$ is negative, at least for low values of $\lambda$, but the motion is outward. One can therefore imagine manipulating, through a judicious choice of interface, both of these quantities to exploit kinetics as a means to control the microstructure. We note that the evidence of the previous section suggests that, at least crudely, it is the value of these quantities immediately prior to arrest that are most important for determining the microstructure: for the late arrest scenario $\lambda_{a}\sim 0.6$, the compressive regions that exist early in the relaxation appear to have little influence on the final state.

VI Conclusion

This work aimed to distinguish the relative contribution of geometry and kinetics in determining the microstructure of droplets produced by arrested coalescence. Relaxation of the shape drives the transition to a final rigid structure, while motion of the particles, here modelled as diffusive, facilitates reordering of the structure to achieve higher packing fraction. The importance of each processes was shown to be quantified by a single parameter $\dot{A}/ND$ which compares the rate of change of area experienced by each particle to its diffusion constant. As $\dot{A}/ND\to 0$, the microstructure closely resembles the static packing scenario where the system distributes defects to match the Gaussian curvature. At finite $\dot{A}/ND$, kinetics tends to reduce the overall order of the arrested structure, and tends to blur out modulation of the order as the structure can no longer be fully relaxed. Disordered regions induced by the Gaussian curvature can, however, remain even in the kinetically dominated regime because, as the surface relaxes, regions of compression and expansion serve to trap them. The kinetic influence therefore, perhaps surprisingly, depends on geometric properties of the hypersurface including mean curvature but these are extrinsic in origin, i.e. depending on the embedding.

The remarkably rich influence of kinetics suggests the possibility of exploiting it as a means to control the microstructure of colloidosomes. One can speculatively envisage designing a shape using the formula (7) that incorporates compressive regions prior to arrest and could selectively lock in disordered regions; these might become targets for further coalescnce events in multistage assembly. The extent to which one is free to choose the mean curvature and the local velocity of the surface remains an open problem, however. The Young-Laplace equation Eq. (2) implies that locally these are not independent quantities. That they are different in the present work shows that constraints such as the overall volume play an essential role, and the consequent restrictions on the design space remain to be illuminated.

A second direction that should be pursued is the connection of these arrested shapes to jamming, a transition to rigidity as a function of density that occurs in particulate media(Torquato and Stillinger, 2010; Liu and Nagel, 2010). While the kinetically arrested structures observed here are not jammed, in that they may possess unconstrained collective motions of particle, the static packings we use as a comparison are because this is explicitly enforced. Two of the present authors recently proposed that rigid structures formed as a result of shape evolution form a new class called* “metric jamming”* (Burke and Atherton, 2016) where the final state is rigid both with respect to perturbations of the particles and the manifold on which they are embedded. Analysis of the arrested coalescence problem along these lines may help determine the longevity of the undoubtedly metastable arrested structures, as well as provide tools to determine the mechanical properties of these structures.

VII Methods

**Static packings—**Particles are initially dispersed with their center of mass on the surface at zero radius, diffused by Brownian motion according to the Langevin equation,

[TABLE]

where $\mathbf{\eta}_{i}$ is a random step drawn from Gaussian distribution along the tangent plane, and $D$ is the diffusion constant such that the variance of stepsize for Brownian motion in time $t$ is $2Dt$. We may therefore define a characteristic diffusion time scale $\tau_{d}$ that gives a standard deviation of stepsize equal to the particle diameter,

[TABLE]

As the particles diffuse, their radii $r$ are increased (inflation moves) very slowly, with $\frac{\delta r}{\sqrt{2Dt}}\sim 10^{-4}$ in unit time. Collective motions that undo overlap are found at each stage by gradient descent on an artificial potential,

[TABLE]

that penalizes overlap. The simulation is halted when no further move is possible without inducing overlaps.

Generically, packings produced by this algorithm need not be rigid, i.e. there may exist collective motions of particles that can unjam the system and allow further relaxation of the surface. We therefore adapt(Burke and Atherton, 2016) a linear programming approach (Donev et al., 2004) to identify these collective motions, execute them, and restart the packing simulation. Before applying the linear program, the configuration is conditioned by minimizing an artificial soft repulsive potential imposed between all pairs of particles; this tends to push the particles away from one another. This process is repeated until a rigid final state is achieved.

**Dynamic simulations—**A second algorithm was used to understand how the relaxation process affects the final structure. For these simulations, particles are initially dispersed by random sequential deposition with a fixed particle radius $r$. During the simulation, diffusion moves are made as before in Eq. (8). During relaxation moves, the particles are constrained to the surface with overlaps prevented to first order using Lagrange multipliers(Vesely, 2013). After each relaxation step, Eq. (10) is minimized to remove overlaps. If not all overlaps could be undone, the timestep is reduced. The algorithm halts when the timestep of relaxation is smaller than a threshold $\delta t$.

Acknowledgements.

This material is based upon work supported by the National Science Foundation under Grant No. DMR-1654283.