Timescales of emulsion formation caused by anisotropic particles

Particle stabilized emulsions have received an enormous interest in the recent past, but our understanding of the dynamics of emulsion formation is still limited. For simple spherical particles, the time dependent growth of fluid domains is dominated by the formation of droplets, particle adsorption and coalescence of droplets (Ostwald ripening), which eventually can be almost fully blocked due to the presence of the particles. Ellipsoidal particles are known to be more efficient stabilizers of fluid interfaces than spherical particles and their anisotropic shape and the related additional rotational degrees of freedom have an impact on the dynamics of emulsion formation. In this paper, we investigate this point by means of simple model systems consisting of a single ellipsoidal particle or a particle ensemble at a flat interface as well as a particle ensemble at a spherical interface. By applying combined multicomponent lattice Boltzmann and molecular dynamics simulations we demonstrate that the anisotropic shape of ellipsoidal particles causes two additional timescales to be of relevance in the dynamics of emulsion formation: a relatively short timescale can be attributed to the adsorption of single particles and the involved rotation of particles towards the interface. As soon as the interface is jammed, however, capillary interactions between the particles cause a local reordering on very long timescales leading to a continuous change in the interface configuration and increase of interfacial area. This effect can be utilized to counteract the thermodynamic instability of particle stabilized emulsions and thus offers the possibility to produce emulsions with exceptional stability.


I. INTRODUCTION
Particle stabilized emulsions play an important role in pharmaceutical, food, oil and cosmetic industries [1].The particles are adsorbed at the interface between two immiscible fluids and as such stabilize the emulsion.The stability of the emulsions depends on several parameters like particle coverage at the interfaces and the wettability of the particles.It was found that the particle coverage at the interface is the most important parameter for stabilizing emulsions [2].The colloidal particles act in a similar way as surfactants.In both cases the free energy of the interface is reduced.However, the fluid-fluid interfacial tension is not being modified by particles [3].
Several types of particle stabilized emulsions are known including the bicontinuous interfacially jammed emulsion gel (bijel) and the more widely known Pickering emulsion.The Pickering emulsion was discovered in the beginning of the 20 th century independently by Pickering and Ramsden [4,5].It consists of discrete particle covered droplets of a fluid immersed in a second fluid.The bijel was predicted in 2005 by simulations and experimentally realized for the first time in 2007 [6,7].It consists of two continuous phases.The choice of control parameters such as particle concentration, particle wettability and ratio between the two fluids determines if a bijel or a Pickering emulsion is obtained [8,9].There are * j.harting@tue.nlmany kinds of particles/colloid types which can stabilize an emulsion.I.e., next to spheres [10,11], the colloidal particles can also be of more complex nature and include anisotropic shapes [12], magnetic interactions [13,14], or anisotropic Janus style properties [15].
The influence of the particle shape on the stabilization of Pickering emulsions was studied experimentally with prolate and oblate ellipsoids, e.g. in Ref. [16].As the degree of the particle anisotropy increases, the effective coverage area increases.In this way they are more efficient stabilizers for emulsions than spherical particles.Furthermore, the rheological properties of the emulsion vary with changing aspect ratio because the coverage of the fluid interfaces and the capillary interactions differ.
In Refs.[8,[17][18][19][20] the adsorption of a single particle at a flat interface is studied in absence of external fields.The stable configuration for elongated ellipsoids is the orientation parallel to the interface [8].This state minimizes the free energy of the particle at the interface by reducing the interfacial area [17,19,20].If the particle shape is more complex like e.g. the super-ellipsoidal hematite particle [21], several equilibrium orientations are possible.
Furthermore, if particles are adsorbed at an interface they generally deform the interface.This deformation can be caused for example by particle anisotropy [22], external forces such as gravity or electromagnetic forces acting on the particles [23,24], or non-constant interface curvature [25].This deformation leads to capillary interactions between the particles.In case of ellipsoids at a flat interface it is a quadrupolar potential [26], which leads to spatial ordering [27].
In general, particle stabilized emulsions are thermodynamically unstable and just kinetically stable.The energetic penalty for creating the interface is much higher than the entropic increase.While thermodynamic stability for emulsions has been reported in some special cases, one can generally assume that this requires the interplay of several effects such as particle interactions due to charges, amphiphilic interactions (Janus particles) or additional degrees of freedom [28][29][30].
Due to the short timescales and limited optical accessibility, the dynamics of the formation of emulsions has only found limited attendance so far [31].The focus of the current article is to study the influence of the geometrical anisotropy and rotational degrees of freedom of ellipsoidal particles on the time development of fluid domain sizes in particle-stabilized emulsions.To obtain a deeper understanding of the individual contributions to the stabilization and formation process due to the particles we investigate model systems involving either a single particle or particle ensembles at a simple interface.We will demonstrate that the rotational degrees of freedom of ellipsoids can have an impact on the domain growth and might be a suitable way to generate particle stabilized emulsions with exceptional long-term stability.
This article is organized as follows: the simulation method is introduced in section II.Dynamic emulsion properties are studied in section III.Sections IV and V discuss a single particle and a particle ensemble at a flat interface, respectively.Section VI describes the behavior of a particle ensemble at a spherical interface.We finalize the paper with a conclusion.

A. The lattice Boltzmann method
For the simulation of the fluids the lattice Boltzmann method is used [32].The discrete form of the Boltzmann equation can be written as [3] where f c i (x, t) is the single-particle distribution function for fluid component c with discrete lattice velocity c i at time t located at lattice position x.The D3Q19 lattice with the lattice constant ∆x for three dimensions and with nineteen velocity directions is used.∆t is the timestep and is the Bhatnagar-Gross-Krook (BGK) collision operator [33].
The density is defined as ρ c (x, t) = ρ 0 i f c i (x, t) where ρ 0 is the proportionality factor of the density.τ c is the relaxation time for the component c and is the third order equilibrium distribution function.
is the speed of sound, u c = i f c i (x, t)c i /ρ c (x, t) is the velocity and ζ i is a coefficient depending on the direction: ζ 0 = 1/3 for the zero velocity, ζ 1,...,6 = 1/18 for the six nearest neighbors and ζ 7,...,18 = 1/36 for the next nearest neighbors in diagonal direction.The kinematic viscosity can be calculated as In the following we choose ∆x = ∆t = ρ 0 = 1 for simplicity.In all simulations the relaxation time is set to τ c ≡ 1.

B. Multicomponent lattice Boltzmann
There are different extensions for the lattice Boltzmann method to simulate multi-component and multiphase systems [34][35][36][37][38].An overview on different methods for multi-component fluid systems and the treatment of fluid-fluid interfaces is given in Ref. [39].In this paper, the method introduced by Shan and Chen is used [34].Every species has its own distribution function following Eq.(1).To obtain an interaction between the different components a force is calculated locally and is included in the equilibrium distribution function.it is summed up over the different fluid species c and x , the nearest neighbors of lattice positions x. g cc is the coupling constant between the species and Ψ c (x, t) is a monotonous weight function representing an effective mass.For the results presented here, the form is used.To incorporate F c (x, t) in f eq i we define The macroscopic velocity included in f eq i is shifted by ∆u c as As we are interested in immiscible fluids we choose a positive value for g cc which leads to a repulsive interaction.This interaction has to be strong enough to obtain two separate phases but it should not be too high in order to keep the simulation stable.Here, we use the range of 0.08 ≤ g cc ≤ 0.14.

C. Nanoparticles
Particles are simulated with molecular dynamics where Newton's equations of motion F = m upar and D = J ωpar (10) are solved by a leap frog integrator.F and D are the force and torque acting on the particle with mass m and moment of inertia J. u par and ω par are the velocity and the rotation vector of the particle.
The particles are also discretized on the lattice.They are coupled to both fluid species by a modified bounceback boundary condition which was originally introduced by Ladd [9,[40][41][42][43].This changes the lattice Boltzmann equation as follows: where c i is the velocity vector pointing to the next neighbor.C depends linearly on the local particle velocity, ī is defined in a way that c i = −cī is fulfilled.A change of the fluid momentum due to a particle leads to a change of the particle momentum in order to keep the total momentum conserved: If the particle moves, some lattice nodes become free and others become occupied.The fluid on the newly occupied nodes is deleted and its momentum is transferred to the particle as A newly freed node (located at x) is filled with the average density of the N FN neighboring fluid lattice nodes x iFN for each component c, Hydrodynamics leads to a lubrication force between the particles.This force is reproduced automatically by the simulation for sufficiently large particle separations.If the distance between the particles is so small that no free lattice point exists between them this reproduction fails.
If the smallest distance between two identical spheres with radius R is smaller than a critical value ∆ c = 2 3 the correction term is given as [43]: µ is the dynamic viscosity, rij a unit vector pointing from one particle center to the other one and u i is the velocity of particle i.To use this potential for ellipsoidal particles Eq. ( 15) is generalized in a way proposed by Berne and Pechukas [8,44,45].We define σ = 2R and = 3πµ 8 σ.Both are extended to the anisotropic case as with and ôi the orientation unit vector of particle i. R and R ⊥ are the parallel and the orthogonal radius of the ellipsoid.Using Eq. ( 16) we can rewrite Eq. ( 15) and obtain F is a dimensionless function taking the specific form of the force into account and in this example it is . The lubrication force (including the correction) already reduces the probability that the particles come closely together and overlap.For the few cases where the particles still would overlap we introduce the direct potential between the particles which is assumed to be a hard core potential.To approximate the hard core potential we use the Hertz potential [46] which has the following shape for two identical spheres with radius R: r is the distance between particle centers.For larger distances φ H vanishes.K H is a force constant and is chosen to be K H = 100 for all simulations.To use this potential for ellipsoidal particles Eq. ( 18) is generalized in a similar way as the lubrication force.Using Eq. ( 16), σ = 2R and = K H σ 5 2 we can rewrite Eq. ( 18) and obtain φH is a dimensionless function taking the specific form of the potential into account and in this example it is φH (x) = (1 − x) 5/2 .The Shan-Chen forces also act between a node in the outer shell of a particle and its neighboring node outside of the particle.This would lead to an increase of the fluid density around the particle.Therefore, the nodes in the outer shell of the particle are filled with a virtual fluid corresponding to the average of the value in the neighboring free nodes for each fluid component: ρ c virt (x, t) = ρ c (x, t).This can be used to control the wettability properties of the particle surface for the special case of two fluid species which will be named red and blue.We define the parameter ∆ρ and call it particle color.For positive values of ∆ρ we add it to the red fluid component: For negative values we add its absolute value to the blue component: In Ref. [8] it is shown that there is a linear relation between ∆ρ and the three-phase contact angle θ.

III. EMULSIONS
In this section the different types of particle stabilized emulsions and the effect of the particle shape on some of their properties are discussed.We find two different types of emulsions in our simulations, namely the Pickering emulsion (Fig. 1, left) and the bijel (Fig. 1, right).The choice of parameters (such as particle contact angle, particle concentration, fluid-fluid ratio, particle aspect ratio) determines the type of emulsions.Parameter studies for emulsions have been discussed in Refs.[9] and [8] for spherical and ellipsoidal particles, respectively.In the current publication we limit ourselves to anisotropy effects on the time dependence of the emulsion formation.We use the following particle shapes (m = R /R ⊥ is the particle aspect ratio.R and R ⊥ are the parallel and orthogonal radius of the particles, respectively): prolate ellipsoids (m = 2; Fig. 1, top), spheres (m = 1; Fig. 1, center) and oblate ellipsoids (m = 1/2; Fig. 1, bottom).For m = 1/2 we choose R = 5∆x and R ⊥ = 10∆x.For the other values of m the radii R and R ⊥ are chosen as such that the particle volume is kept constant, resulting in R ≈ 12.6∆x and R ⊥ ≈ 6.3∆x for m = 2 as well as R = R ⊥ ≈ 7.9∆x for spheres.The interaction parameter between the fluids (see Eq. ( 6)) is chosen as g br = 0.08 which corresponds to a fluid-fluid interfacial tension of σ = 0.0138.The particles are neutrally wetting (contact angle θ = 90 • ) and the particle volume concentration is chosen as C = 0.24.The simulated systems of volume S have periodic boundary conditions in all three directions and a side length of L S = 256∆x = 32R .Initially, the particles are distributed randomly.At each lattice node a random value for each fluid component is chosen so that the designed fluid-fluid ratio is kept (1:1 for the bijels and 5:2 for the Pickering emulsions).When the simulation evolves in time, the fluids separate and droplets/domains with a majority of red or blue fluid form.
The average size of droplets/domains L(t) can be determined by measuring Here, is the average domain size in direction i.
is the fluctuation of φ which is the Fourier transform of the order parameter field ϕ = ρ r − ρ b .In this publication, the time is given in simulation timesteps, which can be converted to physical units.We use Eq. ( 4) and Eq. ( 5) to relate the kinematic viscosity to ∆x and ∆t.By assuming ν = 10 −6 m 2 /s, the kinematic viscosity for water, R = 125nm and R = 7.9∆x (this is the value used for the spherical particle, see above) we fix the chosen resolution  22)) for m = 1 and m = 2.At first view, a steady state is reached after about 10 5 timesteps.L(t) is larger for bijels than for Pickering emulsions, which is due to the measurement being based on the Fourier transform of the order parameter.Ellipsoids are able to stabilize larger interface areas than spheres leading to smaller L(t).
of the simulation.Thus, we obtain ∆x = 15.8nm and ∆t = 4.2 × 10 −11 s and a total system size of L S ≈ 4µ.The interfacial tension is then σ = 3.14 × 10 −08 N/m.Larger system sizes can be reached with the same computational effort by compromising on the resolution.
The time development of L(t) for the three different particle types (prolate, spherical and oblate (m = 2, 1 and 1/2) and for Pickering emulsions and bijels is shown in Fig. 2. We can identify three regimes: in the first few hundred timesteps the initial formation of the droplets/domains starts.Then, the growth of droplets/domains is being driven by Ostwald ripening.At even later times, droplets/domains grow due to coalescence.When two droplets unify, the area coverage fraction of the particles at the interface is increased because the surface area of the new droplet is smaller than that of the two smaller droplets before.At some point the area coverage fraction of the particles is sufficiently high to prevent further coalescence.The state which is reached at that time is (at least kinetically) stabilized and one obtains a stable emulsion.The values for L(t) are larger for bijels than for Pickering emulsions.This can be explained by the way we calculate L(t) (see Eq. ( 22) and related text) using a Fourier transformation of the order parameter field.
It can clearly be seen that anisotropic particles are more efficient in interface stabilization than spheres since they can cover larger interfacial areas leading to smaller fluid domains (note that the simulation volume is kept constant).However, the difference in L(t) for m = 2 and m = 1/2 is small.This can be understood as follows: if a neutrally wetting prolate ellipsoid is adsorbed at a flat interface, it occupies an area A P,F (m > 1) = m 1/3 A p,s , where A p,s is the occupied interface area for a sphere with the same volume.This corresponds in the case of m = 2 to the occupied interface being larger by a factor of 1.26 as compared to spheres.For an oblate ellipsoid the occupied interface area is A P,F (m < 1) = m −2/3 A p,s which for m = 1/2 is by a factor of 1.59 larger than the area occupied by spheres.Since in emulsions the interfaces are generally not flat, these formulae can only provide a qualitative explanation of the behavior of L(t): If the interface curvature is not neglectable anymore, we loose some of the efficiency of interface stabilization, which is more pronounced for m < 1.This explains why the value of L(t) for m = 1/2 is only slightly smaller than for m = 2.
It seems that L(t) reaches a steady state after some 10 5 timesteps for both types of emulsions and for all three values of m.However, if one zooms in one can observe that L(t) develops for a longer time period if the particles have a non-spherical shape.As will be demonstrated below, the reason for this phenomenon is the additional rotational degrees of freedom due to the particle anisotropy.Furthermore, the time development of L(t) for emulsions stabilized by prolate particles requires more time than that for the oblate ones.If a particle changes its orientation as compared to the interface or a neighboring particle this generally changes the interface shape.In this way the domain sizes are influenced, leading to changes of L(t) -an effect which is not observed for m = 1.Fig. 3 and Fig. 4 depict a zoom-in of the time development of L(t) for Pickering emulsions and bijels with m = 2 and m = 1/2, respectively.One observes that L(t) decays in all four cases.The kink in Fig. 3 after about 2.8 million timesteps is due to the coalescence of two droplets of the Pickering emulsion.A substantial difference is the range of the decay.It is larger for the bijel since it consists of a single large interface whereas the Pickering emulsion consists of many small interfaces.The large interface in the bijel is much more deformable.This explains the larger range of the decay of L(t) for the bijel.The fluctuations are of the same order for Pickering emulsions and bijels.Furthermore, the range of the decay is larger for m = 2 than for m = 1/2.The time of reordering is much shorter for m = 1/2 as compared to m = 2.These effects can be explained by the presence of additional rotational degrees of freedom for the anisotropic particles.While oblate particles have only a single additional rotational degree of freedom as compared to spheres, prolate particles show an even more complex behavior due to their second additional rotational degree of freedom.
In this section we demonstrated that particle anisotropy causes additional timescales to influence the growth of domains in particle-stabilized emulsions.In the following sections we discuss model systems in order to obtain a deeper understanding of this effect.We will restrict ourselves to prolate particles with m = 2. Furthermore, the high resolution of the particles in the current section was only chosen to be able to sufficiently resolve the oblate objects.In order to reduce the required computational resources, we use smaller particles in the model systems studied below (R = 8∆x and R ⊥ = 4∆x).It has been checked carefully that the reduced particle size does not have a qualitative impact on the results.

IV. SINGLE PARTICLE ADSORPTION
In the previous section we demonstrated that there is an additional time development of the average domain size L for emulsions stabilized by anisotropic particles.In the following sections we relate this behavior to the orientational degree of freedom of the particles at the interface.To obtain a more basic understanding of the additional timescales some simple model systems are dis- ), m = 2 and σ ≈ 0.041.A particle is placed as such that it just touches the undeformed interface.The dashed lines denote the adsorption trajectories, the solid lines show the points where the particle touches the undeformed interface.The circular points depict the stable and the metastable point.The square points are related to the snapshots describing the adsorption process in the inset.For initial particle orientations of ϑ(t = 0) = 0 • the particle ends in its stable configuration orientated parallel to the interface.
cussed.The simplest possible example is the adsorption of a single particle at a flat fluid-fluid interface.To characterize the particle orientation we introduce the angles ϑ and φ. ϑ is the angle between the particle main axis and the y-axis, where the y-axis is oriented perpendicular to the flat fluid-fluid interface.φ is the angle between the particle main axis and the x-axis, where the x-axis is orientated parallel to the interface.ξ is the distance between the particle center and the undeformed interface in units of the long particle axis.In this section we consider the case of neutral wetting (θ = 90 • ) and restrict ourselves to an aspect ratio of m = R /R ⊥ = 2.The fluid-fluid interaction parameter is set to g br = 0.1 corresponding to an interfacial tension of σ ≈ 0.041.We use a cubic system with 64 lattice nodes in each direction.A wall is placed at the top and bottom in y-direction.Periodic boundary conditions are applied in the x-and z-direction.In order to obtain a flat interface the system is filled with two equally sized cuboid shaped lamellae with an interface orthogonal to the y-axis.The lamellae are mainly filled with red and blue fluid, respectively.The initial majority and minority species are set to ρ maj = 0.7 and ρ min = 0.04.For this study a particle is placed so that it just touches the (undeformed) fluid-fluid interface.This is done for different initial orientations of the particle.The inset of Fig. 5 shows snapshots of a typical adsorption process.In the beginning the particle is oriented almost orthogonally to the interface.In the first ca.2000 timesteps the particle moves towards the interface without changing its orientation considerably.Then, the particle rotates and reaches its final orientation after 3600 timesteps.The outer plot of Fig. 5 shows a ϑ − ξ diagram of the adsorption.The points where the particle just touches a flat interface for the different orientations are marked with solid lines.The dotted lines indicate the adsorption trajectories.Each black square is related to one of the snapshots in the inset of Fig. 5. Almost all dashed lines end in the upper circle which corresponds to the equilibrium point where the free energy function has a global minimum.Just the cases with an initial value of ϑ(t = 0) = 0 • end at the metastable point at ϑ = 0 • as shown by the circle at the bottom of Fig. 5.This metastable point might not be found in experiments: on the one hand fluctuations will cause a rotation of the particle towards the stable points and on the other hand, it is impossible to place the particle exactly at ϑ = 0 • .Fig. 6 and 7 depict the dynamics of the particle adsorption and the influence of the initial particle orientation ϑ with respect to the flat interface.Fig. 6 shows the time development of ϑ for different values of ϑ(t = 0).For ϑ(t = 0) = 0 • and ϑ(t = 0) = 90 • (upper and lower lines) the orientation remains unchanged and the adsorption at the interface causes only a translational particle movement.The lines for the three other simulation runs start at ϑ(t = 0) = 22.5 • , ϑ(t = 0) = 45 • and ϑ(t = 0) = 67.5 • .All of them go in the 'wrong' direction during the first few 10 2 timesteps and end at ϑ = 90 • corresponding to the stable point, but the time needed for reaching this value differs.Furthermore, in all cases during the first timesteps, ϑ decreases but then it increases up to this final value.The time t e the particle needs to reach the final orientation of ϑ = 90 • depending on ϑ(t = 0) is shown in the outer plot of Fig. 7. Due to the discretization of the particle on the lattice, its orientation shows small deviations from the theoretical final value.Therefore, we measure t e as the time when the angle reaches 98% of the theoretical final angle.The particle oscillates arround this final value but these oscillations are very small and their magnitude falls below the threshold for the measurement of t e .t e increases with decreasing ϑ and diverges for ϑ → 0. This divergence can be understood using the inset of Fig. 7.If the starting angle ϑ(t = 0) comes closer to ϑ = 0 • (corresponding to the metastable case where the particle never flips) the capillary forces causing the particle rotation become smaller and vanish.
We have seen that anisotropy of particles causes additional timescales in the development of the domain sizes in the emulsions, because of orientational ordering.This timescale is of the order of 10 6 LB timesteps.In this section we have shown that the adsorption of a single particle at an interface and its orientational ordering takes of the order of 10 3 timesteps and depending on the initial particle orientation towards the interface.We can identify one extra timescale where the particles rotate towards the interface.This timescale plays a role in the beginning of the emulsion formation (during droplet formation and droplet growth) when the particles come in contact with the interfaces.However, this timescale does not yet explain the full time development.We require additional model systems to obtain a full understanding of the additional timescales.Thus, we consider many particles at a flat interface as well as at a single droplet in the following sections.

V. PARTICLE ENSEMBLES AT A FLAT INTERFACE
After having studied the adsorption of a single particle we discuss the behavior of a many-particle ensemble at a flat interface.What is the influence of the hydrodynamic interaction between many particles on the timescales involved in emulsion formation?For the case of the singleparticle adsorption the particle orientation towards the interface (ϑ) is an important parameter.For prolate particles, also the mutual orientation (φ) of the particles is important and one has an additional degree of freedom leading to particle orientational ordering.To characterize the ordering of the particles we use two order parameters and two correlation functions.
Measures for global ordering effects of the particles are the orientational order parameters S and Q.We define the uniaxial order parameter S [47,48] as where denotes the averaging over particles.Originally S is an order parameter for studying liquid crystals which indicates the phase transition from the isotropic to the anisotropic/nematic phase.Here, the parameter S is used as a measure for the orientation of the particle ensemble towards the interface.If all particles are oriented orthogonal to the interface we have S = S ⊥ = 1 (see top right of Fig. 8).The orientation of all particles parallel to the interface leads to S = S = −0.5 (see top left of Fig. 8).
The biaxial order parameter Q [48] is defined as The parameter Q is a measure for the mutual orientation of the particles oriented parallel to the interface.If all particles lying parallel to the interface are oriented in the same direction it is Q = Q aniso = 1.5.Q = Q iso = 0 means that the particles oriented parallel to the interface have a two-dimensional isotropic ordering.
The local ordering effects are investigated by using two correlation functions.The discretized form of the pair correlation function g(r) is defined as where N is the number of particles, r and r ij are the distance from a reference particle and the distance between the two particle centers of particle i and j in units of R , respectively, and g n is a normalization factor chosen such that g(r) → 1 for r → ∞. g(r) gives a probability to find a particle at a distance r from a reference particle.It is a measure for the ordering of the particle centers and ignores the orientation.As a measure for the local orientational ordering effects the angular correlation function [49] is defined as (in the discrete form) with l = 1 in order to have the appropriate values of h for a given value of ϑ discussed below.h(r) gives a measure for the average orientation of particles at distance r from a reference particle.If the particles at distance r from the reference particle are all oriented parallel to the reference particle we have h(r) = 1 (see right and left configuration in the bottom of Fig. 8) and an orthogonal orientation leads to h(r) = −1 (see central configuration in the bottom of Fig. 8).In the following we use smoothed versions of g and h, where we average over neighboring data points.The flat interface considered in this section is periodic in two dimensions parallel to the interface and each period has a size of A I = L 2 I , with L I = 512 = 64R .The system is confined by walls 40 lattice units distant from the interface in the third dimension.The particle coverage fraction for N particles adsorbed at the interface is defined as χ(ξ, ϑ) = N A P (ξ,ϑ) A I .
A P (ξ, ϑ) is the area which the particle would occupy on a hypothetical flat interface and depends on the distance between the particle center and the undeformed interface and the particle orientation relative to the flat interface and ξ is the distance between particle center and undeformed interface.In the following we relate the coverage fraction to the case of ξ = 0 and ϑ = 90 • (χ I ) or ϑ = 90 • (χ F ) corresponding to the initial state and the equilibrium state for θ = 90 • (see previous section).This leads to χ I = N A P,I A I and χ F = N A P,F A I with A P,I = πR 2 ⊥ and A P,F = πR R ⊥ .Initially, the particles are oriented almost orthogonally to the interface (see top right of Fig. 8).The initial value for the polar angle is chosen as ϑ ≈ 0.6 • for all particles, whereas φ and the particle positions are chosen randomly.Analogously to the case of the single-particle adsorption the particle flips to an orientation parallel to the interface (see Fig. 9).Fig. 10(a) shows the time development of S for different values of χ I (χ I ≈ 0.08 (squares), χ I ≈ 0.38 (circles), χ I ≈ 0.46 (upward pointing triangles) and χ I ≈ 0.52 (downward triangles)) and the time development of Q for χ I ≈ 0.38 (diamonds).The parameter Q starts at 0 and ends at a small value (Q final ≈ 0.05 Q aniso ) far away from the value of total ordering.A similar behavior is found for all values of χ I .Fig. 9 shows that there are smaller domains where particles are oriented in the same direction.But every domain has a different preferred particle direction which might lead to small but still finite values of Q.Another reason for this effect is the finite system size and finite particle number which change the parameter as follows [49]: Q ∞ is the value of the biaxial order parameter that the corresponding system with an infinite amount of particles would have.The parameter S starts for all values of χ I with a value of S ⊥ = 1, corresponding to the initial configuration.For lower values of χ I the parameter S reaches S = −0.5, corresponding to the case that the particles flip completely.For higher values of χ I the final value of the parameter is S final > S = −0.5.This corresponds to the case where some particles cannot flip completely to the equilibrium orientation because there is insufficient space.The final values of S (obtained after 10 5 timesteps) are shown in the outer plot of Fig. 10(b) as a function of χ I .We find a transition point at χ I,C ≈ 0.42 corresponding to χ F,C ≈ 0.84.If all particles are oriented parallel to the interface the system corresponds practically to a two dimensional system of ellipses.However, the value of χ I,C found is below the value of the closest packing density for a two-dimensional system of ellipses with m = 2, which is χ F,max ≈ 0.91.Such a system was also studied in Ref. [49] with Monte Carlo simulations.For the case of an ellipse with an aspect ratio m 2d = 2 a transition point of χ 2dmc ≈ 0.78 from isotropy to a solid phase was found.The solid phase de-  24) and Eq. ( 25)) for m = 2, θ = 90 • , σ ≈ 0.041.Q(t) is shown for a single value of χI only since it stays at a value of approximately 0 for all χI .S(t) is shown for different values of χI .In case of highly packed interfaces, i.e. for large values of χI , not all particles are able to fully align with the interface.For larger values of χI S needs a longer time to get into the equilibrium than shown here.(b) Outer plot: The final values of the order parameter S are plotted for different particle densities χI .As shown in 10(a) a transition from a fully ordered to a disordered state can be found at a critical value of χI,C ≈ 0.42.Inset: the time the order parameter S (defined in Eq. 24) requires to reach the final value (time which particles need to flip).For small values of χI t f is independent of χI but above a critical value of χI = χI,C t f increases with increasing χ by almost one order of magnitude.
scribes a state where the particle centers as well as the orientations are ordered.We do not reach the limit of the solid phase.This suggests that hydrodynamic interactions and absence confinement in the third dimension still play a dominant role.The biaxial order parameter in the MC system grows up to Q ≈ 1 (see Fig. 11 in Ref. [49]) corresponding to a global anisotropic state with a quite high degree of ordering for χ F > χ 2dmc .This effect is not observed in our system.The reason for this difference is the method used to reach this state.A twodimensional system of ellipses was studied in Ref. [49] wheres we simulated three-dimensional ellipsoids which form an effective two-dimensional system by flipping to the interface.We can see that in the many-particle system and for small and moderate χ I about 10 3 timesteps are required for the particles to flip which is the same order of magnitude as in the case of the single particle adsorption for small values of χ I .The inset in Fig. 10(b) shows the time the order parameter S needs to reach its final value.This corresponds to the time required for the whole particle ensemble to be flipped completely (χ I < χ c ) or to reach the semi-flipped state for χ I > χ c .For χ I < 0.38 t f stays almost constant at about 4500 timesteps.In this regime the distance between the particles is sufficient so that the influence of hydrodynamic interactions on the flipping behavior can be neglected.For higher values it increases very sharply and hydrodynamic interactions between the particles must not be neglected anymore.Furthermore, the time needed to flip completely for the very dense systems (jammed state) is about one order of magnitude larger.
The biaxial order parameter does not show any global ordering but the snapshot in Fig. 9 shows some local ordering effects.Hence, we need other ways to characterize the local ordering effects and utilize the two local correlation functions g(r) and h(r) defined above.The particles have a contact angle of 90 • , so there are no capillary interactions between them in the final state when all of them have flipped completely and the system has reached an equilibrium.However, there are dipolar interface deformations and thus the interactions during the flipping process of the particles and for χ > χ which causes capillary interactions at this time.After flipping there are still some capillary waves going through the system, leading to interactions between the particles.The pair correlation function g(r) is shown in Fig. 11(a) for three different values of χ I (χ I ≈ 0.23, χ I ≈ 0.31 and χ I ≈ 0.38) after 10 5 timesteps.The first peak is pronounced in all three cases.The distance r of this peak decreases for increasing χ I as well as the degree of ordering.For the highest χ I a depletion region leading to a minimum after the peak is pronounced.To obtain a measure of the local orientational ordering effects we investigate the orientational correlation function h(r) as shown in Fig. 11(b) for the same 3 values of χ I .The first two positive peaks and the first negative peak can be explained with the drawings in the bottom of Fig. 8.The first positive peak is due to a side-to-side alignment of two particles.Fig. 9 shows several domains of side-to-side alignment.The first negative peak comes from an alignment where the particles are oriented perpendicular to each other and the second positive peak comes from a tip-to-tip alignment or second nearest neighbors of side-to-side orientation.The degree of translational and orientational ordering increases with  26)) (b) orientation correlation function h(r) (defined in Eq. ( 27)).In both cases the ordering increases with increasing χI .
increasing χ I .After having discussed the correlation functions we investigate the time development of g(r) in order to understand the time development of the average domain size L(t).Fig. 12 shows g(r) at different times between 10 4 and 10 6 timesteps.The first peak decreases but at later times the following peaks are more pronounced.Thus, the degree of ordering increases.After 4 • 10 5 timesteps this development has almost come to an end.The reason for this remaining development is the particle reordering.The particles form domains where they align parallel to each other.These domains become larger with time.
In this section we have shown shows that the presence of many particles at an interface leads to two additional timescales in the reordering.The first one is the rotation of the particle towards the interface.The particle rotates towards its final orientation parallel to the interface.For lower values of χ I this process does not depend on χ I and is not different from the single particle adsorption.For larger values of χ I the time needed to come to its final orientation increases.Hydrodynamic as well as excluded volume effects become more important.Above a critical value not every particle reaches its 'final' orientation.
The reordering of h (corresponding to g in Fig. 12) can also be observed.The first 2 peaks get more pronounced after several 10 5 timesteps as compared to the state after 10 4 timesteps shown in Fig. 11(b).

VI. PARTICLE ENSEMBLES AT A SPHERICAL INTERFACE
In the previous chapter the behavior of particle ensembles at a flat interface was discussed.However, in emulsions the interfaces are generally not flat.Pickering emulsions usually have (approximately) spherical droplets and a bijel has an even more complicated structure of curved interface.The simplest realization of a curved interface is a single droplet and as such is studied in this section.
The simulated system is periodic and each period has a size of L S = 256 lattice units.The droplet radius and the number of adsorbed particles are chosen to be R D = 0.6L S ≈ 76.8 and 600, respectively.In the beginning of the simulation the particles are placed orthogonal to the local interface tangential plane.As we have seen already for the case of flat interfaces the particles flip to an orientation parallel to this tangential plane.This state is shown in Fig. 13 after 2 • 10 5 timesteps.A preliminary comparison between flat and spherical interfaces has already been given in our previous contribution [39].The time development of S is shown in Fig. 11(a) in Ref. [39].It has been found that the influence of the interface curvature on the flipping process is larger than the influence of the particle coverage.The time needed for the particles to flip is about a factor two smaller in the case of the Here, we investigate the particle correlation function (see Eq. ( 26)) for the particle ensemble.Fig. 14 shows g for χ I ≈ 0.27 at three different times.After 10 4 timesteps it is still close to the correlation function of the initial condition.After 10 5 timesteps some changes can be seen.The first peak is reduced but the second peak is more pronounced.There is no substantial change between 1•10 5 and 2•10 5 timesteps.Compared to the state at 10 5 timesteps the correlation function shows pronounced peaks at longer distances from the particle (about 6R p ).The particles mostly reorder during the first 10 5 timesteps since at later times only minor changes in the particle order can be observed.Similar to the case of flat interfaces that was discussed in the previous section, the particle ensemble forms domains where the particles are ordered in a nematic fashion.The peaks in the correlation function are more pronounced in the case of droplets than in the case of a flat interface.The reason is given by the capillary interactions between the particles which are much stronger in the case of curved interfaces.In particular, non-zero capillary interactions persist between spheroids even in the case of neutrally wetting particles.
The time development of g at the droplet as discussed in this section differs from the behavior in the case of a flat interface.For the droplet, g arrives at its final structure after about 10 5 timesteps whereas at the flat interface about four times more as many steps are required.In addition, for flat interfaces, g only shows one or two peaks (depending on χ I ), while for the particle covered droplet five peaks are found due to a larger range of ordering of the particles.This is a result of the stronger capillary interactions between the particles due to the interface curvature.
We can understand one of the additional timescales with the behavior of the ellipsoidal particles at a single droplet.The particles reorder and it can be shown that this leads to a small deviation of the shape of the droplet which is (almost) exactly spherical in the beginning [50].A change of the interface shape caused by reordering of anisotropic particles leads to a change of L(t).The reordering of particle ensembles at flat as well as spherical interfaces takes of the order of 10 5 timesteps.This reordering takes place in idealized systems with constant interfaces which do not change their shape considerably.In real emulsions, however, the interface geometry changes substantially during their formation.For example, two droplets of a Pickering emulsion can coalesce.After this unification the particle ordering starts a new.This explains the fact that the additional timescale we find in our emulsions is of the order of several 10 6 timesteps.

VII. CONCLUSION
In this article we have investigated the dynamics of the formation of Pickering emulsions and bijels stabilized by ellipsoidal particles.In contrast to emulsions stabilized by spherical particles, spheroids cause the average time dependent droplet or domain size to slowly decrease even after very long simulation times corresponding to several million simulation timesteps.The additional timescales related to this effect have been investigated by detailed studies of simple model systems.At first, the adsorption of single ellipsoidal particles was shown to happen on a comparably short timescale (≈ 10 4 timesteps).Second, many particle ensembles at flat interfaces, however, might require substantially more time in case of sufficiently densely packed interfaces.Here, local reordering effects induced by hydrodynamic interactions and interface rearrangements prevent the system from attaining a steady state and add a further timescale to the emulsion formation (≈ 10 5 timesteps).Third, this reordering is pronounced in the case of curved interfaces, where the movement of the particles leads to interface deformations and capillary interactions.During the formation of an emulsion, droplets might coalesce (Pickering emulsions) or domains might merge (bijels).After such an event the particles at the interface have to rearrange in order to adhere to the new interface structure.Due to this, the local reordering is practically being "restarted" leading to an overall increase of the interfacial area on a timescale of at least several 10 6 timesteps.With the nanoscale resolution chosen above, this corresponds to physical times of the order of 10 −5 s.
Our findings provide relevant insight in the dynamics of emulsion formation which is generally difficult to investigate experimentally due to the required high temporal resolution of the measurement method and limited optical transparency of the experimental system.It is well known that in general particle-stabilized emulsions are not thermodynamically stable and therefore the involved fluids will always phase separate -even if this might take several months.Anisotropic particles, however, provide properties which might allow the generation of emulsions that are stable on substantially longer timescales.This is due to the continuous reordering of the particles at liquid interfaces which leads to an increase in interfacial area and as such counteracts the thermodynamically driven reduction of interface area.

2 FIG. 2 .
FIG.2.Pickering emulsion and bijel: Time development of the average domain size L(t) (see Eq. (22)) for m = 1 and m = 2.At first view, a steady state is reached after about 10 5 timesteps.L(t) is larger for bijels than for Pickering emulsions, which is due to the measurement being based on the Fourier transform of the order parameter.Ellipsoids are able to stabilize larger interface areas than spheres leading to smaller L(t).

2 FIG. 3 .
FIG. 3. Pickering emulsion: zoom of the time dependent average domain size L(t) for m = 2, m = 1 and m = 1/2.The slow but continuous decrease of L(t) clearly shows the occurrence of additional timescales in the domain growth.The kink in the measurement for m = 2 can be adhered to the coalescence of two droplets.

2 FIG. 4 .
FIG.4.Bijel: Zoom for m = 2 and m = 1/2: Time development of the average domain sizes depicting the impact of the additional timescales.The range of the variation of L(t) is larger as compared to the Pickering emulsion due to the impact of a small deformation on the larger effective interface of the bijel.

FIG. 5 .
FIG.5.Outer plot: ϑ-ξ-plot for neutral wetting (θ = 90 • ), m = 2 and σ ≈ 0.041.A particle is placed as such that it just touches the undeformed interface.The dashed lines denote the adsorption trajectories, the solid lines show the points where the particle touches the undeformed interface.The circular points depict the stable and the metastable point.The square points are related to the snapshots describing the adsorption process in the inset.For initial particle orientations of ϑ(t = 0) = 0 • the particle ends in its stable configuration orientated parallel to the interface.

FIG. 6 .
FIG.6.Time development of the particle orientation ϑ(t) for different initial orientations.For ϑ(t = 0) = 0 • and ϑ(t = 0) = 90 • the particle rotates in the 'wrong' direction in the first timesteps.The time needed to be in the final orientation depends on the initial orientation.

FIG. 7 .
FIG.7.Outer plot: Time te which the particle needs to reach the final orientation (ϑ = 90 • ) for different initial orientation angles ϑ0 = ϑ(t = 0) from ϑ = 0 • to ϑ = 90 • .te diverges if ϑ0 approaches 0 • .The reason for the divergence is the approach of ϑ0 to the orientation of the metastable point, as it is shown in the inset: If the starting angle (middle dashed line) approaches ϑ(t = 0) = 0 (lower dashed line) the time required to reach the equilibrium point diverges.

FIG. 8 .FIG. 9 .
FIG. 8. Top right: Sketch of the initial condition for the many-particle system.All particles are oriented almost orthogonal to the interface corresponding to the initial configuration.Top left: Sketch of the final state.All particles are oriented parallel to the interface.Bottom: Different constellations of mutual orientation of next neighbors.

FIG. 10 .
FIG.10.(a) Time development of the two order parameters S(t) and Q(t) (see Eq.(24) and Eq.(25)) for m = 2, θ = 90 • , σ ≈ 0.041.Q(t) is shown for a single value of χI only since it stays at a value of approximately 0 for all χI .S(t) is shown for different values of χI .In case of highly packed interfaces, i.e. for large values of χI , not all particles are able to fully align with the interface.For larger values of χI S needs a longer time to get into the equilibrium than shown here.(b) Outer plot: The final values of the order parameter S are plotted for different particle densities χI .As shown in 10(a) a transition from a fully ordered to a disordered state can be found at a critical value of χI,C ≈ 0.42.Inset: the time the order parameter S (defined in Eq. 24) requires to reach the final value (time which particles need to flip).For small values of χI t f is independent of χI but above a critical value of χI = χI,C t f increases with increasing χ by almost one order of magnitude.

FIG. 12 .
FIG. 12. Time development of g(r) for χI ≈ 0.38.The second peak is more pronounced at later timesteps.The particles reorder and the ordering increases.The reordering process is almost done after 4 • 10 5 timesteps.