Electrohydrodynamic assembly of colloidal particles on a drop interface

1.   Introduction

Emulsions of particle-coated drops (so called Pickering emulsions ^[1,2]) are widely used in the pharmaceutical, food, personal care and many other industries ^[3]. Colloidosomes, which are microcapsules with shells made of colloid particles ^[4], are used for drug delivery due to the great degree of control of the shell permeability ^[5]. Colloidal particles get trapped at interfaces between immiscible (e.g., oil/water) fluids ^[6] (since the energy to detach a particle adsorbed at an interface exceeds the thermal energy by the thousands) and at high packing density form a shell enapsulating the drop that stabilizes emulsions against coalescence and enables selective permeability ^[4,7,8,9,10]. Colloids at low surface coverage, however, do not form static structures, but instead assemble dynamically ^{[11,12,13,14]}. For example, a uniform electric field was found to induce various patterns such as an equatorial belt, pole-to-pole chains or a band of dynamic vortices ^[11,12]. The latter intriguing phenomenon has not been explained thus far and motivates our study.

The collective dynamics of colloidal particles adsorbed at the interface of a drop in the presence of an electric field presents a challenging problem. The presence of the interface strongly modifies the electrostatic and hydrodynamic interactions between the particles, and introduces new interactions such as the classic capillary attraction (Cheerios effect)^[15]. The latter arises from local deformation of the interface (e.g., if the particle has weight) ^{[16,17,18,19,20]}. The interfacial distortions increase the interfacial area and thus raise the interfacial free energy; one way of minimizing this effect is to bring the particles together so that their menisci overlap. The electrostatic interaction between charged particles in the absence of electric field ^[21] or in the presence of electric field ^[22,23,24] can significantly differ from the particle interactions in a homogeneous medium. Our previous work ^[25] showed that depending on the particles and suspending fluid conductivities, the far-field interaction between interface-trapped particles may get significant contribution from the electric-field induced particle quadrupole and overcome the dipole-dipole interaction. Particle hydrodynamics at interfaces is also highly nontrivial as shown by studies of the motion of an isolated particle ^{[26,27,28,29,30,31]} or simulations of many particles ^[32,33,34].

In previous work ^[25] we used the effective moment method to study the dynamics of particles in a homogeneous fluid but in an applied non-uniform electric field. Here we generalize this approach to determine the collective motion and assembly of colloidal spheres trapped at a drop interface in a uniform applied electric field. Our aim is to use this model to simulate the dynamics seen in recent experiments on particle covered drops ^[11,12] and to systematically show the impact of particle concentration and field strength on the dynamics. The model presented here accounts for the electric field driven flow within the drop and suspending fluid, particle-particle electrostatic interaction, and the particle motion and rotation due to the induced flow and the applied electric field. The impact of particle concentration and electric field strength on the collective motions of the particles is investigated. We also present simulations illustrating the effect of changing particle coverage and applied field strength to the observed clustering phenomenon near the equator of the drop.

2.   Problem formulation

A drop placed in an electric field polarizes if its permittivity ${\epsilon}_ d$ and conductivity $\sigma_ d$ are different than the suspending fluid permittivity ${\epsilon}_ f$ and conductivity $\sigma_ f$. Finite conductivity, even if very low, enables the passage of electric current and electrical charge accumulates at the drop interface. The electric field acting on this induced surface charge creates shear electric stress that drag the fluids into motion. In the case of a uniform DC electric field of strength $E_0$, the fluid undergoes axisymmetric straining flow about the drop. For a spherical drop of radius $R_0$, the surface velocity is ^[35]

where $\theta$ is the spherical polar angle measured from the direction of the applied uniform electric field (see Figure 1). Here $\mu_f$ and $\mu_d$ are the viscosity of the suspending fluid and drop, respectively. This flow field is directed either from pole to equator or equator to pole depending on the relative magnitude of the drop and suspending fluid charge relaxation time, $\frac{\sigma_d}{\epsilon_d}$ and $\frac{\sigma_f}{\epsilon_f}$.

Particles adsorbed at the drop interface are advected by the flow and if the flow is from pole to equator, i.e., $\frac{\sigma_d}{\epsilon_d} < \frac{\sigma_f}{\epsilon_f}$, particles are expected to accumulate at the drop equator. Hence, the band formation observed in the experiments ^[11,12] is likely driven by the electrohydrodynamic flow. However, particles also polarize in the electric field and interact electrostatically, which typically results in chaining ^{[11,12,36,37]}.

To study the complex interplay of flow advection and dielectrophoretic motion on the collective dynamics of particles, we develop a model under several assumptions. Drop shape remains spherical, i.e., the electric capillary number $Ca = {\epsilon}_f E_0^2R_0/\gamma\ll 1$, where $\gamma$ the surface tension. We only consider particle motions tangential to the drop surface, i.e., in the normal direction all forces exerted on a particle are balanced. This is a reasonable assumption given the large surface energy needed to remove a particle from the interface, $\sim \gamma \pi a^2\sin^2\Theta$, where $\Theta$ is the contact angle at the interface of the three-phase line between the two fluids and the particle, and $a$ is the particle radius. Since the particles are small, interfacial deformation due to particle weight is negligible (Bond number $Bo = a^2\Delta \rho g/\gamma \ll1$, where $\Delta \rho$ is the density difference, and $g$ the acceleration due to gravity). The problem geometry is sketched in Figure 1.

Hereafter, all variables are nondimensionalized using the radius of the drop $R_0$, the field strength $E_0$, a characteristic applied stress $\tau_c = {\epsilon}_ f E_0^2$, and the properties of the suspending fluid. Accordingly, the time scale for the electrohydrodynamic flow is $t_ {\mathrm {ehd}} = \mu_ f/\tau_c$. The flow velocity scale is $u_c = R_0 \tau_c/\mu_ f$.

We utilize the model of $N$ neutral particles subjected to an applied electric field in a homogeneous medium developed by us ^[25]. The evolution of the induced electric dipole and quadrupole of particle $i$ ($i = 1, ..., N$) are

and

where

Here $sym$ denotes $A_{ij}^{sym} = A_{ij}+A_{ij}^\mathsf{T}$ and superscript $T$ denotes transpose. The interparticle distance vector from particle $i$ to particle $j$ is defined as ${{\bf R}}_{ij} = {\bf x}_i- {\bf x}_j$ where ${\bf x}_i$ is the position of the $i^{th}$ particle. The vector ${\bf{\Omega}}_i$ is the rotation rate of particle $i$. The constants in the above equation are defined by

and

The interface at which the particles are trapped separates fluids with different properties. We assume that the particles are instead in a homogeneous medium with effective properties $\mu_m = (\mu_ d+\mu_ f)/2$, $\sigma_m = (\sigma_ d+\sigma_ f)/2$ and ${\epsilon}_m = ({\epsilon}_ d+ {\epsilon}_ f)/2$. The parameters with the subscripts $cm$ are sometimes referred to as the Clausius-Mossotti factors while the parameters with the subscripts $mw$ are the Maxwell-Wagner times which measure in the absence of rotation the times for the dipole and quadrupole moments to relax toward a steady state. Note here the evolution of higher order moments can be obtained in a similar way. This analysis assumed that the ratio $\gamma = a/d$ of the particle radius $a$ to mean inter-particle distance $d$ is small and of the same order as the ratio $d/R_0$ of the mean inter-particle distance to the drop radius $R_0$. In our earlier work ^[25] we used this assumption to ensure that the accuracy of the above electrostatic force model for the particle dynamics entered at the same order of accuracy as the hydrodynamic model below. This requires the assumption that the length scale over which changes in the local electric field strength varies along the particle is given by $R_0$. The latter assumption is certainly reasonable given the the local electric field varies from being tangent to the drop interface near the equator to normal to the interface at the poles of the drop.

Particle $i$ translates with velocity ${\bf u}_i$ and rotates with rate ${{\bf{\Omega}}}_i$ in response to the electrostatic forces, the global flows, and the flows generated by neighboring particles ^[25].

where $\widehat{ {\bf R}}_{ij}$ is the unit vector in the ${{\bf R}}_{ij}$ direction and from Eq (2.1), ${ {\bf u}}^\infty_i = B \sin(2 \theta_i) \vec{t}_\theta$ is the electrohydrodynamic flow induced about the drop. The rotation vector is given by

with

The dimensionless electric torque is

and ${{\bf F}}^{rep}_{ij}$ is the assumed interparticle replusive force. Our earlier work ^[25] assumed a polynomial in the separation distance between particles with a sharp cut off for large separations. Here we have chosen a continuous form of the repulsive force which is easier to program and accomplishes the same task. The non-dimensional coefficients in the above equations are

Here $E_R$ is a reference applied electric field.

Particle rotation is induced by the flow and electric field. In a uniform applied electric field, particles in a homogeneous fluid can rotate above a threshold electric field due to the Quincke effect ^[38]. The critical field strength for an isolated particle as been calculated for both uniform ^[39,40] and nonuniform electric field strengths ^[25]. For a uniform field it is given by

Our aim is to simulate the experimental results of Ouriemi and Vlahovska ^[12,41] and to make more general statements about the effects of particle number and electric field on a partially coated drop. To this end, we focus on polyethylene particles on a silicon oil drop in castor oil. With this in mind, we set $\epsilon_p = 2.25 \epsilon_0$, $\epsilon_d = 2.8 \epsilon_0$, and $\epsilon_f = 4.7 \epsilon_0$, where $\epsilon_0$ is the permittivity of free space. In addition we set $\mu_d = 0.05$ $Pa \, s$, $\mu_f = 0.69$ $Pa \, s$, $\sigma_d = 3.6 \times 10^{-12} \, S/m$, $\sigma_f = 3.8 \times 10^{-11} \, S/m$, and we set $\sigma_p \sim 0$. This allows us to set $\sigma_{cm} = -1/2$ and $\sigma^\prime_{cm} = -1/3$. If more conductive metals (silver or aluminum) are considered we set $\sigma_{cm} = 1$ and $\sigma^\prime_{cm} = 1/2$. In addition, we have set $E_R = 200$ $V/mm$, $a = 50 \mu$m, $R_0 = 2.5 mm$, and $\beta_R = 1$. We have introduced the damping factor $\chi$ into the definition of the parameter $B$. One can think of this as representing the effect of the applied flow on the particle dynamics for our complex problem beyond what can be expected for a single particle in a homogeneous flow. In our calculations we have set $\chi = 1/20$, more will be said later about this choice.

Since the particle motions are restricted to a spherical surface, we need to finally project the translational vector velocity in the tangential direction. Using spherical coordinates $\theta$ (polar angle measured from the $z$ axis) and $\phi$ (azimuthal angle), the trajectory equation (in dimensionless coordinates) is given by,

where $\hat{\theta}$ and $\hat{\phi}$ are unit vectors in the $\theta$ and $\phi$ directions, respectively. This completes the description of our model.

3.   Simulation results

The model proposed in the previous section consists of a system of $2N$ ordinary differential equations for the particle positions along the spherical drop interface, $\theta_i$, $\phi_i$, $i = 1, ..., N$, the $3N$ ODE's for the components for each particle dipole moment $\bf{P}_i$, and $5N$ ODE's for the independent components of the quadrupole moment $\bf{Q}_i$, giving a total of $10N$ ODE's to solve. Note that $\bf{Q}_i$ is symmetric with zero trace. In addition, the rotation vector $\bf{\Omega}_i$ and the electric torque $\bf{\Gamma}_i$ must be calculated for each particle. Each ODE at particle $i$ includes a contribution from each of the other particles. Our solution approach was to solve this system using the Runge-Kutta-Fehlberg (RKF) adaptive time stepping method on a GPU parallel computing system. This proved to be a very efficient approach which allowed for shared memory, and the simultaneous solution of the ODEs on each of the GPU's cores. The error tolerance was set to $10^{-4}$ with an initial time step of $10^{-4}$.

The low conductive particle case will be considered in Sections 3.1-3.3. In this case Qunicke motion is possible since $\epsilon_{cm} - \sigma_{cm} = -0.154 + 0.5 = 0.346 > 0$. In Section 3.4 we consider the high particle conductivity case where $\epsilon_{cm} - \sigma_{cm} = -0.154 -1.0 = -1.154$ and by Eq 2.9 there is no critical electric field strength for Quincke motion.

3.1.   Introductory examples

All parameter values are given in the above sections except for the strength of the applied electric field $E_0$ and $N$ the number of particles. Solutions as a function of each of these parameters will be given below. $E_0$ will be given in dimensional units of $Volts/mm$. Initially the particles are randomly distributed on the surface of the spherical drop. An example of an initial distribution is given in Figure 2a. Different initial data give slight differences in the overall quantitive dynamics but not in the qualitative dynamics, e.g., how a solution with a given set of parameter values will look at an instant in time. Besides the particle positions the initial values of the dipole and quadrupole moments were set to zero. All runs were ran out to ${t}_f = 15000$ dimensionless time units. Other times were tested but this dimensionless time appeared to give reasonable and consistent predictions for the range of parameters considered here. Note that the scaling implies that the dimensional time is proportional to ${t}$ divided by the applied electric field strength squared, i.e., $E_0^2$. Hence results presented below for $E_0 = 400 V/mm$ are presented at a dimensional time four times earlier than results for $E_0 = 200 V/mm$. Our choice of ${t}_f$ was taken so that the results presented for the largest applied fields appeared to have stabilized into a quasi-steady state (there are still random fluctuations) at the time the results are presented.

Figure 2b shows the particle positions on the drop surface at time ${t}_f$ for $N = 512$ particles and an applied field of $E_0 = 200 V/mm$. Since the strength of the applied field is greater than the critical field for Quincke rotation, i.e., $E_0 > E_c$, we expect the particles to rotate. The sign of the scalar $\zeta = {{\bf{\Omega}}}_i\cdot {\bf \hat n}_i$ indicates the rotation of each particle relative to the normal of the interface. To represent this rotation we color each particle either red ($\zeta > 0$) for a clockwise rotation or blue ($\zeta < 0$) for a counter-clockwise rotation. Note that the particles have evolved into a oscillatory belt like structure around the equator of the drop. Although particles with different rotations are scattered about the band, there appears to be a clustering of particles with the same rotation about the normal to the drop interface. A time evolution of the initial data evolves in 3 stages: (1) the applied flow field forces all the particles to move towards the equator and form a belt like structure. The initial particle rotation is determined by Eq 2.6. Since the initial dipole and quadrupole moments are assumed to be zero, the initial rotation in this model depends on the random initial location of each particle and the inter-particle repulsive force given in Eq 2.7. Because of this the particle rotations appear to be randomly distributed around the equator; (2) once the uniform belt is formed there is an oscillation of the belt about the equator with a random distribution of particle rotations; (3) the particles begin to cluster into regions with similar rotations as illustrated in Figure 2b. This is a very dynamic belt which is continuously evolving both on the particle scale with particles moving between clusters and changing the sign of $\zeta$, plus evolving on the scale of the drop radius with cluster moving around the diameter of the drop. Different initial data could result in a different specific picture of the particles at a given instant in time, but the qualitative dynamics over long times is similar. The amplitude and period of oscillation is clearly the most striking visual observation for the figure. How the strength of the field and particle number impact these quantities will be investigated below.

Figure 2c shows the particle positions on the drop surface at time ${t}_f$ for $N = 512$ particles and an applied field of $E_0 = 100 V/mm$. The strength of the applied field is below the critical field for Quincke rotation. Because of this the particles appear to form periodic and line structures near the equator in their steady state. The dynamic pictures show some particle motion with particles occasionally moving between neighboring strictures, but the overall shape is retained in time. The particles still have a colored rotation because we apply the color due to the sign of $\zeta$ and not its magnitude, but it is small and can be ignored. Notice that where clusters of particles exist there is a hexagonal like packing. The width and the interparticle spacing is dependent of the applied electric field strength $E_0$ which directly controls the strength of the applied flow field, and the number of particles $N$. It might be expected that a stronger applied flow field would compress the particles, even if there was significant rotation.

In our model the strength of the applied flow field is controlled by the damping factor $\chi$. We have set $\chi = 1/20$ in all reported results because with this value our computational results come closest to reproducing the experimental results reported in the work of Ouriemi and Vlahovska ^[12,41]. To illustrate this in Figure 3 we plot the particle positions on the drop surface at time ${t}_f$ for $N = 512$ particles and an applied field of $E_0 = 200 V/mm$, but with a large damping factor of $\chi = 1$. Notice that unlike the results in Figure 2b for $\chi = 1/20$, the particles are now strongly confined to a belt about the equator of width approximately 3 particles in depth. This belt is in a (quasi-) steady state with little particle motion. The particles rotate but because of the strength of the flow there is little spacial motion. This steady shape with a hexagonal particle packing is similar to the low applied field strength case presented in Figure 2c but unlike that case the particles are forced into a more narrow width belt. Although we admit that the choice of $\chi$ is somewhat arbitrary, chosen to reproduce as best as possible the predictions of a given set of experiments, we do claim that with this choice of damping parameter, the overall qualitative behavior of a partially coated drop with $E_0$ and $N$ is determined. Additional research is needed to relate $\chi$ to the experimental parameters.

3.2.   Impact of particle coverage and electric field

Examples of changing the number of particles $N$ along the interface are given in Figure 4a for $N = 256$ and in Figure 4b for $N = 1024$ both with $E_0 = 200 V/mm$. In Figure 4a we find that decreasing the particle density allows for clusters of particles with similar rotations to occur. This is similar to the experimental observations of Ouriemi and Vlahovska ^[12,41] for low particle coverage. Comparing to the $N = 512$ case shown in Figure 2b we might have expected a continuous belt but the decrease in particle number and the fact that we have not changed the applied force from the flow field, since it depends in $E_0$, has allowed for clusters of the same rotation, and height similar to those in Figure 2b to exist. Increasing the particle number to $N = 1024$ appears to have two effects. The particles still appear in a continuous band as in the $N = 512$ case, but the wavelength of the band oscillation appears to have increased, and clustering of the particles along the sides of the bands is more noticeable. The distribution of particles now appears to be evenly distributed. Both of these cases are very dynamic, with the clusters in the $N = 256$ under constant rotation, and the band in the $N = 1024$ case doing an oscillatory motion.

As the number of particles continues to increase, the integrity of the belt structure degrades. An example is given in Figure 4c for $N = 2048$ and $E_0 = 200 V/mm$. At this field strength, oscillations and particle circulation was observed for the lower values of $N$, but now, because of the large number of particles a steady state structure can be observed in regions. As the particles appear to collect in a region the hexagonal like structure appears in spots but fault lines are observed where there is a slight rotation of the basic structure. These fault lines are expected to be caused by the curvature of the drop surface since a uniform packing of the particles is not allowed as we increase from the equator to the poles. We also see the thinning of the belt in regions with an active transport of particles. This is just a snapshot in time but what is clear is that a steady state oscillatory like belt structure is not observed. It is possible that such a structure could exist if the calculation were continued further in time, but we have not observed it. Because of this, we limit any general conclusion below to the behavior of the particle belts to cases where $N$ is less than or equal to $1024$. A more improved model may be necessary to make general conclusions about higher density partially coated drops.

Increasing the electric field strength can be expected to decrease the amplitude of the oscillation because it increases the applied flow strength. This is illustrated in Figure 5 where $N = 512$ and $E_0 = 400 V/mm$ and should be compared to the results in Figure 2b for $E_0 = 200 V/mm$. There also appears to be a slight change in the particle belt period. These results are typical for our model with other values of particle number $N$. Care needs to be taken in increasing the applied electric field strength in our model. It has been observed experimentally that as the applied field strength increases, the drop becomes unstable and can deviate significantly from the steady spherical shape assumed here. Examples can be found in the experimental results of Ouriemi and Vlahovska ^[12,41] where both steady shape deformations and wobbling drops are illustrated. In the Ouriemi and Vlahovska study of silicon oil drops in castor oil, $\gamma = 4.5 mN/m$. They considered partially coated drops for $0 < Ca < 5$, where for the small values of $Ca$ they observed steady state spherical drop shapes, while for slightly larger values the drops became oblate and prolate steady shapes up to $Ca \approx 1.5$. For larger values of $Ca$ the drops became unstable and started to wobble. In our analysis, using the values for silicon oil drops in castor oil we have that for $E_0$ between $100 V/mm$ and $400 V/mm$, that $0.23 < Ca < 3.70$. Hence while we can expect the spherical drop approximation to be good for the $E_0 = 100 V/mm$ ($Ca = 0.23$) or $E_0 = 200 V/mm$ ($Ca = 0.92$) cases discussed above, for $E_0 = 400 V/mm$ ($Ca = 3.70$) the drop is probably unstable and we present these results primarily to illustrate what to expect with the particle dynamics as the applied electric field is increased.

3.3.   Impact of external field strength

The above results give us a qualitative picture of how varying the applied electric field and the number of particles impact the dynamics of the belt of particles along the equator of the drop. To get a more quantitative behavior we examine next how both applied electric field and particle number impact the belt width and the number of peaks (the wavelength) along the equatorial belt.

The amplitude of the oscillation $A_\theta$ can be simply defined as the angle between the maximum distance of all particles along the $z-$direction minus the minimum $z$ value. Clearly the answer will depend on the time selected for this measurement. Our approach was to average this measurement over an interval of saved time snapshots, approximately one unit in time. Since the Runge-Kutta-Fehlberg (RKF) is an adaptive time stepping method with the time step variable, the results (a snapshot) were usually saved at times slightly larger than one unit, and hence there were less that 15000 particle time snapshots saved per run. To compute $A_\theta$ the first ${t}_s = 5000$ snapshots were ignored, and the rest up to ${t}_f$ were included in the averaging. Increasing ${t_s}$ did have some effect on the average $A_\theta$ but the variation was within the error presented in the plots below. An average over a minimum of 20 runs were done for each $E_0$ and $N$ value for different random initial data. The error bars in the plots below are a result of these different cases.

In Figure 6 we plot $A_\theta$ as a function of $E_0$ for both $N = 512$ and $N = 1024$. The plots were made using the Matlab function boxplot with the central mark indicating the median, and the bottom and top edges of the box indicating the $25th$ and $75th$ percentiles, respectively. Although some variation exists in the data, the trend is obvious. As noted earlier, $E_0 = 100 V/mm$ is below the Quincke threshold $E_c$ and we get a tight packing of the particles about the equator. As $E_0$ is increased above $E_c$ the belt begins to oscillate and $A_\theta$ suddenly increases. Further increase in $E_0$ decreases the amplitude of oscillation. This is expected since the the magnitude of the applied flow field increases with $E_0$. We also find that increasing $N$ increases $A_\theta$, which may be expected since more particles are trying to collect about the equator.

The wavelength of the belt oscillation is another quantity of interest. Although we could visually estimate this value from the snapshot of the particle distribution at time ${t}_f$, we decided to find an average value in a manner similar to how the average $A_\theta$ was calculated. The same data sets and times were analyzed as for $A_\theta$. At each time snapshot, a curve for the average shape of the belt was found as follows. The drop was divided into 80 sections in the $\phi$ direction of size $2 \pi / 80$. In each section the average particle height was found. If no particles were within a section, the average was set to $0$. These points were Fourier transformed and the wave number with the maximum amplitude was identified, $A_w$. As in the calculation of $A_\theta$ this was an average value determined for all time snapshots greater than $5000$, and these values were averaged over 20 runs with random initial data. The results as a function of $E_0$ and for $N = 512$ and $1024$ as plotted in Figure 7. The low value of $E_0 = 100 V/mm$ was not plotted since it formed an approximately uniform belt.

Figure 7 suggests that the wavelength of oscillation is insensitive to $E_0$ but it does depend on $N$. In particular, a smaller value of $N$ gives a higher wavenumber, or smaller wavelength. This is consistent with the visual results presented in Figures 2 and 4.

3.4.   Very conductive particles

An example of a highly conductive particle case is presented in Figure 8 at time ${t}_f$ for $N = 512$ and $E_0 = 200 V/mm$. In this case Qunicke motion is not predicted and we do not observe it. We do find particle chaining along the equation. This is an expected phenomena since now the polarization of each particle is parallel to the direction of the applied electric field. Depending on the applied field strength increasing particle density either lengthens the chain length or adds to the number of chains along the equator. There is some observed particle dynamics at these long times with chains growing or shortening. Some of this dynamics is due to numerical noise in the calculations, and some is due to the quasi-stability of a particular chain.

4.   Conclusions

We developed a model to simulate the collective dynamics of colloids trapped at the interface of a drop placed in a uniform DC electric field. We assumed that the particles were constrained to move on a sphere, representing the drop, and replaced the two-phase drop-suspending fluid system by an effective medium with properties that were the average of the drop and suspending fluid. The model accounted for the electric field driven flow within the drop and suspending fluid, particle-particle electrostatic interaction, and the particle motion and rotation due to the induced flow and the applied electric field. The electrostatic interactions are computed by approximating the particles by dipoles and quadrupoles. The hydrodynamic interactions are accounted by reflections accurate to $O(d/R_0)^4$ in the ratio between the interparticle separation $d$ and the drop radius $R_0$ ^[25]. Despite our modeling simplifications, the model captures the experimentally observed particle assemblies such as chains, bands and dynamic vortices around the equator. Specifically we found that the model predicts the formation of chains in the case of conducting particles and an undulating band around the equator in the case of dielectric particles. In the case of non-conducting particles, we showed that in the presence of Quincke rotation, the amplitude of the undulations of the observed equatorial particle belt increases with particle concentration but decreases with electric field strength, and that the wavelength of the undulations appears independent of the applied field strength. Our simulations also show that with increasing particle coverage or applied field strength isolated particle clusters become continuous belts with possible rotating clusters within the belts. Given that the strength of the hydrodynamic Taylor flow increases with electric field strength, the belt formation is not unexpected as the field strength is increased.

Although the model seems to qualitatively capture the observed dynamics of the experiments, improvements are needed for quantitative predictions. For example, a better modeling of the electrohydrodynamic force interaction between particles which accounts for the change in both fluid and electric properties across the fluid interface could impact specific predictions. Allowing for drop deformation due to the applied field could have a significant impact on the predictions with increasing electric field strength. Also missing here is the impact of the three phase contact line at the particle-fluid-fluid intersection. Future work will focus on these and other improvements in the modeling.

Acknowledgements

This work was supported in part by NSF grant DMS-1716114, and CBET-1704996.

Conflict of interest

The authors declare there is no conflicts of interest.