An Algorithm for Emulsion Stability Simulations: Account of Flocculation, Coalescence, Surfactant Adsorption and the Process of Ostwald Ripening

The first algorithm for Emulsion Stability Simulations (ESS) was presented at the V Conferencia Iberoamericana sobre Equilibrio de Fases y Diseño de Procesos [Luis, J.; García-Sucre, M.; Urbina-Villalba, G. Brownian Dynamics Simulation of Emulsion Stability In: Equifase 99. Libro de Actas, 1st Ed., Tojo J., Arce, A., Eds.; Solucion’s: Vigo, Spain, 1999; Volume 2, pp. 364–369]. The former version of the program consisted on a minor modification of the Brownian Dynamics algorithm to account for the coalescence of drops. The present version of the program contains elaborate routines for time-dependent surfactant adsorption, average diffusion constants, and Ostwald ripening.


Introduction
Emulsions are dispersions of two immiscible liquids, kinetically stabilized by the action of a surface-active substance known as a surfactant. The role of the surfactant in emulsion stability is decisive. Even the external phase of the emulsion resulting from the mixing of oil (O) and water (W) depends on the surfactant solubility [1][2][3][4][5]. In the absence of surfactants, drops aggregate rapidly as a consequence of the van der Waals force. Surfactants adsorb to the interface of the drops creating a repulsive barrier that decelerate aggregation and Ostwald ripening. They also favor the occurrence of immobile interfaces, delaying the drainage of the intervening film between flocculated drops [6].

OPEN ACCESS
Conversely, surfactants lower the interfacial tension of these O/W/O films, favoring the appearance of surface oscillations and holes. Depending on their interfacial properties, these films can drain and rupture or remain stable for long periods of time [7][8][9][10].
In order to understand the role of the surfactant in the stability of emulsions, simulations can be a valuable tool. However, there are severe computational restrictions for the simulation of emulsions apart from the large number of processes involved. First, the number of surfactant molecules in a typical system is very high even at dilute concentrations. Hence, it is not possible to simulate the movement of surfactant molecules explicitly along with the movement of the drops. Second, the time step of the simulation has to be very small in order to sample appropriately the potential of interaction between the drops. Third, drops of small size exhibit Brownian motion due to their thermal interaction with the solvent [11]. This has to be taken into account into their equation of motion. Fourth, it is not possible to simulate the behavior of drops using a constant potential as it is done in molecular dynamics. The potential of interaction changes with time as a consequence of surfactant adsorption and the progressive decrease of interfacial area of the emulsion due to the coalescence of drops. This demands the frequent calculation of the interfacial properties that determine the forces between the particles.
In this review we concentrate on the evolution of oil-in-water (O/W) emulsions composed of nondeformable drops. The next section introduces the technique of Brownian Dynamics simulations and the most relevant aspects of the classical theory of Derjaguin-Landau-Verwey-Overbeek (DLVO). Section 3 describes the algorithm of Emulsion Stability Simulations (ESS) in detail. In section 4 some illustrative results of ESS are discussed. Following, the modifications of the former algorithm required for the simulation of deformable drops are outlined. The paper finishes with a brief conclusion and the bibliography.

Brownian Dynamics (BD) Simulations
The movement of one small particle (1 nm -10 μm) diffusing in a quiescent media is the result of external and random forces. Random forces represent the effect of millions of collisions occurring between the solvent molecules and the particle surface. In the absence of external forces the energy for the displacement of a spherical particle is provided by the solvent, which at the same time takes away some energy from the particle in the form of friction. The whole process is a manifestation of the Fluctuation-Dissipation theorem that is expressed concisely in the form of the Stokes-Einstein equation [11][12]: (1) In this equation, D and f stand for the diffusion tensor and the resistance tensor of the particle (or a set of particles), and k B is the Boltzmann constant. In the case of a spherical particle the diffusion tensor is diagonal and its three components are equal (D). According to Equation (1), the absolute temperature of the reservoir (T ) is kept constant while the movement of the particle occurs, although there is a continuous exchange of energy between the particle and the solvent. This stochastic process can be described analytically in terms of the Langevin equation [13]. The solution of Langevin's equation for one spherical particle provides explicit expressions for its mean square displacement and its diffusion constant (D). The diffusion constant turns out to be proportional to the temporal correlation of the random fluctuations of the particle movement. The overall effect of the interaction between a suspended particle and the surrounding molecules is a random drift in the trajectory of the particle, whose displacement shows a Gaussian distribution with zero mean and a standard deviation equal to t D Δ 2 along each co-ordinate axis. Here t Δ is the lapse of time considered and D is the diffusion constant of a solitary particle [13]. The diffusion tensor can be evaluated from the resistance tensor: for a particle moving through the liquid as a consequence of external, hydrodynamic or interparticle forces. In turn, the resistance tensor can be obtained from the drag force experienced by the particle when it moves through the fluid at a constant velocity v [14]: The specific form of the diffusion constant of a sphere depends on the velocity of the liquid at the particle surface. If the surface is rigid and smooth, the velocity of the fluid becomes zero at the surface. In this case the diffusion constant comes out to be R T k D D B η π 6 0 = = , where: η and R are the shear viscosity of the solvent and the radius of the particle, respectively.
In the case of liquid drops, porous spheres, and bitumen drops, the diffusion constant can be expressed as: is a correction function depending on the characteristics of the O/W interface [14].
The addition of more particles to the system increases the complexity of the problem considerably. First, the random movement of the particles must be connected in such a way that they fulfil the Stokes-Einstein relation. Second, the movement of each particle generates fluxes (disturbances) in the solvent which affect the movement of the surrounding particles and its own. Thus, it is necessary to account for hydrodynamic interactions between the particles. Third, the particles interact with deterministic forces other than external forces. Thus, their movement is a combination of deterministic hydrodynamic and random forces.
One of the most widely used algorithms for Brownian Dynamics simulations is the one of Ermak and McCammon [15]. In this formalism the position of a particle at time t+ Δ t, ( ) t t r i Δ + is equal to: where the superscript "0" indicates that the variable is evaluated at the beginning of the time step. Subscripts i and j run over the particle coordinates ( ). F j is the sum of interparticle and external forces acting on direction j. D ij are the components of the diffusion tensor. The gradient of the diffusion tensor (second term on the right hand side) can usually made equal to zero selecting a tensor that depends on the distances between the particles and not on the particles coordinates. The third term on the right hand side stand for deterministic contributions and the fourth term correspond to the random displacements. A random deviate has the form: Here X j stands for a random variable sampled from a Gaussian deviate generator: , where ij δ is the Kronecker delta. The weighting factors are given by: According to the Ermak and McCammon equation (4) is compatible with a Fokker Planck description of the problem in the phase space. In the past, our group carried out the implementation of the above-mentioned algorithm with different tensors including Oseen, Rotne-Prager, and Batchelor's [16] finding several limitations of this technique for emulsion stability calculations. Equations (5)- (7) describe a particular method for taking the square root of the diffusion tensor ( t D Δ 2 ). However, the methodology suggested (Cholesky decomposition [17]) only works with particles of equal radii, and dilute systems. Other methods found in the bibliography like the QR decomposition also fail in simulations of polydispersed concentrated systems. This failure is caused by the assumption of pair wise additive hydrodynamic interactions. These schemes overestimate the hydrodynamic fluxes generated between the one central particle and its neighbours. They do not take into consideration that the flux coming from a particle located in the second coordination layer is partially screened by the inner neighbours. Dickinson suggested that the assumption of pair wise hydrodynamic interactions might even lead to negative diffusion constants [18]. When an average diffusion constant is used instead of a tensor, its value is de-coupled from the random deviates. Using this approximation and 0 D D = , we quantified the kinetic energy of micronsize particles produced by the random fluctuations of the Box-Muller algorithm [17]. The kinetic energy was computed from the time step of the simulation and the displacement of each particle in the absence of deterministic and external forces. It was observed that fluctuations as high as 10,000 k B T must be allowed in order to reproduce a mean square displacement equal to: . Cut off thresholds in the value of the deviates corresponding to kinetic energies of 1000 k B T and 100 k B T fail to reproduce the correct value of 2 r . It was evident that values of the displacement corresponding to the outskirts of the Gaussian distribution are necessary in order to simulate the Brownian movement of the particles correctly [19]. The above considerations are profoundly related to the outcome of ESS calculations and its discussion is not a mere technicality. The most famous theory of colloidal stability, Derjaguin-Landau-Verwey-Overbeek's DLVO theory [20] is based on the problem of diffusive passage over a potential barrier.
The presence of a strong repulsive force between the particles generates a potential barrier and two minima at each side of the barrier (Figure 1).

Figure 1.
Classical DLVO potential. The curve corresponds to the interaction potential between two 3.9-μm drops of dodecane (A H = 5.03 x 10 -21 J) suspended in water. The drops are partially covered with sodium dodecyl sulfate (C s = 10 -4 M, z s = -0.2057) using a homogeneous surfactant distribution (Section 3). On the one hand, primary minimum flocculation occurs at very small distances of separation and is assumed to be irreversible. Irreversible meaning that the aggregates formed do not separate if one lowers the ionic strength of the solution, increasing the repulsive force. This behaviour is due to the strong van der Waals force experienced by the particles at short distances of separation. On the other hand, secondary-minimum flocculation is usually reversible for small particles and could be "irreversible" for micron size drops [21]. Irreversible meaning in this case that the minimum is deep enough to prevent the particles from moving in and out of the potential well.
The occurrence of primary minimum flocculation depends on the diffusive passage of the particles over the potential barrier. According to Chandrasekhar [22] and Kramer [23] the probability of "jumping" over the barrier decreases exponentially with the size of the barrier. In order to account for the effect of the barrier on primary minimum flocculation, Fuchs [24] defined a "Stability ratio" W. The complete formula of W was deduced later by McGown and Parfitt [25]: Here f k stands for the flocculation rate. It can be either fast ( fast f k ) or slow ( slow f k ) if the potential of interaction between the particles is only due to attractive forces (V A ) or caused by a combination of attractive and repulsive (V R ) contributions (V T = V A + V R ). In Equation (8) the dependence of the friction on the distance between the particles (r) has been remarked. Numerical evaluations of Equation (8) lead Prieve and Ruckenstein [26] to a very useful relation between the height of the repulsive barrier (ΔV) and the stability ratio: Fuchs demonstrated that the theory of Smoluchowski [27] for fast flocculation could be conveniently modified in order to incorporate the average effect of a repulsive barrier through the stability ratio. Thus, the change in the total number of aggregates of a dispersion per unit volume (n), can be written as: where n 0 = n(t = 0). According to Equations (8)-(10) a repulsive barrier produces large values of W which retard the attainment of primary minimum flocculation. A measure of colloid stability towards flocculation is given by the half lifetime of the dispersion t 1/2 . This is the time required for a decrease of n 0 /2 in the initial number of aggregates: In the case of non-deformable drops, primary minimum flocculation implies coalescence. Drops that jump over the repulsive barrier necessarily coalesce since there is no other repulsive force to prevent it. Due to the irreversible nature of the coalescence process the correct simulation of the diffusion tensor and the random deviates is critical.

Equation of Motion
Emulsion Stability Simulations start from a cubic box that contains N drops randomly distributed. The particles move with an equation of motion similar to the one of Brownian Dynamic Simulations [28][29][30]: In Equation (12) subscripts i and j refer to particles i and j. The displacement of particle i during , is the result of two contributions. The second term on the right hand side of Equation (12) represents the effect of deterministic forces acting on particle i. It is composed of inter In the case of non-deformable particles is calculated following the methodology described in Ref. [30]. At every step during the simulation the space around particle "i" is divided in three spherical regions. If at least one neighbour particle reaches the internal region:  (3)). The first correction term ( ) 1 ( corr f ) takes into account those factors that change the expression of the diffusion constant of a particle at infinite dilution (D) [14]. And: In Equation (13) φ stands for the local volume fraction of oil around a central particle "i", and where i R is the radius of particle i, and 0 R a radius of a reference. The third term on the right hand side of Equation (12) gives the random deviates of the moving particles. The stochastic vector auss G r stands for a set of random numbers generated by the Box-Muller algorithm. The characteristic mean square displacement of the Brownian movement is obtained multiplying each random deviate by In ESS, non-deformable drops coalesce if the distance between their centres of mass gets smaller than the sum of their radii: When this occurs a new drop is created at the centre of mass of the coalescing drops. The radius of the new drop results from the conservation of volume:

Interaction Forces and the surface excess
In ESS the attractive force between the drops is determined by the effective Hamaker constant of two oil drops separated by water [32][33]. This constant (A H ) is often of the order of 10 -21 J for hydrocarbons and lattices. In the case of Bitumen some old experimental evidence indicated a value of A H ~ 10 -19 J, but recent evaluations suggest a much lower value (A H ~ 10 -21 J).
The van der Waals potential for two spherical drops of different radius (R 1 , R 2 ) is equal to [32]: There is always an attractive force between two drops of similar composition suspended in aqueous media. However, the repulsive force depends on the characteristics of the O/W interface and the type of surfactant adsorbed. Oil drops often exhibit an electrostatic surface potential between -10 mV (octacosane [34]) and -35 mV (benzene [35]) when suspended in water. This charge originates from the preferential adsorption of OHions to the oil/water interface. If this is the case, an initial value of the surface charge per unit area (σ 0 ) can be introduced as input in order to reproduce the total value of the surface charge σ T in the absence of surfactants.
When ionic surfactants are present, they add their charge to the initial surface charge of the drop (σ 0 ): Here N i , A i , z s , e stand for the number of surfactant molecules attached to drop i, the area of the drop: , the effective valence of a surfactant, and the unit of electrostatic charge (1.6 x 10 -19 Coul.). Γ stands for number of surfactant molecules per unit area at the oil/water interface. The interfacial area of one surfactant molecule at maximum packing (A s = 1 max − Γ ) is typically of the order of 50 Å 2 [36]. It depends on the adsorption time and on the way the surfactant partition between the immiscible phases and the interface. The effective charge of an ionic surfactant molecule (q = z s e) can be calculated from the zeta potential (ζ) of a drop saturated with surfactant molecules. In the most usual situation z s is varied until the experimental value of ζ is reproduced. This calculation is done assuming that σ 0 = 0. In this case, the effective charge of the surfactant contains the contribution coming from the hydration layer. It is also possible to use a finite value of σ 0 to reproduce the surface potential of an oil drop in water, and then vary z s to fit the value of ζ. The value of z s results from apportioning the effective charge of a drop to a discrete number of surfactant molecules (0.17 ≤ z s ≤ 0.23 [37]).
The estimation of q requires knowledge of the maximum number of surfactant molecules in each drop. This number is equal to: The current version of the program contains four analytical expressions for the calculation of the electrostatic force [38][39][40][41][42]. The formalism of Sader et al. [39][40][41] was used in previous simulations of non-deformable drops. According to these authors, the electrostatic potential between two spherical drops is given by Equation (21): where: In Equations (21)-(25), 0 ε is the permittivity of vacumm, ε the dielectric constant, κ -1 the Debye length, and T k e B P 0 the reduced electrostatic potential of the particle at its surface. The relation between the surface charge and the surface potential is given by: where: Notice that knowledge of σ T and the ionic strength of the solution allow the numerical evaluation of the surface potential from Equation (26), as well as the rest of the variables (Equations (22)-(25)) of the potential (Equation (21)). Since the ionic strength is set by the experimental conditions, and σ T depends on the surface excess, only Г is required in order to calculate the electrostatic potential between suspended drops.
Surfactants can be ionic, non-ionic, zwitterionic, and have a complex molecular structure. Hence, the interaction potential between the drops can vary amply. The program of ESS includes expressions for van der Waals, electrostatic, oscillatory, depletion, hydration, and steric forces.
In the case of non-ionic surfactants the stabilization force between emulsion drops is assumed to be steric [43][44]. Steric forces are composed of an osmotic (mixing) contribution and an elastic one: where h stands for the minimum distance between the surfaces of the drops. The mixing contribution is generated by the cross linking of the hydrophilic chains of the surfactants of each drop. The elastic force comes from the drastic modification of the surfactant conformation at close distances of separation between the particles. We developed expressions for both types of contributions modifying the equations reported by Vincent [45] and Bagchi [46]. On the one hand, Vincent used simple analytical formulas to describe the distribution of polymer monomers perpendicular to a planar interface. Following he used the Derjaguin approximation [47] to evaluate the expression of the free energy of mixing deduced by Flory and Krigbaum [48], for two spheres stabilized with polymer molecules. On the other hand, Bagchi calculated the exact volume of overlap between the polymer layers of two colliding drops, but used the Flory-Huggins [49] expression for the calculation of the free energy. Lozsán et al. [43] showed that the methodology of Vincent with the exact calculation of the volume of overlap is equivalent to the methodology of Bagchi if the free energy is evaluated with the expression of Flory-Krigbaum. The ESS code has several expressions for the calculation of the steric interactions [43] including the expression of De Gennes [50]. An example of these expressions is given by Equation (29) for distances between δ < h < 2δ [43][44]: In Equation (29) V w stands for the molar volume of the continuous phase, χ is the Flory-Huggins solvency parameter, and i φ stands for the volume fraction of surfactant molecules in the interfacial layer of drop i. In order to estimate i φ it is usually assumed that the hydrophobic chain of the surfactant is dissolved in the oil phase, and only its hydrophilic chain lies in the outer region of the drop. If the density of hydrophilic chains in the interfacial layer is assumed to be constant, i φ can be expressed in terms of the surface excess [44]: where: ρ 2 stand for the density of the hydrophilic chain.
As in the case of the electrostatic potential, the steric potential depends on some parameters and is a function of the surface excess.

Surfactant Distribution (Evaluation of the surface excess)
The value of the surface excess of a surfactant at the interface of emulsion drops cannot be measured directly. It is extrapolated from the variation of the interfacial tension in systems with a macroscopic O/W interface. Unfortunately, emulsions are constantly evolving, continuously changing their drop size distribution and total interfacial area, A T : As a result, the surface excess is not constant and the interaction potential between the drops changes as a function of time.
Ionic surfactants are not soluble in the oil phase. They can migrate to the oil phase in the form of inverse micelles if the salt concentration of the system is unusually high. In the typical case, ionic surfactants adsorb to the O/W interface before they form micelles. The amount of surfactant required for the complete coverage of the interface can be larger or equal to the critical micelle concentration (CMC). It varies with the volume fraction of oil (φ) and A T . The larger the values of φ and A T , the higher the surfactant concentration required for the complete coverage of the drops.
In the case of non-ionic surfactants the situation is more unpredictable. The partition of non-ionic molecules depends on the affinity of its lipophilic and hydrophilic moieties for the oil and water phases. For example, the partition coefficient of alkylphenol oligomers between water and n-alkanes is equal to [51]: where m is the number of ethylene oxide units in the surfactant molecule, SACN is the number of carbon atoms in the alkyl chain of the surfactant, and ACN is the number of carbon atoms in the nalkane molecule. The adsorption of surfactant molecules depends on time and on several formulation variables like the number of methylene groups of the oil molecule (Equation (32)), the salt concentration, the presence of alcohol molecules in the system, etc. Furthermore, adsorption can be reversible or irreversible [52][53].
The routines of surfactant distribution attempt to recreate the most common experimental situations. The strategy of ESS is to apportion the surfactants to the interfaces of the drops in such a way that it reproduces the variation of the surface excess in the experimental system. Consequently, only the movement of the drops is considered explicitly in Equation (12).
Some of the routines for surfactant distribution have a formal theoretical background [54][55]. Some others are only practical approximations to the very complex problem of surfactant diffusion and adsorption [56][57]. These routines resemble the cases of Homogeneous and Non-homogeneous surfactant distributions. The effect of non-homogeneous distributions results from an incomplete mixing of the emulsion components. Most recent studies concern the cases in which the surfactant is evenly distributed among the surfaces of the drops. Three examples of these strategies are given below.

Homogeneous Surfactant Distribution with fast and irreversible surfactant adsorption
This is the simplest methodology. It consists in ascribing to each drop a number of surfactant molecules ( i s N , ) proportional to its interfacial area.
Here T s N , stands for the total number of surfactant molecules in the system. Equation (33)  If the number of drops decreases as a consequence of coalescence, Equation (33) can be used to recalculate the number of surfactant molecules adsorbed to the remaining drops. The number of surfactant molecules of each drop increases as the calculation progresses because the total interfacial area of the emulsion decreases as the drops coalesce. However, the interfacial area of one surfactant molecule cannot be lower than its minimum cross-sectional area at the O/W interface ( s A ). Hence, . If the number of surfactants in the system surpass the amount required for the complete coverage of all the drops, the value of i s N , is set equal to max ,i s N .

Homogeneous Surfactant Distribution with time-dependent surfactant adsorption
In the case that the adsorption is basically controlled by the diffusion of surfactant molecules from the bulk to the subsurface [58], the surface excess is a function of the total surfactant concentration C T , the diffusion constant of the surfactant D s , and the time: The mechanism of surfactant adsorption can be very involved, including barriers of adsorption, reorientation at the surface, etc. However, Liggieri et al. [59] demonstrated that most process of mixed adsorption kinetics can be reformulated in terms of a diffusion controlled mechanism (Equation (34)) if D s is substituted by an "apparent" diffusion constant (D app ). Furthermore, this equation is also compatible with the findings of Hua and Rosen [60][61] for a large number of surfactants with different molecular structures. According to these authors, the surface tension of most surfactants shows four characteristic regions of change known as: the induction period, the fast-fall region, the mesoequilibrium region and the equilibrium region. The time required for the surface pressure to drop to half its value at mesoequilibrium, is found to follow Equation (35). This equation allows estimating the value of D app required for the evaluation of Equation (34).
The routine of time dependent surfactant adsorption uses D app , and C T as input, estimating the value of the surface excess from ( ) Here ( ) t A s stands for an "effective" area per surfactant molecule at the oil/water interface. As ( ) t Γ increases, ( ) t A s decreases, increasing the number of surfactant molecules adsorbed to each drop.

Equilibrium Surfactant Distribution (Gibbs)
Equilibrium isotherms are only attained after long periods of time. In this routine we assume that: a) the adsorption of surfactant is very fast, and b) equilibrium adsorption is obtained instantaneously.
Gibbs isotherm only applies to that range of surfactant concentration between the beginning of the decrease of the interfacial tension C T = C 1 , and the Critical Micelle Concentration (CMC). When C T ≥ C cmc the maximum number of surfactants per drop is already adsorbed: (20)).
When C T ≤ C 1 the change in the interfacial tension often shows a small maximum. We disregard this feature of the experimental curve. Instead we approximate the low concentration limit either by a) using an additional Gibbs isotherm between C T = 0 and C T = C 1 ; or b) calculating an average surface excess including the point C T = 0, γ = γ 0 .
In terms of finite differentials the equation of Gibbs is equal to: where c = 1 or 2 depending on the type of surfactant and the ionic strength of the solution. Assuming proportionality between the suggested value of the tension and the number of surfactant molecules adsorbed: Equation (38) is very convenient because it expresses the tension in terms of the surfactant population of each drop. According to Equations (37)- (38): Equation (39) gives the number of surfactants attached to drop i as a function of the total surfactant concentration in the system. Here, 1 γ stands for the value of the interfacial tension at C 1 . The value of Notice that there might be cases in which the number of surfactant molecules in the system might not be enough to cover the drops according to Equation (39). In those cases a homogeneous surfactant distribution is applied (Equation (33)) despite the initial selection of the Gibbs routine.

Ostwald Ripening
Although a large number of numerical techniques are available for the simulation of Ostwald ripening [62][63][64][65][66] they are based on population balance equations and cannot be incorporated into the algorithm of ESS. In order to simulate the process of Ostwald ripening we included the algorithm of De Smet et al. [67][68][69] in our ESS code. The fundamental equation of this method is derived from the Fick's law and the Kelvin's equation assuming α << R i : interface. R c is the critical radius of the emulsion equal to [70]: Here N is the total number of drops. Parameter α is the so called capillary length, defined as: In Equation (42) V m is the molar volume of the oil and R the gas constant. Defining: (43) and: A simple equation for the exchange of oil molecules is obtained: At any step of the simulation there exists a critical radius R c of the emulsion. Particles with R i < R c , dissolve while particles with R i > R c grow. Particles with the same radius as the critical radius R i = R c preserve their size. The number of molecules exchanged by particle i is equal to the product ( ) ( ) In ESS the value of M(t) is set once the time step of the simulation is chosen, ( ) M t M = . When the smallest particle i = small contains fewer molecules than the number it should lose according to Equation (45), Recently we implemented a new procedure to avoid a substantial decrease in the number of particles [71]. The simulation starts from a given Drop Size Distribution (DSD) of N 0 drops, and evolves until it reaches N(t = t') = N 1 = 200. At this point, a new Drop Size Distribution (DSD) with N = N 0 drops is built. This can be achieved approximating the probability distribution of particle sizes by the relative number of drops of each size existing at time t = t': is the number of drops with radius R i existing at time t'. Thus, the number of drops of , is equal to: Use of Equations (47)-(48) produces a new DSD with N 0 drops that exactly matches the old one. The auxiliary code also calculates the new size of the simulation box, which is required in order to preserve the initial volume fraction of oil with the new set of particles. Once the input file is modified and the new distribution read from an external file, the program generates a new set of co-ordinates for the new particles. As in the beginning of the simulation the particles are distributed at random avoiding overlap. At this point the code resumes the calculation of the main cycle (see below).

The Algorithm of ESS
The algorithm for ESS is shown in Figure 2 in the form of a flowchart.  At the beginning of the simulation the code reads or generates a drop size distribution and a set of co-ordinates for each particle. Additionally, the time step(s) of the simulation (single or double) must be specified. A combination of a small time step and a large one is used for the calculation of dilute dispersions [72]. For this purpose a distance of closest approach (d min ) must also be selected. The maximum range of the interaction potential usually approximates this distance (50 nm ≤ d min ≤ 100 nm). A double time-step calculation implies the use of a longer time step (Δt L ) when the particles are far apart from each other, r ij > d min . If the distance between two particles becomes equal or lower than d min , all particles are returned to their previous positions, and a small time step (Δt S ) is used in . This technique is efficient in those cases where the repulsive potential does not prevent coalescence. Otherwise the drops aggregate making r ij always lower than d min . This causes the continuous use of Δt S after the first floc is formed. After the creation of the DSD, the program distributes the surfactant molecules among the drops. There are 11 routines to cover the most common experimental situations.
At this point the interfacial tension of each drop can be determined. The calculation of the diffusion constants is also executed here, because some of its expressions depend on interfacial parameters. First, the program assigns the diffusion constant of Stokes to all particles according to their radius. Second it makes the corrections necessary to account for the effect of the interfacial properties and hydrodynamic interactions on the diffusion constant ( ) ).
Following, the forces are calculated and the drops are moved according to Equation (12). If the Ostwald ripening mode is selected, the oil molecules are exchanged at this point, and the program follows the path indicated by the thin solid arrows in Figure 2. Otherwise the program follows the path indicated by the thin broken-line arrows. In the latter case, the program checks for coalescence after the particles are moved. In the former case the program checks for coalescence after the exchange of oil molecules takes place.
If coalescence occurs, the total interfacial area of the emulsion changes along with the volume of one (or several) drop(s). Therefore the surfactant population had to be redistributed among the drops. The same occurs if Ostwald ripening happens. However, in the latter case a minimum number of drops must be maintained. Therefore, the program checks if the remaining number of drops is higher than a fixed value (N(t = t') > N 1 ). If this is the case, the program proceeds with the redistribution of surfactant molecules. Otherwise, the program stops. An auxiliary program is used to build a new drop size distribution with N = N 0 = N (t = 0) drops. This program also calculates the new size of the box necessary to keep the volume fraction of oil constant with the new set of particles. The new distribution is read from a file by the program of ESS. The co-ordinates of the new particles are then generated and the main cycle of ESS continues.
If the Ostwald ripening process is not selected, the program neither exchange oil molecules nor does it re-builds the drop size distribution when N(t = t') ≤ N 1 . In this case the code follows different routes depending on the outcome of the coalescence check. If coalescence occurs it is necessary to reassign the surfactant population. If it does not occur there are two possibilities. If one of the routines of time-dependent adsorption is used, it must recalculate the surfactant population anyway. If this is not the case, the program proceeds with the evaluation of the hydrodynamic interactions and a new cycle begins.
It is important to remark at this point that the change of the number of particles during the simulations is only equal to the variation of the number of aggregates in the absence of a repulsive force. This was demonstrated quantitatively in Ref. [73].
Whenever a repulsive barrier is present, the change in the number of aggregates as a function of time has to be calculated at the end of the simulation using a different program [21,[73][74]. This additional code uses as input the positions of the particles produced during the simulation and a fixed flocculation distance. The flocculation distance is usually approximated by the position of the secondary minimum of the interaction potential. When these computations are finished, statistical data regarding aggregates, flocs, radius of gyration, etc, is obtained.

Results
The general results of the simulations can be classified into three categories depending on the total surfactant concentration of the system [74].

Low surfactant concentration
If the surfactant concentration is not enough to stabilize the initial drop size distribution, drops coalesce as soon as they collide with each other. Hence, the change in the number of particles is equal to the change in the number of aggregates [73]. In this case, the systems follow the dynamics of Smoluchowski, in the sense that the number of aggregates changes as predicted by Equation (10)  Notice that Equation (49) was intended to explain the phenomenon of irreversible flocculation of solid particles. Danov et al. [75] demonstrated that Equation (49) also applies to emulsions subject to the simultaneous processes of flocculation and coalescence. For this demonstration no distinction was made between aggregates formed by the collision of smaller flocs, and those of the same size resulting from the partial coalescence of flocculated drops. Similarly, no distinction was made between single particles formed through coalescence, and those initially present at t = 0. These are the same assumptions we use for determining the variation of the number of the aggregates as a function of time. Danov et al. also supposed a constant collision kernel, but did not make any specific hypothesis regarding the coalescence rates.
The fact that Equation (49) holds for coalescing emulsions is related to the suppositions followed by Smoluchowski for determining the flocculation rate. He pictured the case in which one particle was fixed in space while the others collided with it as a consequence of a gradient of concentration. This gradient was established as soon as t > 0 between the collision radius of the fixed particle and the bulk of the liquid. However, the density of particles at the collision radius was assumed to be equal to zero during the whole aggregation process. Hence, the fixed particle acted as a perfect sink. Every particle that collided with the fixed particle was assumed to disappear at the moment of the collision. Moreover, the collision efficiency of a cluster composed of single particles was estimated using the same procedure, but employing the average hydrodynamic radius of the aggregate to approximate its collision radius. As a result of these assumptions, the process of flocculation described by the theory of Smoluchowski "includes" the possible occurrence of instantaneous coalescence. Instantaneous meaning that the time required for coalescence is negligible in comparison to the time required for flocculation.
Inhomogeneous surfactant distributions as well as reversible adsorption, also lead to fast aggregation and to the coalescence of drops, favouring the validity of Equation (49) [54]. The same occurs if the surfactant does not adsorb rapidly to the O/W interface (D app < 10 -12 m 2 /s) [55].  In general, the initial interfacial area of an emulsion is higher than the one that can be stabilised with the surfactant concentration available. In this case a Smoluchowskian drop in the number of particles occurs at short times. In the preparation of emulsions, chemical methods usually employ a larger surfactant concentration than mechanical ones. This favours the formation of a repulsive barrier as soon as the drops are formed. In this case, the dynamic of aggregation is not expected to follow Equation (49).

Intermediate surfactant concentration. Insufficient surfactant molecules to prevent the initial coalescence of drops
According to ESS of non-deformable droplets, the variation of the number of aggregates in the presence of an appreciable repulsive barrier conforms to a remarkably simple expression [73]: where A, B, k 1 and k 2 are constants and B = 1 -A. The second term in Equation (52) results from the consideration of the coalescence rate as a first order process, depending on the number of flocculated doublets. We used a Smoluchowskian term to represent flocculation and an exponential term to represent coalescence. Hence, k 1 and k 2 were formerly ascribed to the flocculation rate ( slow f k ) and the coalescence rate (k c ), respectively [73]. Coefficients A and B measure the extent of these two processes during the simulation. The curve described by Equation (52) generally shows a pronounced initial decrease followed by a much slower change in the number of aggregates per unit volume ( V N n a a = ). See Figure 4. The initial decrease corresponds to the combined processes of flocculation and coalescence (FC period). During the FC interval, the total interfacial area of the emulsion decreases. This causes a substantial redistribution of surfactant molecules amongst the remaining drops if a homogeneous surfactant distribution is assumed. As a result there is a progressive increase in the repulsive potential between the drops, which slows down both destabilisation processes. Equation (52)  distributions. This is remarkable since the structure of the aggregates and their spatial distributions changed considerably from one system to another. We have also compared the prediction of Equation (52) with some unpublished experimental data on dodecane-in-water nano-emulsions (φ = 0.20, T = 25 C) obtaining good results. For this purpose the average radius of the emulsion was measured as a function of time. The number of particles was estimated dividing the total volume of oil by the average volume of a drop at each time. Based on the results of the previous section (4.1), we approximated the number of aggregates by the number of particles during the FC period, and used Equation (52) to fit the results. The fitting is shown in Figure  5 along with the experimental data (green circles). Notice that the experimental curve shows a sharp decrease in the FC rates after a pronounced initial drop in the number of particles. Equation (52) also appears to fit the data beyond the FC interval.   (52) around t = 0, showed that these two functions are very similar. In fact, it is possible to exchange the pair of coefficients corresponding to each process (A, k 1 ) ↔ (B, k 2 ) and obtain a fitting of similar quality. This demonstrated that the processes of flocculation and coalescence are mixed in the kinetic rates of Equation (52). Moreover, it was also observed that Equation (52) was able to fit a terminal aggregation stage in those systems in which the number of particles stabilises after a certain period of time.
It was also found that the interfacial area of the emulsion after ~200 seconds was proportional to the number of surfactant molecules. The slope of this curve provides a measurement of the area per surfactant molecule that is required (~ 236 Å 2 ) in order to avoid an initial pronounced drop in the number of particles. Using the effective charge of a surfactant molecule (0.21e), the corresponding value of the surface charge can be obtained (σ = 13.95 mCoul/m 2 ). The electrostatic surface potential can also be calculated (Ψ 0 ~44 mV) using the relation between σ and Ψ (Equation (26)). Unexpectedly, these values do not correspond to a positive magnitude of the potential energy for the referred systems. Instead, they correspond to the appearance of a shoulder at negative values of the potential energy (see Figure 6-7). This shoulder originates the secondary minimum (for which 0 along with a small potential barrier. Yet, this barrier occurs at negative values of the interaction potential where an attractive force between the particles exists.    From the variation of n a vs. t (or N a vs. t) in the systems studied (see Figure 4) it is clear that a considerable stabilisation is reached shortly after the FC interval. At this point, several factors may favour the decrease of the flocculation rate. The number of remaining drops after the FC period is smaller and the mean free path between the drops is larger than it was at the beginning of the simulation. Moreover, the diffusion constant of the remaining drops is smaller since it is inversely proportional to the radius of the drops. However, these factors are likely to produce a progressive decrease of the FC rate as is observed in Figure 8. Instead a sharp change in the behaviour of n a was observed in most systems studied in Ref. [73].
We believe that the drastic change in the slope of n a vs. t is basically caused by the increase of the repulsive force between the remaining particles of the emulsion. In this regard it should be noticed that when the particles execute Brownian motion, their kinetic energy is lost in a short period of time. As a result, the total energy of the drops should approach closely the potential energy curve. Hence, the absence of a driving force at the secondary minimum ( 0 considerably the movement of small drops. In this situation a small repulsive force could be enough to prevent primary minimum flocculation (coalescence).  A more conventional effect of the repulsive barrier was observed in the simulation of hexadecanein-water nanoparticles (R i = 200 nm) stabilized with nonyl phenol ethoxylated (NPE) surfactants. In these calculations the repulsive potential was assumed to be steric, and the interfacial area of the surfactant molecules was kept constant during the whole simulation. These calculations showed that there is an analytical relationship between the height of the repulsive barrier of the interaction potential and the stability ratio: This equation is very similar to Equation (9) except for the method of evaluation of V Δ and W.
Equation (9) was deduced from the systematic evaluation of Equation (8) for the case of solid particles suspended in water. Instead Equation (53) was deduced from emulsion stability simulations of liquid drops exposed to flocculation and coalescence. In the present case V Δ was measured from the base of the secondary minimum up to the top of the potential barrier. In the former case the height of the barrier was measured from V = 0. Moreover, Equation (8) Equation (54) results from the inverse proportionality between the half-lifetime and the flocculation rate suggested by Equation (11). However, it does not depend on the particular form of the curve of n a vs. t.

High surfactant concentrations
Equation (8) does not impose any restriction to the value of W: the higher the potential barrier the higher the stability ratio. However, according to ESS of non-deformable drops [44], the number of particles of the emulsion is preserved when the repulsive barrier between the drops is higher than T k B 7 . 12 . In this case the dynamics of the system depend on the depth of the secondary minimum of the interaction potential [21]. According to Verwey and Overbeek [20] a 25-k B T barrier is necessary for a 7-day stabilization of a concentrated emulsion (n 0 ~ 10 14 m -3 ). This corresponds to W = 10 9 . Much higher values of the stability ratio (10 10 -10 17 ) are found in the literature published between 1917 and 1940 (see Refs. [21,76] and citations therein). These high stability ratios mostly correspond to sols of very small particles (between 20 nm and 200 nm). Much lower values of W (1 ≤ W ≤ 100) are found for latex suspensions (see Table 1 in Ref. [21]).
The wide range of stability ratios reported during the past century is partially caused by the inadequate definition of W. The first expression of W deduced by Fuchs [24]  ). However, even when the correct expression of W is used (Equation (8)), its value rarely surpasses 1000 for particles between 0.5 μm and a few microns. These low values of W do not agree with the fact that most of these particles exhibit a high repulsive barrier against primary minimum flocculation (Equations (9)). We studied this problem in Ref. [21] using several particle sizes and interaction potentials. Varying the surface charge of the particles and the ionic strength of the solution, it was evident that particles between 10 nm and 100 nm usually show very shallow secondary minima or no secondary minima at all. Micron size particles instead show deep secondary minima.
The simulations indicated that the systems with barrier heights higher than 20 k B T preserve the number of particles (coalescence does not occur during the extent of the simulation). In these cases the depth of the secondary minimum determines the evolution of the system. On the one hand, shallow minima do not lead to aggregation [21,74]. Hence, the total number of aggregates (n a ) oscillates frequently around the number of particles indicating the sporadic formation of doublets and their quick dissolution. On the other hand, deep secondary minimum promotes fast flocculation. In this case irreversible secondary-minimum flocculation occurs. The random fluctuation of the number of aggregates previously observed with shallow minima does not occur. The particles do not go in and out the secondary minimum as it happened in the former case. The curves of n a vs. t show a progressive decrease as a function of time, as it is usually exhibit in the case of fast aggregation. It is important to realise that deep secondary minimum often occur at high ionic strength, where the repulsive barrier and the secondary minimum are closest to the particle surface.
Recently Kuznar and Elimelech studied the deposition of micrometer-sized particles on a single layer of packed glass beads using a flow-cell [77]. The colloidal particles used had a mean diameter of 4.1 μm. The diameter of the glass spheres was 1 mm. The shapes of the DLVO potential between the particles and the glass beads (Figs. 6-7 in Ref. [77]) were very similar to those reported by our group in Figure 1 of Ref. [21]. In the experiments of Kuznar and Elimelech the depth of the secondary minimum varied from -0.62 k B T at 1 mM KCl to -21.7 k B T at 30 mM. The repulsive barrier of the DLVO potential changed from 5029 k B T (at 10 mM) to 1089 k B T (at 30 mM). Notice that huge values of W result from introducing these potential barriers into Equations (9) or (53). Hence, these high repulsive barriers should prevent the deposition of particles over the glass beads. However, this behaviour was only observed at 0 mM KCl where the highest repulsive barrier occurred. At 100 mM the repulsive barrier was completely screened by the counter ions, and large amounts of latex particles were deposited. These two findings agree with the theory of DLVO. However, between 10 mM and 30 mM increasing amounts of latex particles were deposited. This behaviour completely contradicts the predictions of DLVO. Moreover, when the ionic strength was decreased after the deposition of the particles, most of them were washed out, but some of them remained attached to the glass beads. In a previous article [78] Tufenkji and Elimelech reported deviations from the classical colloid filtration theory in the presence of repulsive DLVO interactions. According to these authors, the fast deposition of microbial particles on porous columns was caused by the combined effect of favourable and unfavourable colloidal interactions. The authors related the deposition rates to the influence of deep secondary minimum in the aggregation process, proposing a dual-deposition model that includes slow and fast aggregation rates.
The experimental evidence described above strongly supports the predictions of Ref. [21]. Hence, it is possible to obtain fast secondary-minimum flocculation in systems of micron-size particles even in the presence of a high repulsive barrier. Figure 9 shows unpublished experimental data regarding the aggregation of 96-nm sulfonated latex particles. This latex was synthesised by Ms. K. Rahn and Dr. A. Lozsán at the University of Almería, Spain, under the tutorage of Dr. M.S. Romero-Cano. The experimental values of W shown in Figure 9 (blue symbols), were evaluated from the initial slope of the absorbance (Ab) of latex suspensions as a function of time [79]: In Equation (55) subscript "csalt" stands for the salt concentration of the solution (ionic strength). Theoretical values of W were estimated following a methodology similar to the one reported in Ref. [80]. In the present simulations we take advantage of the fact that primary minimum flocculation implies coalescence for non-deformable droplets. Hence, we approximated the flocculation time between two latex particles in the experimental set up, by the coalescence time between two nondeformable drops with the same interaction potential. The drops were initially located 20 nm apart from each other. The average coalescence time -at each ionic strength-was estimated from 100 simulations, which only differ in the initial value of the random number generator. These calculations are the simulation analogue of the numerical evaluation of Equation (8). While many-particle calculations give information about the evolution of the emulsion and the structural characteristics of the aggregates formed, two-particle simulations are restricted to the evaluation of the time required for primary-minimum flocculation. However, these calculation are fast and provide all the information necessary for the evaluation of W.
As shown in Figure 9, the value of the stability ratio tends to 1 at high ionic strengths and steeply increases below 570 mM NaCl. The values of W produced by the simulations vary non-monotonously around the experimental measurements. This could be caused in part by the low number of calculations used to compute the average values of the coalescence time. Notice also that we could not evaluate the stability ratio corresponding to 400 mM NaCl, because in this case the repulsive barrier was higher than the coalescence threshold of T k B 7 . 12 ( ). The calculation corresponding to 400 mM NaCl has been running since several months ago, and not a single coalescence event has been observed. In our view, the agreement between theory and experiment shown in Figure 9 is reasonable. This suggests that the coalescence times determined by ESS for systems of non-deformable particles should be meaningful.

Time-dependent adsorption
On the one hand a 0.1 g/L solution of sodium dodecyl sulfate (SDS) causes a drop of 15 mN/m in the interfacial tension of a dodecane/water system in a period of 100 s [81]. On the other hand a diluter solution of SDS (C T = 0.014 g/L) achieves the same fall in only 6 s in the presence of 0.1 g/L NaCl. Which chemical condition is more effective for stabilising a dodecane-in-water emulsion considering that the electrostatic interaction is screened in the latter case?
In order to throw some light into this problem, we computed the evolution of a set of oil-in-water emulsions with volume fractions between 0.10 ≤ ≤ φ 0.40, apparent diffusion constants in the range 10 -12 m 2 /s ≤ ≤ app D 10 -9 m 2 /s, and two surfactant concentrations C T = 1 x 10 -4 M, and C T = 5 x 10 -4 M [55]. These parameters cover an ample range of experimental situations. Non ionic surfactants usually show diffusion limited adsorption [58,82]. A typical value for their diffusion constant is 10 -10 m 2 /s, unless their molecular weight is very large. Ionic surfactants show kinetically limited adsorption [59]. This means that the molecules approaching the interface experience the electrostatic repulsion of those previously adsorbed. Hence, the apparent diffusion constant of these molecules could be substantially lower than 10 -10 m 2 /s. Furthermore, according to Rosen et al. [60,61] a concentration of 5 x 10 -4 M is the lowest surfactant concentration required to achieve a 1-second surface tension reduction that does not change much if the surfactant concentration is increased. In these simulations the surfactant was assumed to be ionic. Therefore, the interaction potential between the drops was supposed to be DLVO. (57) Table 1 shows the values of t msa as a function of app D and C T . The value of app D depends on the molecular structure of the surfactant and other environmental parameters as the salt concentration, the temperature of the system, and the type of oil.  If the surfactant concentration is high C T = 5 x 10 -4 M, and its diffusion constant fast app D = 10 -9 m 2 /s, the initial drop size distribution of the emulsion is preserved [55]. In the opposite case, app D = 10happen between these extreme situations. They show an initial decrease in the number of drops due to the early destabilisation of the former DSD. Destabilisation occurs faster when the mean free path between the drops is small (high volume fraction, inhomogeneous spatial distribution of drops, etc.). As time passes, surfactant adsorption progressively increases until it reaches a point in which the repulsive potential generated by the surfactant is enough to preserve the current DSD. At this point the kinetic rates of coalescence and flocculation markedly change. The variation of number of aggregates as a function of time follows the prediction of Equation (52) for a completely different reason as the one discussed before. The repulsive force between the drops increases with time, not as a result of the surfactant redistribution among the drops, but due to the augment of the surfactant adsorption as a function of time. Notice also that in the case of non-deformable drops the diffusion constant decreases with R i and the repulsive forces augment. This might also cause a progressive stabilization of the system as a function of time.
In our calculations, C T = 5 x 10 -4 M is the lowest surfactant concentration able to generate a substantial repulsive potential for app D ≥ 10 -10 m 2 /s [55]. However, it must be kept in mind that in the referred simulations the total surfactant population is distributed among the surfaces of the drops. Hence, the mass balance does not include the surfactant solubility in the water phase. As a result it is expected that the surfactant concentration required for the attainment of the fastest stabilisation of the emulsion should be higher than 5 x 10 -4 M.
As Equation (34) shows, a slow surfactant adsorption can be compensated with a high surfactant concentration. However, a low surfactant concentration can only be partially compensated with the increase of app D . That is the case of the systems with C T = 10 -4 M and φ ≥ 0.20 [55]. A very low surfactant concentration does not generate a high repulsive force even if instantaneous adsorption occurs.
We also observed a marked effect of the hydrodynamic interaction in the most concentrated system (φ = 0.40). This system was expected to experience the highest destabilisation due to the short mean free path between the drops. However, this was partially the case, since the hydrodynamic interaction slowed down the flocculation of the drops considerably, producing a substantial lag time between the beginning of the simulation and the first coalesce event.
It should be remarked that the effect of the ionic strength was not considered in Ref. [55]. On the one hand, a high surface coverage generates a substantial repulsive force between the drops at low ionic strength (Figure 1). On the other hand, a high ionic strength (~ 1 M NaCl) can screen completely the electrostatic interaction, eliminating any effect of dynamic adsorption on emulsion stability, for the case of an ionic surfactant. However, according to the results of Figure 9, an ionic strength as high as 0.55 M might not be not enough to eliminate the repulsive barrier between particles of nanometer size ( T σ = 39 mCoul/m 2 , A H = 4.27 x 10 -21 J, W = 1.22).
Furthermore, it should be noticed that the value of the Hamaker constant used in the simulations of Ref. [55] was meant to represent the attraction between bitumen drops stabilized with a cationic surfactant [83]. A more typical value of A H for oils (~10 -21 J) should decrease the flocculation rate, diminishing the effect of the time dependent adsorption even more.
Taking into account the results of the simulations [55] and all the limitations outlined, the effect of dynamic adsorption is only expected to be relevant at low surfactant concentrations, C T ≤ 10 -4 M, mostly between 0.20 ≤ ≤ φ 0.30. More precision requires further calculations employing the specific characteristics of the experimental system. These include not only the exact values of φ , C T , and app D , but also the initial DSD of the emulsion, the ionic strength of the aqueous solution, and the solubility of the surfactant (or its adsorption isotherm).

Ostwald Ripening
The Liftshitz-Slezov-Wagner (LSW) theory of Oswald ripening, predicts a linear variation of the cube of the average radius of a dispersion as a function of time ( 3 c R vs. t) [84][85]. It also envisages a left-skewed drop-size distribution with a cut-off radius of 1.5 c R . These predictions are obtained assuming that: a) The particles are fixed in space. b) The system is infinitely dilute (implying the absence of interactions). c) The molecules of the internal phase are transported from one particle to another by molecular diffusion. d) The concentration of oil is the same through the whole system except in a direct neighbourhood of the particles, where it is given by the Kelvin equation.
LSW does not provide analytical expressions for the variation of c R during the transient period of ripening, nor does it consider the effects of flocculation and coalescence that simultaneously occur. At long times (in the asymptotic limit) LSW predicts a constant value for the Oswald ripening rate, V OR , equal to: Sakai et al. [86] studied the evolution of the DSD for several alkane/water systems in the absence of a surfactant. For a dilute dodecane in water emulsion ( ≈ φ 2.29 x 10 -4 ) a OR V rate equal to 4.0 x 10 -26 m 3 /s was found. This rate is three times larger than the theoretical estimation (1.3 x 10 -26 m 3 /s). Moreover, the DSDs of all measured systems were skewed to the right, in contradiction with the predictions of LSW. Using freeze-fracture electron microscopy, Sakai et al. [87] also reported direct observations of flocculation and coalescence of metastable benzene droplets under surfactant-free conditions. They observed small drops with diameters at 30-100 nm immediately after sonication, as well as aggregates of medium size (200-500 nm) composed of small droplets.
In Ref. [71] we revisited the evolution of a dodecane-in-water nanoemulsion in the absence of surfactant molecules. The research included simulations and experimental results. The main purpose of the investigation was to determine the effect of flocculation and coalescence in the temporal variation of 3 c R .
The first set of simulations followed the original algorithm of De Smet [67]. This code does not consider the movement of the drops. The initial DSD was Gaussian with According to these calculations, c R does not increase prior to 800 s as a consequence of molecular diffusion.
During this transient period the average radius of the emulsion slightly decreases due to the exchange of oil molecules between the drops. Moreover, c R does not increase until the first small particle dissolves. Hence, any increase of ( ) We carried out similar calculations with the ESS program including Ostwald Ripening (ESS+OR simulations). For these simulations 1000 non-deformable drops were used. The potential was completely attractive, characterized with a Hamaker constant, A H = 5.02 x 10 -21 J. When the movement of the particles is incorporated, the simulations showed that dt R d 3 changes linearly with time as a result of flocculation and coalescence: dt R d 3 = V FC = 3.11 x 10 -22 m 3 /s (r 2 = 0.9979). This rate is four orders of magnitude higher than the theoretical value of OR V and cannot be attributed to Oswald ripening. Moreover, right-skewed distributions of particles were obtained [71]. This finding motivated the measurement of c R at short times. Dilute O/W dispersions of dodecanein-water (φ = 2.3 x 10 -4 ) were prepared using an aqueous solution of NaCl (0.5 M) to screen the electrostatic charge produced by the preferential adsorption of hydroxyl ions. Immediately after sonication, the light scattered by the emulsion was measured at 90º using a BI-200SM goniometer (Brookhaven Instruments). The mean diffusion coefficient was derived from the intensity autocorrelation function using a cumulant analysis [88]. Figure 10 shows the result of one of these measurements. A value of 2.6 x 10 -23 m 3 /s was found for dt R d 3 . This value differs in one order of magnitude from the one of the simulations, but is three orders of magnitude away from the theoretical value of ripening. It also suggests the existence of a small repulsive barrier that hinders flocculation and/or coalescence even at high ionic strengths. In Section 3 we mentioned that in the absence of a strong repulsive force, the expression for the total number of aggregates of a dispersion (Equation (49)) is the same for suspensions and emulsions except for some minor restrictions regarding the counting of the aggregates. In the same reference [75] Danov et al. demonstrate that the formula for the concentration of i-particle aggregates is also equal to the one deduced by Smoluchowski if the coalescence time is infinite, or what is the same, if only flocculation occurs: The average size of the particles in the system can be calculated substituting Equation (60)

Time (s)
Smoluchowski supposed that the dispersion was initially composed by particles of equal radius. Hence, some prescription was necessary to connect the size of the aggregates with the initial particle radius ( 0 R ). For example: For coalescing emulsions l should be equal to 3. However, only certain values of l allow obtaining an analytical expression for R . A value of l = 1, consistent with the formation of linear aggregates, produces a very simple formula: If the aggregates formed are not linear and/or coalescence occurs, Equation (63) is not valid and constitutes a very rough approximation to the mean particle radius. However, if we assume the validity of Equation (63) where F V stands for the velocity of flocculation. The short-time limit of Equation (64) provides a very useful relation between F V and f k : Notice that a real exponent p in Equation (63): , still produces the same short- ( ) The data shown in Table 2 corresponds to the ESS + OR simulations referred above. As expected V F is proportional to k f . However, the quotient between V F and 3 0 0 R n k f does not appear to be related to any of the exponential factors suggested by Equations (65)- (66). Hence, it appears that the effect of the coalescence on the value of t d R d 3 is significant. Yet, there appears to be a proportionality between the value of the flocculation rate and dt R d 3 .

Modifications of ESS for the Calculation of Deformable Drops
The problem of drop deformation is formidably complex. Deformation results from the interplay between hydrodynamic and interaction forces. Hence it depends on all the parameters that affect these forces. As soon as deformation occurs, the geometrical part of the interaction potentials changes. Analytical forms for the potentials of truncated spheres are available in the literature, but they are function of the radius of the film between flocculated drops and its width. Moreover, the process of deformation itself involves several stages including the formation of a dimple, its evolution to a plane parallel film and the destabilisation of the film either by surface oscillation or by the formation of holes [91]. This complexity is additional to the one resulting from the simulation of all destabilisation processes that occur in emulsions. Hence, some drastic approximations are necessary. The methodology outlined in this section mostly corresponds to the thesis of Dr. Toro-Mendoza [92,93].

Regions of deformation
The present version of the ESS code has routines for non-deformable and deformable drops, but cannot simulate both types of drops during the same calculation. Thus, if the mode of deformable droplets is selected, it is assumed that the deformation of the drops occurs independently of the energy required for this process. Both types of simulation use the same equation of motion (Equation (12)) and follow the algorithm illustrated in Figure 2, but they differ in the way of calculating the diffusion constants, and the forces. Moreover, they differ in the criteria employed for the coalescence of drops.
In the mode of deformable droplets, the drops move as spheres until they reach the initial distance of deformation h 0 . When this occurs, a model of truncated spheres is used [6,38,[94][95]. The transient time between the formation of the dimple and the surface oscillations that lead to the formation of a thin liquid film between flocculated drops is not taken into account. As soon as ij r < 0 h R R j i + + the code calculates the dimensions of truncated spheres that are compatible with h 0 and ij r . Using the resulting film width and film radius, the program computes the potential of interaction between truncated spheres. In the absence of other coalescence mechanisms different from film drainage, the drops coalesce if the separation between their surfaces reaches a critical distance h = h crit . In order to calculate the forces and diffusion constants three regions of approach are defined: Region I: The distance of separation between the centres of mass of the drops r ij, is greater than . In this case the drops maintain their spherical shape. Consequently: Region II: It covers the range of distances between the beginning of the deformation film r ≠ 0, and the attainment of the maximum film radius: . Within these limits, the closest distance of separation between the surfaces of the drops is constant h = h 0 [6], and: Region III: The film already attained its maximum radius, , and it progressively drains until reaching a critical distance of approach. Within these separations: Notice that the drops behave as spheres in Region I (r ij > 0 h R R j i + + ) even in the case in which the mode of deformable drops is selected. This means that the potential of interaction and diffusion constant in Region I correspond to spherical particles. When the distance of separation between the surfaces of the drops is less than h 0 (Regions II and III), the expressions of the potentials corresponding to truncated spheres are employed [38]. As explained in the next section, two additional potentials appear during the evolution of the film (Region II). They take into account: (a) the increase of interfacial area between a sphere and a truncated spheroid (extensional potential), and (b) the change of curvature of the interface (bending elasticity potential). These potentials contribute to the total potential of interaction and its force within Region II, but they assume a constant value in Region III and do not contribute to the value of the force in this zone.

Interaction forces between deformable droplets
As soon as deformation occurs, the analytical forms of the interaction potentials change. Analytical forms of the potentials for truncated spheres are available in the literature. However, they are expressed as a function of the radius of the film between flocculated drops and its width [6,38,95]. For example, the geometrical part of the van der Waals potential for truncated spheroids is equal to: where: , film r is the radius of the thin liquid film between flocculated drops, and h is the closest separation between the surfaces of the drops. Different types of interaction potentials for truncated spheres are found in the literature [38,95]. The researchers of the University of Sofia deduced most of them. They include van der Waals, electrostatic, oscillatory, depletion, steric, etc. These mathematical equations are expressed in terms of h and film r . Use of Equations (67)-(72), allows to recast those potentials in terms of r ij . Hence, it is possible to differentiate the potentials with respect to r ij using a package of symbolic algebra. In this way, the force is obtained for each deformation zone. Additionally, two new contributions to the free energy appear. The first one is called dilatational or extensional energy and is caused by the increase of the interfacial area of the particle as it looses its spherical shape. The analytical form for the extensional potential corresponding to the change between a sphere and a truncated spheroid is [38,[94][95]: where i film r , is the radius of film formed by the truncated sphere i.
The variation of the interfacial curvature requires an additional amount of energy. This contribution can be positive or negative depending on the value of the spontaneous curvature of the interface: where R c,0 is the radius of curvature of the surface adopting its lowest free energy configuration. The energy of interfacial bending [38,[94][95] is equal to: where H stands for the actual curvature of the interface: and B 0 is the interfacial bending moment of a flat interface: Constant B 0 is related to the bending moment in the theory of Helfrich [96]: k b . Theoretical estimations suggest a value of 5 x 10 -11 N for B 0 . However, this value produces elastic potentials higher than 600 k B T for micron-size drops. In current calculations we adjust the value of B 0 in order to reproduce the expected bending energy of a 100 nm drop (31 k B T) [6]. This suggests a value of B 0 of 1.6 x 10 -12 N.
The potentials of surface deformation occur at close separation distances, within or nearby the same region of influence of the other potentials of interaction like electrostatic, steric, etc. In the absence of other repulsive potentials, deformation (dilatational + bending) generates two minima in the interaction potential separated by a repulsive barrier, similar to the ones exhibit by the DLVO potential (see Figure 1 and Figure 11). The presence of other repulsive contributions changes the smooth form of the total potential, increasing its value at short distances. This may also produce additional peaks and minima. Whatever the case, the passage of the particles over the closest repulsive barrier implies coalescence.
It is not yet possible to deduce a general analytical expression to relate h, film r , and ij r for all regions of deformation. Such an expression will round the sharp peaks of the total potentials shown in Figure 11, eliminating their present discontinuities. The acute changes of slope exhibited, result from the segmentation of the potential into three Regions (Equations (67)- (72)). However, all contributing potentials are differentiable within each domain and should not generate ill-define forces at the discontinuities. It should be noticed, that Equations (74)-(78) rely on the interfacial tension, which in turn depends on the surface excess (Γ). Moreover, both free energies require the evaluation of the radius of the O/W/O film between flocculated drops ( i film r , ). This radius depends on the properties of the interfacial layer. Hence, in order to evaluate the forces in Equation (12) the number of surfactant molecules adsorbed to the interfaces of the drops has to be previously determined.
Accurate estimation of h 0 and h crit is very difficult since they result from a balance between hydrodynamic and interaction forces. On the one hand, a formal prescription for their evaluation is given in Ref. [97]. That type of calculations is too expensive in computational time, and cannot be incorporated into ESS. On the other hand, emulsions are polydisperse systems and the value of film r depends on the size of the flocculating drops and their interfacial tension. Figures 3 and 4 of Ref. [97] show the value of h 0 resulting from free energy calculations. We obtained an empirical form for h 0 in terms of the radius of a drop and its interfacial tension, h 0 (R i ,γ), fitting the referred curves with an empirical polynomial:

Calculation of the diffusion constant for deformable droplets
The form of the diffusion constant employed at every distance of separation corresponds to the shape of the approaching drops. Within Region I the same formalism used for non-deformable drops is employed [30]. For zones II and III several expressions are available in the program. They correspond to the rates of thinning T V of thin liquid films reported by Gurkov and Basheva [102]. These velocities depend on the force between the particles F, h and film r . Thus, they can be transformed into effective diffusion constants using the Einstein relation: , where f is the friction coefficient (2)).
Most equations compiled by Gurkov and Basheva [102] require knowledge of the properties of the surfactant and/or the interfacial layer (Gibbs elasticity, viscosity, etc). Dr. Toro-Mendoza suggested Equation (81) based on a slight modification of the methodology previously used by our group for spherical drops [30].
Thus, the analytical form of ) 2 ( corr f proposed by Honig et al. [31] for spherical drops (expression in parenthesis in Equation (81)), used for int r r ij ≤ is modified as soon as the drops change their spherical shape to form truncated spheres with a plane parallel film. Hence, the additional friction generated by the creation of two planar disks between flocculated drops, decreases the diffusion constant beyond the estimation of Honig et al.

Coalescence criteria for deformable drops
The process of coalescence of deformable droplets is involved. It is known that thin liquid films not necessarily drain until reaching h crit [7][8][9][10][103][104][105]: Depending on the properties of its interfacial layers, films can collapse through the formation of holes [9][10] or due to the enhancement of surface oscillations [7][8]105]. Stable Newton black films can also be formed, prolonging the stability of the flocs for months. In the case of small drops of nanometer size it seems unlikely that deformation takes place due to their high internal pressure.
By default, the mechanism of coalescence of deformable drops is the one of film drainage. As the drops approach, a plane parallel film forms and thins until h crit is reached. At this point, a new spherical drop is generated following Equation (17).
It is known that thin liquid films of macroscopic radii experience surface oscillations. In some circumstances the interfacial waves at each O/W interface may grow until they touch. In this case, a channel between the oil phases of the drops could form causing their coalescence. This mechanism is additional to the one of film drainage. Other possible process like the condensation of holes in common black films are also known to promote film rupture [9][10]. As in the case of surfactant adsorption, we use an algorithm which mimics the actual process. We do not attempt to simulate these processes in detail.
If the additional mechanism of surface oscillations is selected, a random number is assigned to the surface of each drop every time a pair of particles enters Regions II or III above. The random number varies between -1.0 and 1.0 units.
The amplitude (height) of each capillary wave ( λ ) at each interface is estimated as the product of the referred random number times the value of h crit .
Equation (84) was deduced assuming van der Waals interactions only. A more general expression requires knowledge of the second order differential of the free energy in terms h under certain restrictions. This differential changes abruptly for the case of deformable droplets. The program has the additional option to introduce the value of Vrij τ as input.
The value of the tension in Equation (84)  Coalescence occurs whenever the total height of the surface oscillations is greater than h crit . The total height of the surface oscillations is approximated by: Equation (86) takes into account the fact that the surface oscillations increase exponentially with time [7][8]105].
It is clear from above, that all coalescence mechanisms are defined in terms of doublets. If aggregation of multiple particles occurs, Equation (86) is applied to each film formed between the flocculated particles.
The force routine is divided into two parts to account for non-deformable and deformable droplets. In the case of deformable droplets, the program estimates the value of h 0 for each drop according to Equation (79). Using this value and the distance of separation between the centres of mass (r ij ), the region of approach can be determined (Regions I, II and III). The analytical expression for the force corresponding to the actual region of approach between each pair of drops is used.

Preliminary Results
The incorporation of deformable droplets is very recent [92][93]. We began the study of deformable droplets with the calculation of the coalescence time between two particles under the influence of van der Waals, extensional and bending potentials. All calculations began from a separation distance equal to h = h 0 (5 nm ≤ ≤ 0 h 15 nm). The particle radius was changed between 100 nm and 100 μm, for γ = 1 mN/m. Within these limits, the diffusion constant used (Equation (81)) markedly decreases as a function of the particles radius [93]. The coalescence times calculated for each particle radius correspond to the average value of 1000 simulations which only differ in the initial value of the random number generator. Only the mechanism of film drainage was considered for the coalescence of the drops. For the 96-nm particles of Figure 9 The coalescence time between two non-deformable droplets increases steadily with the size of the particles. The coalescence time of deformable droplets increase between R i = 100 nm and R i = 5 μm (range 1), and between R i = 10 μm and R i = 100 μm (range 3). However, it decreases between 5 μm and 10 μm (range 2), generating a small peak in the curve of t c vs. R i . This particular behaviour is due to the opposite trends of interaction forces and hydrodynamic effects. Within range 1 the diffusion constant decreases abruptly causing an increase of the coalescence time as a function of the particle size. From 5 μm to 10 μm the diffusion coefficient also decreases but in a much lower rate. Hence the attractive force prevails decreasing t c as a function of R i . Within range 3, the interaction is still attractive. However, the values of h 0 and For φ = 0.10 the value of FC k for non-deformable droplets (3 x 10 -17 m 3 /s) is two orders of magnitude higher than the one of deformable droplets (5 x 10 -19 m 3 /s). For φ = 0.30 the difference reduces to only one order of magnitude (3 x 10 -16 m 3 /s and 4 x 10 -17 m 3 /s, respectively).

Conclusions
This review describes in detail the algorithm of Emulsion Stability Simulations developed by our group. It was not until very recently that we completed the routines for deformable droplets and Ostwald ripening. Consequently, it is now that we are ready to study the role of the surfactant molecule in systems exhibiting the concurrent occurrence of most destabilization processes. This is why we are confident that the most outstanding results from ESS are yet to be produced. However, we had shown that the previous studies related to the behaviour of non-deformable drops are very insightful. In particular they had remarked the role of the secondary minimum in the flocculation of drops, suggested a new threshold for the coalescence of drops, and demonstrated that the cube of the average radius of an emulsion can change linearly with time as a consequence of flocculation and coalescence.