# Numerical analysis of Pickering emulsion stability: insights from ABMD simulations

###### Abstract

The issue of the stability of Pickering emulsions is tackled at a mesoscopic level using dissipative particle dynamics simulations within the Adiabatic Biased Molecular Dynamics framework. We consider the early stage of the coalescence process between two spherical water droplets in decane solvent. The droplets are stabilized by Janus nanoparticles of different shapes (spherical and ellipsoidal) with different three-phase contact angles. Given a sufficiently dense layer of particles on the droplets, we show that the stabilization mechanism strongly depends on the collision speed. This is consistent with a coalescence mechanism governed by the rheology of the interfacial region. When the system is forced to coalesce sufficiently slowly, we investigate at a mesoscopic level how the ability of the nanoparticles to stabilize Pickering emulsions is discriminated by nanoparticle mobility and the associated caging effect. These properties are both related to the interparticle interaction and the hydrodynamic resistance in the liquid film between the approaching interfaces.

^{†}

^{†}thanks: Corresponding author: francois.sicard@free.fr.

## I Introduction

Pickering emulsions Pickering (1907), i.e., particle-stabilized emulsions, have been used in various
applications including biofuel processing Drexler et al. (2012), encapsulation for drug delivery Angelova et al. (2011),
and food preservation Shchukina and Shchukin (2012); Puglia and Bonina (2012).
When nanoparticles (NPs) are used as emulsifiers, their characteristics might have critical effects
on properties such as interfacial tension Miller et al. (2006), droplet size,
and emulsion stability Wang et al. (2011); Zahn et al. (1997); Cheung (2010); Weeks and Weitz (2002); Sonnenburg et al. (1991); Peng et al. (2009); Cheng and Grest (2012).
Although the factors affecting Pickering emulsion stability are still a key topic of research,
the density of NPs on the droplet surface is known to be among the most important parameters X-C.Luu et al. (2013).
Indeed, the interfacial tension of the droplet depends strongly on the NP coverage X-C.Luu et al. (2013),
which also depends on the NP affinity to the interface, i.e., the NP desorption energy.
When the NP density on the droplet is large enough, the three-phase contact angle discriminates
between NPs that are effective at stabilizing emulsion and those that are not FanSM2012 (); Fan and Striolo (2012).
It is known that particles, once adsorbed, are rather difficult to displace from the interface.
It is this property that makes particles such a good colloidal stabilizer for emulsions and foams.
Besides the structure of NP monolayers, the NP diffusion is also believed to affect the system stability.
Janus nanoparticles are believed to be more effective than homogeneous NPs in the stabilization of emulsion FanSM2012 ().
Indeed, Janus NPs typically have high adsorption energies, and therefore are expected to pack densely at an interface.

In our present work, we study the early stage of the coalescence process between two spherical water droplets in decane solvent. The droplets are stabilized with various Janus NP types. We implement Adiabatic Biased Molecular Dynamics (ABMD) simulations Paci and Karplus (1999); Camilloni et al. (2011); Marchi and Ballone (1999), where the system is forced to coalesce adiabatically. We show that, for a sufficiently dense layer of particles adsorbed on water droplets, the stabilization mechanism can strongly depend on the collision speed achieved during the merging process. For collision velocities sufficiently slow compared to the rate of drainage of the thin liquid film trapped between the coalescing emulsion droplets, the ability of the NPs to stabilize Pickering emulsions is dictated by NP mobility and the associated caging effect. This coalescence process is compared with simulations at sufficiently high collision speed, which suggest the droplet faculty to undergo large elastic deformation before coalescing is the stabilization mechanism.

## Ii Methods and algorithms

### ii.1 MD simulation

The Dissipative Particle Dynamics (DPD) simulation method Groot and Warren (1997) was implemented within the simulation package LAMMPS Plimpton (1995). The procedure and the parametrisation details are fully described in prior work Luu et al. (2013); X-C.Luu et al. (2013). The system simulated here is composed of water, oil (decane), and NPs. One ”water bead” (w) represents 5 water molecules. One decane molecule is modeled as two ”oil beads” (o) connected by one harmonic spring of length and spring constant Groot and Rabone (2001), where is the DPD cutoff distance. The size of the simulation box is , where is the box length along the direction. Periodic boundary conditions are applied in all three directions. The NPs are modelled as hollow rigid ellipsoids and contain polar (p) and nonpolar (ap) DPD beads on their surface. One DPD bead was placed at the NP center for convenience, as described elsewhere Luu et al. (2013); X-C.Luu et al. (2013). Hollow models have been used in the literature to simulate NPs, and hollow NPs can also be synthesized experimentally Calvaresi et al. (2009). We considered spherical and ellipsoidal NP with shape defined by the equation , where , , and are Cartesian coordinates and and are the semi-axes of the ellipsoid. When , spherical NPs are obtained. When , the ellipsoid is oblate; when , it is prolate. All NPs simulated had the same volume, , where is the radius of the equivalent sphere. We imposed nm. All types of beads in our simulations have reduced mass of . The total number of beads on one NP surface changes with the aspect ratio . This allows us to maintain the surface bead density sufficiently high to prevent other DPD beads (either decane or water) from penetrating the NPs (which would be unphysical), as it has already been explained elsewhere X-C.Luu et al. (2013). To differenciate every NPs, we report the nonpolar fraction of the NP surface beads. For example, indicates that of the beads on the NP surface are nonpolar. The interaction parameters shown in Table 1 are used herein. These parameters were adjusted to reproduce selected atomistic simulation results, as explained in prior work Luu et al. (2013). By tuning the interaction parameters between polar or nonpolar NP beads and the water and decane beads present in our system, it is possible to quantify the effect of surface chemistry on the structure and dynamics of NPs at water-oil interfaces.

All simulations were carried out in the NVE ensemble Luu et al. (2013). The scaled temperature was , equivalent to K. The DPD time scale can be gauged by matching the self-diffusion of water. As demonstrated by Groot and Rabone Groot and Rabone (2001), the time constant of the simulation can be calculated as , where is the DPD time constant, is the simulated water self-diffusion coefficient, and is the experimental water self-diffusion coefficient. When (cf. Tab. 1), we obtained . For Partington et al. (1952), we finally obtain ps. In the simulations discussed herein, the size of the droplet was fixed. At the beginning of each simulation, the solvent (oil) beads were distributed within the simulation box forming a cubic lattice. Two water droplets of radius were generated by replacing the oil beads with water beads within the volume of two spherical surfaces separated with a distance of .

A number of spherical or ellipsoidal NPs were placed randomly at the water-decane interface with their polar (nonpolar) part in the water (oil) phase to reach the desired water-decane interfacial area per NP. Following previous work Luu and Striolo (2014). the NPs considered in this study are oblate with an aspect ration when spherical and when ellipsoidal. Note that nonspherical NPs have been considered as it has been reported they can be more efficient in stabilizing emulsions than spherical NPs. Indeed, it has been reported that the minimum amount of particle needed to stabilize an emulsion decreases as the particle aspect ratio increases Mejia et al. (2012); Dendukuri et al. (2006). The initial configuration obtained was simulated for timesteps in order to relax the density of the system and the contact angle of the nanoparticles on the droplet. The system pressure and the three-phase contact angles did not change notably after 5000 simulation steps.

144 NPs | 160 NPs | 165 NPs | 170 NPs | ||
---|---|---|---|---|---|

S-55JP | |||||

S-60JP | |||||

S-70JP | |||||

E-70JP | |||||

### ii.2 Droplet characteristics

It is well established that the stability of Pickering emulsions is chiefly determined by the size, topology, and concentration of particles at the interface, as well as their wetting properties FanSM2012 (); X-C.Luu et al. (2013); Luu and Striolo (2014); Qi et al. (2014); French et al. (2015). In the following we consider a system made by two identical spherical water droplets immersed in oil, and stabilized by a sufficiently dense layer of NPs. By sufficiently dense layer of NPs, we mean the NP concentration must be dense enough that the liquid film drainage time, i.e., the time interval from droplet contact to fusion, is significant. This can be tackled in a first instance considering steered simulations with a range of pulling rates representative of shear rates that can be reached in a conventional mechanical emulsification Mason and Bibette (1996) or ultrasonication processes Delmas et al. (2011). In the following, we consider a limited range of NP concentrations and types (i.e. contact angles) as well as different shape (spherical and ellipsoidal), as initial conditions. Table 2 summarizes the thermodynamic (three-phase contact angle, NP orientation angle) and the geometrical properties Vymetal and Vondrasek (2011) (radius of gyration, asphericity, and relative shape anisotropy of the droplet) for different NP types and concentrations that we consider throughout this study.

Three-phase contact angle. To estimate the three phase contact angle on the droplets we calculate the fraction of the spherical NP surface area that is wetted by water, and we obtained as

(1) |

where is the area of the NP surface that is wetted by water and is the radius of the NP. The ratio is obtained by dividing the number of NP surface beads (ap or p), which are wetted by water, by the total number of beads on the NP surface ( for spherical NP). One surface bead is wet by water if a water bead is the solvent bead nearest to it. One standard deviation from the average is used to estimate the statistical uncertainty.

Orientation angle. To characterize the ellipsoidal NPs structure on a droplet, we focus on the orientation of their longest axes (the axes) with respect to the radial directions from the center of the droplet Luu and Striolo (2014). The orientation angle is defined as

(2) |

where is the unit vector along the NP axis and is the unit vector representing the radial direction from the center of the droplet. When , the correspondent NP is parallel (perpendicular) to the droplet radial direction, and therefore perpendicular (parallel) to the local liquid-liquid interface. Note that we define the angle so that it is never larger than . One standard deviation from the average is used to estimate the statistical uncertainty.

NPs | NPs | NPs | NPs | |

S-55JP | ||||

S-60JP | ||||

S-70JP | ||||

E-70JP |

Radius of gyration, asphericity and relative shape anisotropy. The description of the geometrical properties of complex systems by generalized parameters such as the radius of gyration or principal components of the gyration tensor has a long history in macromolecular chemistry and biophysics Vymetal and Vondrasek (2011); Solc (1971). Indeed, such descriptors allow an evaluation of the overall shape of a system (spherical, oblate, or prolate) and reveal its symmetry. Considering, e.g., the following definition for the gyration tensor,

(3) |

where the summation is performed over atoms and the coordinates , , and are related to the geometrical center of the atoms, one can define a reference frame where can be diagonalized:

(4) |

In this format we obey the convention of indexing the eigenvalues according to their magnitude. We thus define the radius of gyration , the asphericity , which measures the deviation from the spherical symmetry, and the relative shape anisotropy , which reflects both the symmetry and dimensionality of a system Vymetal and Vondrasek (2011).

As discussed in the Electronic Supporting Information (ESI), we consider NP types and concentrations such that they are strongly effective at preventing droplet coalescence, and such that the shape of the droplet be spherical before contacting the neighbour droplet. In the simulation discussed herein, we consider spherical water droplets in decane solvent stabilized with 160-165 NPs.

### ii.3 ABMD algorithm

Adiabatic biased molecular dynamics Marchi and Ballone (1999); Camilloni et al. (2011); Paci and Karplus (1999) (ABMD) is an algorithm developed to connect any two points in the conformational space of a given system. The method is based on the introduction of a biasing potential, which is a function of chosen coordinates of the system and which is zero when the system is moving toward the desired target point, while disfavoring motions in the opposite direction. This is similar to what happens in a ratchet-and-pawl system, which undergoes random thermal fluctuations, while the pawl allows the ratchet to move only in one direction Camilloni et al. (2011). In ABMD the chosen coordinate plays the role of the ratchet and the biaising potential that of the pawl. In this respect, the biasing potential does not exert any work to direct the system toward the target conformation, as it would happen when pulling the system with a force; on the contrary the system makes moves toward the target conformation under the driving effect of the potential of the force-field alone. The biasing potential is implemented as

(5) |

where

(6) |

is the distance along the coordinate of the actual configuration of the system with respect to a target value , is a damping constant, and

(7) |

is the minimum distance reached until time . The algorithm is defined by the choice of the coordinate and of the damping constant . The ABMD algorithm is routinely implemented by the biophysical community to explore free-energy paths associated to protein folding and unfolding Paci and Karplus (1999); Provasi and Filizola (2010); Camilloni et al. (2011). In such applications, the most probable sequence of events generated by the ABMD simulations is expected to coincide with the minimum free-energy path independently on the damping constant. The outcome of a set of ABMD simulations is a number of trajectories mapping the free-energy minimum of the system, in which time and the statistical weight of the conformations along the trajectories are unphysically modulated. One of the unique features of this method is that the transformation connecting the two points of interest in the configurational space occurs following natural fluctuations. Additionally, the algorithm minimizes the sum of the potential and kinetic energy of the particles and of the bias potential. This provides means of controlling the quality and adiabaticity of the simulation.

In this study, we use the efficiency of the ABMD algortihm within a different perspective, i.e., to understand the complex dynamic processes of the early stage of droplet coalescence. The delicate issue is the choice of the damping constant in Eq. 5 that drives the mechanism under study. According to the classical theory of fluid particle coalescence Chan et al. (2011), we aim at tackling at a mesoscopic scale the mechanism associated to the drainage of the liquid films of the continuous phase between two approaching droplets. These liquid films should drain to allow the dispersed fluids to make contact and fuse. The coalescence mechanism strongly depends on the velocity of the two approaching droplets, which can be responsible for large deformation of the droplet shape FanSM2012 (). In the following, we consider ABMD parameters that allow us to study coalescence while keeping the shape of the droplets unaltered (i.e., the radius of gyration of the droplets). This means we shall use a value of small enough so as to provide the correct sequence of events associated to coalescence, but large enough to make this happen in a computationally reasonable time. The second ABMD parameter of choice is the target value , which is set on the distance between the two water droplet before coalescence occurs. is set equal to the radius of the isolated droplets, and determined by the positions of the NPs. These data are reported in Table. 3. Biased simulations were performed using the version 2.1 of the plugin for free-energy calculation, named PLUMED Bonomi et al. (2009); Tribello et al. (2014).

### ii.4 Pulling algorithm

Pulling dynamics, also named steered molecular dynamics (SMD) is an algorithm which takes inspiration from single-molecule pulling experiments Grubmuller et al. (1996) and forces a system to evolve away from its initial equilibrium condition. This is intended to accelerate transitions between different energy minima and to estimate the free-energy difference between two (or more) states via the Jarzynski equality Park et al. (2003). Pulling simulations have become very popular and are extensively applied in studying many biophysical processes Isralewitz et al. (2001); Lu et al. (1998) and phenomena in chemical physics FanSM2012 (). In SMD, the Hamiltonian of the system, , is modified into a new Hamiltonian, , which contains a term which depends on time via a harmonic spring potential centered on a point which moves linear with time, i.e.,

(8) |

where is the chosen coordinate of the configuration of the system that is biased, , the pulling velocity, and the spring constant. For an appropriate set of SMD parameters (, ), the system follows closely the center of the moving harmonic spring.

## Iii Results and discussion

In the following, we consider two types of numerical simulations, (1) Pulling simulations where the two droplets are moved towards each other using one harmonic spring with a force constant of tethered to their center of mass, at constant relative velocity (pulling rate) or per time step, and (2) ABMD simulations where we consider the distance between the center of mass of the two droplets as the biased coordinate , a damping constant and a target value .

In Fig.1 (left panel) we show simulation snapshots obtained during the two coalescence simulations: pulling simulations (left column) and ABMD simulations (right column). In both cases we consider water droplets stabilized by 160 NPs of type 55-JP in decane solvent. The interfacial area per NP equals .

Visual inspection of the simulation snapshots highlights the fundamental difference between the two processes. For the ABMD simulation, the shape of the water droplets remains unaltered, while in pulling simulation a clear alteration of the shape of the two coalescing droplets is observed as the distance between the center of mass of the droplets decreases. This behaviour is quantitatively shown in Fig.1 (right panel) where we show the temporal evolution of the distance between the centers of mass and the radius of gyration of the two coalescing droplets as a function of simulation time during ABMD and pulling simulations. Before the two droplets come into contact, the radius of gyration remains constant as the distance between the two droplets decreases. When the distance is such that the two droplets start interacting with each other via NP-NP interaction ( ), increases during the pulling simulation, meaning the droplets deform as the distance imposed decreases. keeps increasing until a minimal distance is reached ( ), when the two droplets coalesce.

If one considers the ABMD simulation, the distance between the two droplets reaches a plateau as soon as the two droplets come into contact, at around . This plateau is slightly lower than the value . This difference indeed can be explained by the fact that ABMD simulation drives the system such that the NPs surrounding the droplet interact with the solvent shell that surrounds the water molecules of the second droplet and thus avoid NP-NP interactions that are more energetically costly. The shape of the two droplets remains unaltered until the two droplets coalesce.

It is necessary to further quantify the difference between these two coalescence mechanisms at a mesoscopic scale. In Fig. 2 we show detailed simulation snapshots obtained during both pulling and ABMD simulations. Considering first qualitatively the pulling simulation, the two droplets start interacting when one NP in the incoming droplet (blue sphere in right panel in Fig. 2) comes into contact with one or more NPs from the opposite droplet (red spheres in Fig. 2). This NP-NP interaction initiates the shape alteration of the interacting droplets. We identify an interaction area as the set of NPs on the neighbour droplet (red spheres in Fig. 2) that defines the first neighbouring shell (i.e., distance between blue NP and red NPs less than ) around the incoming NP before merging. The threshold is set based on observations from prior work Luu et al. (2013). Then the incoming NP shifts towards the center of the interaction area delimited by the red NPs in Fig. 2. This is quantitatively described in Fig. 3 (left panel) where we follow the evolution of the minimal distance between the center of one incoming droplet and the water molecules belonging to the other droplet. The first part of the evolution is representative of the pulling rate, until the distance reaches a plateau around corresponding to where the NP/NP interactions start. Then the incoming NP starts interacting with the solvent in the interacting area and keeps compressing the droplet. This evolution goes on until reaches the minimal solvent shell, ie., , and the incoming NP starts coming into contact directly with the water molecules ( ). Merging of the two droplets occurs when two NPs (each belonging to different droplets) interact with water molecules belonging to the other droplet (cf. arrows (4) and (5) in Fig. 2). A water bridge is thus formed between the droplets, which allows the flow of water molecules.

Let us now compare the behaviour just described with the one observed running the ABMD simulations. Considering the simulation snapshots in Fig. 2 and the quantitative analysis in Fig. 3, we first see that the system does not spend a significant amount of time in the plateau at . Analysis of the trajectories suggest that the NP in the incoming droplet avoids frontal interactions with NPs belonging to the other droplet. Indeed, as prescribed by the ABMD algorithm, the system follows a minimal free-energy path such that the incoming droplet mainly interacts with the solvent molecules located around the interaction area delimited by the red subset of NPs in Fig. 2. The system reaches progressively the plateau corresponding to the minimal solvent shell, i.e., , and the incoming NP eventually comes into contact with the water molecules ( ). The final stages of the process are similar to what is observed during pulling simulations, where a set of two NPs must come into contact with the water molecules of the other droplet to initiate coalescence. Unlike the results observed during the pulling simulations, this evolution happens without an increase of the droplets radius of gyration.

This difference is due to the pulling parameters chosen, which allow the system to sample two different mechanisms during coalescence. The parameters chosen for the pulling simulations are such that the system remains in the advection regime. This can be highlighted with the use of the Péclet number Huysmans and Dassargues (2005)

(9) |

where is the effective advection velocity, a characteristic length, and the diffusion coefficient. Considering the pulling velocity in the center of mass reference frame, , the DPD characteristic length, , and the simulated water self-diffusion coefficient , one obtains . This value for the Péclet number indicates that during the pulling simulations the system is dominated by the pulling velocity. On the contrary, in the ABMD simulations, the evolution follows the natural fluctuations of the system, a regime dominated by diffusion by design.

Our results suggest that when pulling simulations are considered, the coalescence process depends on parameters such as NP contact angles and concentrations. On the contrary, the results obtained from ABMD simulations do not show such dependencies. We illustrate this difference with the example in Fig. 4 (left panel) where we keep the interfacial area per NP constant (160 NPs) and consider spherical 60JP NPs. In this figure the contact point (C) corresponds to the first NP-NP interaction between the two droplets and the merging point (M) identifies the formation of the water bridge between the two droplets. In Fig. 4 we illustrate that both C and M are different in either pulling or ABMD simulations due to the shape deformation of the system. The system seeks to keep the NP contact angles as close as possible to their equilibrium value (data consistent with Table 2). Our simulations suggest that under these conditions, the droplets deform their shape and increase the water-oil interfacial area rather than allowing the NPs to change their contact angle or desorb from the interface. The increased interfacial area, not covered by NPs as shown in Fig. 4, favours coalescence.

The differences between the results obtained between ABMD and pulling simulations also increase when the
surface density of the NPs on the droplets increases. We illustrate this difference with the example
in Fig. 4 (right panel), where we increase the interfacial area per NP (from 160 NPs to 165 NPs) and consider
spherical 60JP NPs. Under these conditions, the NP density becomes sufficiently high that the droplet
cannot create a large NP-free area (as we observed for 160 NPs) to facilitate coalescence.
When the distance between the two droplets decreases, one NP on the incoming droplet
comes in contact with the solvent shell and then the water molecules from the other droplet.
We find it interesting that in this scenario coalescence initializes in an area different from
the initial interaction area (yellow spheres in Fig. 4). This NP is different
from the NP that induces the compression of the neighbour droplet (red spheres in Fig. 4).
The contact point (C) and the merging point (M) are different and the merging does not require the contact of two NPs
originally on the two separable droplets.

In the ESI, we report the results obtained when homogeneous NPs were simulated.
The results are consistent with pulling simulations reported previously FanSM2012 (). The key point is
the difference in the interaction
between the NP on the incoming droplet, and the other droplet. Due to the homogeneous distribution of polar and apolar beads
on the NP surface, the homogeneous NP from the incoming droplet can adsorb simultaneously on the two interfaces, yielding the same
contact angle on both droplets. As a consequence, in both pulling and ABMD simulations, the last step of the coalescence mechanism
differs from the one observed with spherical Janus NPs. When homogeneous NPs are used, the coalescence requires only the contact of one NP to induce
the merging of the two water droplets.

We also considered the impact of different NP shape on the coalescence mechanisms, focusing our attention on ellipsoidal NP of aspect ration (242 beads around the ellipsoidal NP). In Fig. 5 we first consider the qualitative evolution of the system. As in the case with spherical NP, the ABMD simulation does not alter the shape of the droplet, while the pulling simulation induces an increase of the droplets gyration radius. This is visually investigated in Fig. 6 and more quantitatively shown in Fig. 7. The compression of the droplet observed during pulling simulations can be explained in a first instance by the same mechanism we previously described, i.e., the pronounced NP-NP interaction due to the non-diffusive behaviour of the system. The incoming NP interacts with the solvent shell around the water molecules and eventually comes into contact with the water molecules themselves. However, the last step of the coalescence mechanism differs from the one observed with spherical Janus NPs (cf. ESI for visual inspection). When the droplets are stabilized by ellipsoidal NPs, the coalescence process does not require the contact of two NPs (one on each droplet). Due to the longitudinal orientation of the ellipsoidal NP with respect to the interface, contact happens between the apolar NP beads (green spheres in Fig. 6) on the incoming droplet and the water molecules of the second droplet. When this contact is secured, water molecules from the second droplet can travel a short distance and interact with the polar NP beads (purple spheres in Fig. 6), originally exposed to the decane solvent. Once this interaction is established, a molecular bridge is formed between the two droplets that leads to coalescence. It is worth repeating that this coalescence mechanism is facilitated by the longitudinal orientation of the Janus ellipsoidal NPs on the droplet interface.

## Iv Conclusions

In this paper, we considered Adiabatic Biased Molecular Dynamics simulations to investigate the early stage of the coalescence process between two water droplets in decane solvent at the mesoscopic level. The droplets are stabilized with various NPs (Janus and homogeneous) of varying contact angles and shapes, and different NP densities. In ABMD simulations, the system is forced to coalesce adiabatically. This allows us to study how the ability of the NPs to stabilize Pickering emulsions is discriminated by NP mobility and the associated caging effect. We compared how the coalescence process differs when one considers pulling simulations with sufficiently high steering speeds. When ABMD simulations are considered, we showed that the coalescence mechanism does not depend on parameters such as NP contact angles and concentrations. However, it is sensitive to the NP type (Janus NP vs. Homogeneous NP) and shape (spherical vs. ellipsoidal). The coalescence occurs when a water bridge is formed between the two droplets, which drives the subsequent flow of water molecules. The formation of the water bridge depends on the shape of the NPs used as stabilizers. Considering spherical NPs, the water bridge is formed with two NPs (each belonging to different droplets) that interact with water molecules from the other droplet. This differs from ellipsoidal NPs of aspect ratio . Due to the longitudinal orientation of the ellipsoidal NP with respect to the interface, water molecules from the second droplet can interact more easily with the polar NP beads on the incoming droplet originally exposed to the decane solvent. They form a molecular bridge with a single ellipsoidal NP between the two droplets that leads to coalescence.

When pulling simulations are considered, we observed a wider range of coalescence processes. When the system remains in the advection regime, as defined by the Péclet number, the droplet has the ability to undergo large elastic deformations before coalescing. This is due in a first instance to the NP-NP interactions, which are very sensitive to the NP density on the droplets. This is related to the caging effect associated to NP mobility. The collision speed also influences the alteration of the shape of the water droplet, as the coalescence process is sensitive to the time scale associated to the liquid film drainage. These different parameters, associated with the fact that the system tries to keep the NP contact angles as close as possible to their equilibrium value when the shape of the system changes, are responsible for a wide range of coalescence processes.

These results are consistent with the recent microfluidic collision experiment of Zhou et al. Zhou et al. (2016)
studying head-to-head moving water droplets in a microfluidic chip.
Following preliminary investigations on droplet coalescence Liao and Lucas (2010); Leal (2004),
the authors analysed the coalescence process
using the droplet contact time and the liquid film drainage time (i.e., the time interval from droplet
contact to fusion). They defined a coalescence percentage to indicate how coalescence is prone to occur under different
flow conditions. However, as the droplets have the opportunity to detach from each other when droplet velocity is too high,
the coalescence percentage can eventually reach a zero value. The authors highlighted that droplets coalesce
effectively at low velocities and suggested that this is related to the importance of the liquid film drainage between
two approaching droplets and the droplet faculty to undergo elastic deformation.
Our numerical analysis allows us to investigate the early stage of the coalescence process
at a mesoscopic scale and to propose detailed coalescence mechanisms, the observation of which
remains at a scale inaccessible through experiments. Although
the present work only considers frontal droplet-droplet collisions, our results highlight the importance of droplet speed
in modulating the coalescence process.

To go further, it would be interesting to modify the chemistry of the NPs to study their ability to enhance or to perturb emulsion stability. Considering both ABMD and pulling simulations, we showed that the two water droplets start coming into contact through NP-NP interactions. Adjusting the apolar-apolar interaction between NPs might have an impact on the coalescence process, perhaps increasing or shortening the NP-NP interaction time and manipulating the effective interaction with the solvent shell. Similarly, it would be interesting to study the effect of particle diameters on the coalescence process. Given a sufficiently high NP density around the droplets, increasing the diameter of the Janus NPs would also increase the NP-NP interaction time without changing the NP chemistry.

## acknowledgement

We acknowledge X-C. Cuong, V. Garbin, and L. Botto for useful discussions. F.S. thanks M. Salvalaglio for fruitful discussion concerning the ABMD algorithm. Via our membership of the UK’s HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202), this work used the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk).

## References

- Pickering (1907) S. Pickering, J. Chem. Soc., 1907, 91, 2001–2021.
- Drexler et al. (2012) S. Drexler, J. Faria, M. Ruiz, J. Harwell and D. Resasco, Energy Fuels, 2012, 26, 2231–2241.
- Angelova et al. (2011) A. Angelova, B. Angelov, R. Mutafchieva, S. Lesieur and P. Couvreur, Acc. Chem. Res., 2011, 44, 147–156.
- Shchukina and Shchukin (2012) E. Shchukina and D. Shchukin, Curr. Opin. Colloid Interface Sci., 2012, 17, 281–289.
- Puglia and Bonina (2012) C. Puglia and F. Bonina, Expert Opin. Drug Delivery, 2012, 9, 429–441.
- Miller et al. (2006) R. Miller, V. Fainerman, V. Kovalchuk, D. Grigoriev, M. Leser and M. Michel, Adv. Colloid Interface Sci., 2006, 128–130.
- Wang et al. (2011) D. Wang, S. Yordanov, H. Paroor, A. Mukhopadhyav, C. Li, H. Butt and K. Koynov, Small, 2011, 7, 3502–3507.
- Zahn et al. (1997) K. Zahn, J. MendezAlcaraz and G. Maret, Phys. Rev. Lett., 1997, 79, 175–178.
- Cheung (2010) D. Cheung, Chem. Phys. Lett., 2010, 495, 55–59.
- Weeks and Weitz (2002) E. Weeks and D. Weitz, Chem. Phys., 2002, 284, 361–367.
- Sonnenburg et al. (1991) J. Sonnenburg, D. Kremp and R. Sandig, Mol. Phys., 1991, 74, 649–664.
- Peng et al. (2009) Y. Peng, W. Chen, T. Fischer, D. Weitz and P. Tong, J. Fluid. Mech., 2009, 618, 243–261.
- Cheng and Grest (2012) S. Cheng and G. Grest, J. Chem. Phys., 2012, 136, 214702–214709.
- X-C.Luu et al. (2013) X-C.Luu, J. Yu and A. Striolo, J. Phys. Chem. B, 2013, 117, 13922–13929.
- Fan and Striolo (2012) H. Fan and A. Striolo, Soft Matter, 2012, 8, 9533–9538.
- Fan and Striolo (2012) H. Fan and A. Striolo, Phys. Rev. E, 2012, 86, 05610–05120.
- Paci and Karplus (1999) E. Paci and M. Karplus, J. Mol. Biol., 1999, 288, 441–459.
- Camilloni et al. (2011) C. Camilloni, R. Broglia and G. Tiana, J. Chem. Phys., 2011, 134, 045105.
- Marchi and Ballone (1999) M. Marchi and P. Ballone, J. Chem. Phys., 1999, 110, 3697–3702.
- Groot and Warren (1997) R. Groot and P. Warren, J. Chem. Phys., 1997, 107, 4423–4435.
- Plimpton (1995) S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
- Luu et al. (2013) X.-C. Luu, J. Yu and A. Striolo, Langmuir, 2013, 29, 7221–7228.
- Groot and Rabone (2001) R. Groot and K. Rabone, Biophys. J., 2001, 81, 725.
- Calvaresi et al. (2009) M. Calvaresi, M. Dallavalle and F. Zerbetto, Small, 2009, 5, 2191–2198.
- Partington et al. (1952) J. Partington, R.F.Hudson and K. Bagnall, Nature, 1952, 169, 583.
- Luu and Striolo (2014) X.-C. Luu and A. Striolo, J. Phys. Chem. B, 2014, 118, 13737–13743.
- Mejia et al. (2012) A. Mejia, A. Diaz, S. Pullela, Y.-W. Chang, M. Simonetty, C. Carpenter, J. Batteas, M. Mannan, A. Clearfield and Z. Cheng, Soft Matter, 2012, 8, 10245–10253.
- Dendukuri et al. (2006) D. Dendukuri, T. Hatton and P. Doyle, Langmuir, 2006, 23, 4669–4674.
- Qi et al. (2014) F. Qi, J. Wu, G. Sun, F. Nan, T. Ngai and G. Ma, J. Mater. Chem. B, 2014, 2, 7605–7611.
- French et al. (2015) D. French, P. Taylor, J. Fowler and P. Clegg, J. Colloid Int. Science, 2015, 441, 30–38.
- Mason and Bibette (1996) T. Mason and J. Bibette, Phys. Rev. Lett., 1996, 77, 3481–3484.
- Delmas et al. (2011) T. Delmas, H. Piraux, A. Couffin, I. Texier, F. Vinet, P. Poulin and M. C. anf J. Bibette, Langmuir, 2011, 27, 1683–1692.
- Vymetal and Vondrasek (2011) J. Vymetal and J. Vondrasek, J. Phys. Chem. A, 2011, 115, 11455–11465.
- Solc (1971) K. Solc, J. Chem. Phys., 1971, 55, 335.
- Provasi and Filizola (2010) D. Provasi and M. Filizola, Biophys. J., 2010, 98, 2347–2355.
- Chan et al. (2011) D. Chan, E. Klaseboer and R. Manica, Soft Matter, 2011, 7, 2235–2264.
- Bonomi et al. (2009) M. Bonomi, D. Branduardi, G. Bussi, C. Camilloni, D. Provasi, P. Raiteri, D. Donadio, F. Marinelli, R. B. F. Pietrucci and M. Parrinello, Comp. Phys. Comm., 2009, 180, 1961.
- Tribello et al. (2014) G. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, Comp. Phys. Comm., 2014, 185, 604.
- Grubmuller et al. (1996) H. Grubmuller, B. Heymann and P. Tavan, Science, 1996, 271, 997–999.
- Park et al. (2003) S. Park, F. Khalili-Araghi, E. Tajkhorshid and K. Schulten, J. Chem. Phys., 2003, 119, 3559–3566.
- Isralewitz et al. (2001) B. Isralewitz, J. Baudry, J. Gullingsrud, D. Kosztin and K. Schulten, J. Mol. Graphics Modell., 2001, 19, 13–25.
- Lu et al. (1998) H. Lu, B. Isralewitz, A. Krammer, V. Vogel and K. Schulten, Biophys. J., 1998, 75, 662–671.
- Huysmans and Dassargues (2005) M. Huysmans and A. Dassargues, Hydrogeology J., 2005, 13, 895–904.
- Zhou et al. (2016) Q. Zhou, Y. Sun, S. Yi, K. Wang and G. Luo, Soft Matter, 2016, 12, 1674–1682.
- Liao and Lucas (2010) Y. Liao and D. Lucas, Chem. Eng. Sci., 2010, 65, 2851–2864.
- Leal (2004) L. Leal, Phys. Fluids, 2004, 16, 1833–1851.

Numerical analysis of Pickering emulsion stability: insights from ABMD simulations

Supplemental Information

## Appendix A Droplet characteristic

We consider nanoparticle (NP) types and concentrations that are expected to be (1) strongly effective at preventing droplet coalescence, and (2) maintain the shape of the droplet spherical. As described in prior work FanSM2012 (), the first criterion can be related to the interfacial area per NP, . The NPs are not expected to be strongly effective at preventing droplets coalescence when is larger than , independently on whether water or decane droplets are considered. Considering the values for reported in Table 2 in the main text, water droplets stabilized with NPs have an interfacial area per NP . These systems do not show high stability, as quantitatively shown in Fig. S1 where we follow the temporal evolution of the radius of gyration, , of one water droplet stabilized with 144 spherical NPs of type 55JP. The temporal evolutions of are shifted arbitrarily along the ordinate axis for clarity. In both pulling and ABMD simulations, remains constant as the distance between the two droplets decreases. The departure from the plateau indicates the merging of the two droplets. This merging happens without significant increase of , indicating that the resistance to coalescence provided by the NPs is minimal.

The NPs can affect the shape of the droplet at small interfacial area per NP, , because the system tries to keep the NP contact angles as close as possible to their equilibrium value (cf. evolution of the contact angle in Table 2 in the main text). In Fig. S2 we compare water droplets stabilized with either spherical NPs (left panel) or spherical NPs (right panel) of type 55JP. Visual inspection demonstrates the alteration of the droplet shape when NPs are used. In correspondence to the changes described in Fig. S2, the changes in relative shape anisotropy, reported in Table 2 in the main text, are of one order of magnitude of difference.

## Appendix B Water droplets stabilized by homogeneous NPs

We considered the ability of homogeneous (HP) NPs to stabilize Pickering emulsions. As we explained in the main text, due to the homogeneous distribution of polar and apolar beads on the NP surface, one homogeneous NP can adsorb simultaneously on two interfaces yielding the same contact angle on both droplets. The results obtained during the present coalescence simulations are consistent with pulling simulations reported previously FanSM2012 (). In both pulling and ABMD simulations the last step of the coalescence mechanism differs from the one observed with spherical Janus NPs. In Fig. S3 we show detailed sequence of simulation snapshots representing typical collision processes of two water droplets in decane solvent. The droplets are stabilized by 160 spherical NPs of type 75HP (contact angle ). The results are obtained with ABMD simulations. The key point is the interaction between the incoming NP and the second droplet. Indeed, the coalescence process requires only the contact of one NP to induce the merging of the two water droplets.

## Appendix C Coalescence mechanism for ellipsoidal NPs

We show in Fig. S4 detailed sequences of simulation snapshots representing typical coalescence processes of two water droplets stabilized with either 160 spherical NPs of type 60JP (left panel) or 160 ellipsoidal NPs of type 70JP (right panel). The results are obtained using ABMD simulations. When the droplets are stabilized by ellipsoidal NPs, the coalescence process does not require the contact of two NPs from the two merging droplets. This is due to the longitudinal orientation of the ellipsoidal NP with respect to the interface, as explained in the main text.

## References

- (47) H. Fan and A. Striolo, Mechanistic study of droplets coalescence in Pickering emulsions. Soft Matter 2012, 8, 9533-9538.