Mesoscopic electrohydrodynamic simulations of binary colloidal suspensions

A model is presented for the solution of electrokinetic phenomena of colloidal suspensions in fluid mixtures. We solve the discrete Boltzmann equation with a BGK collision operator using the lattice Boltzmann method to simulate binary fluid flows. Solvent-solvent and solvent-solute interactions are implemented using a pseudopotential model. The Nernst-Planck equation, describing the kinetics of dissolved ion species, is solved using a finite difference discretization based on the link-flux method. The colloids are resolved on the lattice and coupled to the hydrodynamics and electrokinetics through appropriate boundary conditions. We present the first full integration of these three elements. The model is validated by comparing with known analytic solutions of ionic distributions at fluid interfaces, dielectric droplet deformations and the electrophoretic mobility of colloidal suspensions. Its possibilities are explored by considering various physical systems, such as breakup of charged and neutral droplets and colloidal dynamics at either planar or spherical fluid interfaces.


I. INTRODUCTION
Electrokinetic phenomena play a fundamental role in both technological and natural systems, from micro-and nanofluidic devices to molecular biological processes [1][2][3]. Their physical description has proven challenging mainly due to the wide range of relevant length-scales involved, from interfacial effects, molecular in origin and thus typically on the nanoscale, to the colloidal and system sizes, usually several orders of magnitude larger. An additional significant challenge is to capture the competing effects of two long-range interactions, namely hydrodynamic and electrostatic. These same factors also complicate the realization of experiments as it is, for example, challenging to resolve the charge distribution at solid-fluid or fluid-fluid interfaces, especially at the nanoscale [4]. In this context, numerical simulations become an attractive alternative to explore those limits where theory and experiments struggle.
Simulations of electrohydrodynamic phenomena vary significantly in methodology, also as a consequence of the wide range of relevant scales. The most common models directly solve the Taylor-Melcher leaky dielectric model, based on the assumption that charge is confined at interfaces in boundaries of negligible length, that is, that the Debye length is small compared to other relevant lengthscales [5,6]. These are, essentially, macroscopic models and simulations, which ignore the kinetics of the ions and volumetric ionic concentrations. On the opposite side of the spectrum, nanoscale electrokinetic problems, increasing in popularity together with the development of nanofluidics, makes full Molecular Dynamics (MD) simulations a viable alternative for studying electrokinetic phenomena [7,8]. Between these two methods lies a variety of mesoscopic models which permit the simulation of larger time and length-scales than MD, while also resolving the kinetics and spacial distribution of ions. Hydrodynamics can be resolved by mesoscale particlebased methods, such as multi-particle collision dynamics, drastically reducing the number of particles needed compared to full-MD studies [9][10][11][12]. Another common approach involves the solution of Navier-Stokes equations using standard CFD techniques, coupled to Molecular Dynamics for resolving either charged colloids and/or individual ions [13]. Full mesoscopic models, which allow the simulation of even longer timescales and a wider range of salt concentrations, treat both the solvents and the ions-through ionic concentrations-at the continuum level. These are usually based on the Nernst-Planck advection-diffusion equation as a model for ion transport, which requires the solution of the Poisson equation to determine the electric field, comprising a system of equations usually referred to as Poisson-Nernst-Planck theory (PNP) [14].
In the following we present a numerical solution of an electrohydrodynamic mesoscopic model of colloidal suspensions in binary fluid mixtures with dissolved ions. We focus on nanoscale flows of negligible inertia and comparable electronic, hydrodynamic and colloidal scales. The model considers the ions (solutes) and the fluids (solvents) at the macroscopic level, described by continuous fields of ion concentrations, mass densities and velocity, respectively. The finite-sized colloids are coupled to the fluids and the ions through proper boundary conditions. The simulation methodology recovers the hydrodynamics by means of the lattice Boltzmann method (LBM) [15], and the ion kinetics via a finite difference discretization of the Nernst-Planck equation, inspired by the link-flux method [14]. Previous LB-based electrohydrodynamic models for binary mixtures have used a freeenergy functional [16] to derive the ion kinetics and the forces that the ions exert on the fluids, as well as the interactions between the different fluids [14,17,18]. Here we present a different approach, deriving the diffusive ion fluxes via the forces exerted on them, and recovering the coupling forces from solvent-solvent and solventsolute microscopic interactions using pseudopotentials of the form proposed by Shan and Chen [19,20]. The colloidal particles' coupling with the binary mixture is implemented using the Ladd methodology [21], extended to binary mixtures [22,23]. Particles interact between themselves via a Hertz potential when in contact, as also via a lubrication force [21,22]. Finally, particle-ion interactions are resolved using a partial-volume discretization that reduces discretization errors, as recently proposed by Kuron et al. [24]. Overall this presents the first description of a model capable of handling the electrokinetics of both binary fluid mixtures and moving particles.
The model is validated by comparing with known theoretical results in several systems encompassing all the method's functionalities. We obtain excellent agreement for dielectric droplet deformation due to Maxwell stresses and distribution of ions at fluid/fluid interfaces. The electrophoretic mobility of colloidal suspensions is also in agreement with previous simulations and experiments. We also show exemplary cases of the possibilities of the model, such as transient dynamics and breakage of droplets in various conditions, colloid wetting properties at fluid/fluid interfaces and the dynamics of colloidcoated droplets.
The paper is structured as follows: in Section II we first describe the hydrodynamic (II A) and electrokinetic (II B) equations, then the interactions between the different species (II C), the Poisson's equation (II D) and finally the colloidal dynamics (II E), each followed by a description of its numerical solution. Section III is dedicated to the study of several systems of increasing complexity, for validation of the model and presentation of the simulation capabilities. Finally conclusions and possible future directions are discussed in Section IV.

A. Hydrodynamics
The lattice Boltzmann method is used to solve two Boltzmann equations with a Bhatnagar-Gross-Krook (BGK) collisional operator [25], where f σ (r r r, ξ ξ ξ σ , t) is the probability distribution function for component σ ∈ (a, b) with particle velocity ξ ξ ξ σ ; τ is the relaxation time (the timescale for collisions to drive the local distributions to thermodynamic equilibrium), taken to be equal for both fluids; Φ Φ Φ σ B the external force andf σ (r r r, ξ ξ ξ σ , t) the local equilibrium distribution function. The macroscopic density and the average macroscopic momentum are recovered via ρ σ = f σ dξ ξ ξ and ρ σ u u u σ = ξ ξ ξ σ f σ dξ ξ ξ, with u u u σ the component macroscopic velocity. When taking the equilibrium distribution as where v v v σ = ξ ξ ξ σ − u u u, with u u u the barycentric velocity and p the pressure, it is known that, following the Chapman-Enskog procedure with the Knudsen number as small parameter, the Navier-Stokes equations are recovered, namely Here ⊗ is the outer product, p the pressure, I the identity matrix, and s s s the deviatoric stress tensor, s s s σ = λ σ (∇ ·u u u)I + η σ (∇u u u + ∇u u u T ), with η the dynamic viscosity and λ the bulk viscosity. The forcing term F F F σ has two main contributions, coming from interactions between components and electrostatic forces, F F F σ ≡ F F F σ I +F F F σ E ; both of these terms are derived further down. As the nanoscale is much smaller than the capillary length, gravity is disregarded.
The LBM consists in the discretization of Eq. (1) in space, time and velocities, such that at a given timestep t and at each site of a regular Cartesian lattice r r r i , the distribution function f σ (r r r, u u u, t) becomes f σ d (r r r i , t), with d the index of the discrete velocity vectors c c c d [26]. Here, we use the usual D3Q19 lattice [15,26] with spacing ∆x, such that particles can travel at each time-step ∆t to the nearest and next-nearest neighbors. Thus Eq. (1) can be written as Notice that we have not included the Φ Φ Φ σ B term, as external forces will be included as perturbations of the velocity in the equilibrium distribution functionf σ , following Shan and Chen [19]. After discretization and expansion to second order, Eq. (2) can be written as The equilibrium velocity of each component is given by The lattice weights for a D3Q19 lattice are B. Electrokinetics The cations and anions dispersed throughout the solvents are considered at the continuum level. The evolution of the concentration of each ion species is given by the advection-diffusion Nernst-Planck equation, where n ± (r r r) is the number density of cations ( + ) and anions ( − ). The diffusive ion flux is with the diffusivity D (which we have assumed to be homogeneous and the same for both species), and β ≡ 1/k B T , with the Boltzmann constant k B and the temperature T . Notice that this temperature only serves as an energy scale, as we have considered isothermal systems with no fluctuations. In all our studies we take β = 1, although for clarity we keep the symbol β where present. The Nernst-Planck equation ignores ion-ion interactions. These are known to be relevant at the nanoscale for moderate ion concentrations, where steric and electric interactions can alter the flow of ions [27][28][29][30]. Therefore our model is strictly valid for low ionic concentrations n ± (r r r) 1/∆x 3 . The force density applied to the ions, F F F ± (r r r), is coupled to the flux through the Smoluchowski relation with mobility Dβ. This implies that the inertial timescale of the ions is much faster than any other resolved timescale.
Therefore ions instantly reach their drift velocity as F F F ± (r r r) varies, their dynamics being dominated by the solvents' viscosity. The force has two contributions, electrostatic and a term coming from the microscopic interactions between solvents and solutes, with q ± (r r r) = z ± en ± (r r r), z ± the valence of cations and ions and −e the electron charge. In simulations we take e as the unit of charge, e = 1, although for clarity we keep the symbol in the expressions. Furthermore, we only consider monovalent salts, z ± = ±1. The total electric potential φ(r r r) is the sum of an internal and external contribution φ(r r r) = φ I (r r r) + φ E (r r r), with the external one usually being a parameter used to control the flow. Analogously, we define equivalent electric fields as E E E(r r r) = −∇φ(r r r) = E E E I (r r r) + E E E E (r r r). The internal contribution is itself a function of the charge distribution via Poisson's equation, where we have defined the total charge density q(r r r) ≡ e(z + n + (r r r) + z − n − (r r r)) + q p (r r r). Here q p (r r r) is the total charge from suspended particles, as will be detailed further below. The permittivity ε(r r r) = ε 0 ε r (r r r), with ε 0 the vacuum permittivity and ε r (r r r) the relative permittivity of the solution, is not homogeneous, as each solvent and the colloidal particles can have different permittivities. At each point, ε(r r r) is either a function of the fluid concentration, or the permittivity of the particles, such that ε(r r r) = ε c (r r r) r r r inside the particle, ε f (r r r) r r r inside the fluid.
As the dependency of the fluid's permittivity on the concentration of each fluid component is a priori unknown, we follow Rotenberg et al. [17] in taking the simplest model, a linear interpolation between the value of each fluid, ε a and ε b , as a function of the local concentration, ε f (r r r) = 1 2 ε a (1 − c(r r r)) + ε b (1 + c(r r r)) , with c(r r r) ≡ (ρ b (r r r) − ρ a (r r r))/(ρ a (r r r) + ρ b (r r r)). The Nernst-Planck equation (8) is numerically solved using a finite difference scheme with different discretization methods for the advective and diffusive terms, inspired by the link-flux method [14,17]. Integrating Eq. (8) over a lattice site cell's volume, which we assume to be always a cube of side length ∆x, results in n ± (r r r i , t + ∆t) − n ± (r r r i , t) ∆t = 1 ∆x × d∈{Q6} −u u u(r r r i )n ± (r r r i+d/2 ) + j j j ± (r r r i+d/2 ) · c c c d . (14) where we have defined r r r i+d/2 ≡ r r r i +c c c d ∆t/2. Here we take d to run only over the D3Q7 lattice directions, as it simplifies the discretization and improves performance. Note that in the D3Q7 lattice c c c d is equivalent to the normal of the surface of the lattice cell in the d direction, which is always a square of area ∆x 2 . Previous research using the link-flux method has always considered a D3Q19 lattice, both to have consistent discretizations for solutes and solvents and to reduce the spurious diffusion produced by advection. Using the methodology detailed in Ref. [14], we have measured that in our case the effective diffusivity stays well below a 2% difference with the nominal diffusivity D, for the comparably low value used here D = 0.01, and the largest possible fluid velocities. Previous researchers have performed careful comparisons of the accuracy of different differential operators, revealing that higher lattice connectivities indeed increase the accuracy and convergence (as a function of the resolution), although total errors are expected to be negligible in our cases of variations on the order of the diffusive interface [31].
For the advection terms we use a first-order upwind discretization, such that d∈{Q6} n ± (r r r i+d/2 )u u u(r r r i ) · c c c d = d∈{Q6} H[u u u(r r r i ) · c c c d ]n ± (r r r i )u u u(r r r i ) · c c c d , (15) with H[x] the Heaviside function. For the diffusive flux we follow a different procedure than in Refs. [14,17], as it has been recently shown that their proposed exponentiation of the fluxes leads to a quadratic error for flows far from equilibrium [32]. We thus use a straightforward finite-difference discretization for the gradients, and a linear interpolation to obtain the value of the fields at the volumes' surfaces. This has been shown to be more computationally efficient while also reducing spurious fluxes [32]. Denoting r r r i+d ≡ r r r i +c c c d ∆t, Eq. (9) results in j j j ± (r r r i+d/2 ) = D n ± (r r r i+d ) − n ± (r r r i ) ∆x c c c d −β q ± (r r r i ) + q ± (r r r i+d ) 2 The form and discretization of the last term, the interaction force, is given in what follows.

C. Couplings
Microscopic interactions between the different solvents and solutes are modeled using the pseudopotential method of Shan and Chen [19], which is among the most popular multiphase/multicomponent LB methods [33]. It is usually preferred to other alternatives due to its simplicity of implementation. The method allows us to determine F F F ± I (r r r) and F F F σ I (r r r) in a consistent manner. Microscopic molecular interactions between species (either solvents or solutes) are captured as local force densities of the form (17) with the indexes α, β ∈ {a, b, +, −}, α = β. As we have neglected ion/ion interactions, we take the coupling constants G +− = G −+ = 0 Let us first focus on fluid/fluid interactions, α, β ∈ {a, b}. The coupling constant G ≡ G ab = G ba determines the strength of attraction (G < 0) or repulsion (G > 0) between the solvents. The pseudopotential ψ σ (r r r i , t) is taken to have the common form The reference density is set to ρ 0 = 1, although we leave the symbol for clarity. In the continuum limit, and to fourth order in derivatives, the Shan-Chen force (17) can be shown to be [34] F This results in a modified equation of state [19] with the speed of sound c s = ∆x/ √ 3∆t. It can thus be seen that, depending on the value of G, Eq. (20) leads to phase separation. In this study we always take G = 4.0∆t −2 ρ −1 0 . Next, we consider the case of fluid/ion microscopic interactions, referred to as solvation. These are the cases α ∈ {a, b}, β ∈ {+, −} in Eq. (17). For the ionic pseudopotentials we take the same form as for the fluids, although as n ± (r r r) 1/∆x 3 they can be simplified, ψ ± (r r r i ) = n ± 0 (1 − exp(−n ± (r r r i )/n ± 0 )) ≈ n ± (r r r i ). (21) with the reference density n ± 0 = 1/∆x 3 . It follows that, in the continuum limit, the forces applied to the fluids from ion interactions take the form We have here assumed, for simplicity, that G b± = 0. This is without loss of generality, as the coupling coefficients can also be negative, and thus hold the possibility of modelling both attraction and repulsive forces from solvent a, the latter equivalent in our scheme to an attractive force to solvent b. Analogously, the forces applied to the ions from their microscopic interaction with the fluids take the form In order to specify the value of G a± , we interpret our model of ion/solute microscopic interactions as an approximate model of solvation. Solvation interactions are captured at the macroscopic level by chemical potentials, µ ± s (r r r), corresponding to the cost in free energy of adding an ion to the solvent or mixture. For example, in the most common case of hydration (solvation in water), µ ± s (r r r) is fundamentally a function of the hydrogen bonds formed between the ions and the solvent molecules [35,36]. Selective solvation in mixtures, when there is a difference in the solvation energies between the two solvents and solutes, can have crucial effects on the global dynamics, as usually the associated energies far exceed k B T [36]. Even though the solvation energies are known to depend on several parameters, here we follow a similar approach as Onuki et al. [36] and take a simple form of a linear dependency with the density of solvent a, The , commonly referred to as Gibbs transfer energy [17,35], is the difference of the respective solvation chemical potentials at the bulk of each component in a phase-separated system. For example, in the case of a planar interface at x = 0 between two immiscible fluids,ρ a a = ρ a (−∞), and ρ a b = ρ a (∞). Analogously ∆ρ a ≡ρ a a −ρ a b . Variations of chemical potentials translate to forces, given by the Gibbs-Duhem relation (at constant temperature) f f f = −∇p = ρ∇µ ρ . Therefore, variations in the solvation chemical potential generate a force of the form −(∆µ ± /∆ρ a )n ± (r r r)∇ρ a (r r r). Comparing with Eq. (23), and noticing that for small densities ψ a ≈ ρ a , it is evident that in order to interpret the pseudopotential interactions as solvation effects, the Shan-Chen force coupling parameter has to be taken as Having determined the forces stemming from microscopic molecular interactions between solvents and solutes, two other forces rest to be considered. The first comes from the relative movement of the ions with the respect to the solvent, and is given by a simple friction coupling between the two, that is The second remaining force has its origin on the dielectric nature of the solvents. Polarization effects give rise to the Kelvin force density, which assuming an isotropic permittivity takes the form F F F E (r r r) = − 1 2 E(r r r) 2 ∇ε + 1 2 ∇(E(r r r) 2 φ∂ε/∂φ) [37]. Using Eq. (13) gives with the average permittivityε = 1 2 (ε a + ε b ). The Kelvin force density is strictly valid in the electroquasistatic limit and for incompressible media [38]. The latter condition is fulfilled considering flows of low Mach number, Ma = u 0 /c s 1, with u 0 the typical flow speed. The electroquasistatic limit is satisfied for small ion fluxes, that is, D ∆x 2 /∆t. In what follows we always take D = 10 −2 ∆x 2 /∆t.
Notice that in Eq. (27) we have disregarded the electrostatic term of the Kelvin force density, as it has already been included in Eq. (26). This is just a consequence of the macroscopic Lorentz force being simply the sum of the individual ion contributions.
In summary, the external force term in the Navier-Stokes equation (4) is given by After some manipulation, and defining the total ionic concentration n(r r r) ≡ n + (r r r) + n − (r r r), Eq. (28) can be written as Finally we present the discretization of the force terms. The forces on the solvent due to microscopic interactions are given by Eq. (17), with the already specified coupling constants and pseudopotentials. For Eq. (26) we use the discretization of the fluxes, Eq. (16), with Therefore, the total force included in the lattice Boltzmann method is given by

D. Poisson equation
A significant challenge in the numerical implementation of electrokinetic models is the efficient and accurate solution of Poisson's equation (11). φ I (r r r) must be determined at each time-step from the charge distribution q(r r r) and the permittivity ε(r r r). Here we use a standard Fourier spectral method, which allows us to find a solution for the potential in O(G log(G)), where G is the number of lattice sites. A significant advantage of this method is the possibility of efficient parallelization. On the other hand, boundary conditions on the electric potential can be hard to implement, although here we consider only periodic systems.
One additional difficulty in our case comes from the space dependent ε(r r r). In that case, the general Poisson equation (11) has to be solved, which can be written as ∇ 2 φ I (r r r, t) = −q(r r r, t)/ε(r r r, t) − ∇ε(r r r, t) · ∇φ I (r r r, t). (32) As ε(r r r, t) is known, only the ∇φ term on the right presents a problem. In order to decouple this term and be able to use the standard discretization of the Poisson equation, we assume a slowly varying potential, such that φ I (t) ∼ φ I (t − ∆t), and thus ∇ 2 φ I (r r r, t) = −q(r r r, t)/ε(r r r, t) − ∇ε(r r r, t) · ∇φ I (r r r, t − ∆t).
(33) The condition of a slowly varying φ I (r r r, t) is essentially a demand for a slowly varying q(r r r, t), a condition already set by the total electric force by assuming low Mach numbers and diffusivities, which imply small advection and diffusion of ions, respectively.

E. Colloidal dynamics
Having specified the model and numerical implementation for the hydrodynamics and ion kinetics, only the dynamics of the colloidal particles rests to be determined. For simplicity, we consider only spherical and rigid particles of radius r p . Particles follow Newton's equation of motion with j = 1, ..., N the particle index, N the total number of particles, m the mass, x x x j the position and G G G G G G G G G j the total force exerted on particle j. For simplicity we take all particles to have the same mass m = ρ p (4/3)πr 3 p , with ρ p = 5ρ 0 , an arbitrary choice which has been seen to have no influence on the results. The total force includes both particle-particle and particle-fluid interactions, as well as an electrostatic term, with V (x x x j , x x x k ) the interaction potential between particles j and k, G G G G G G G G G j fl the force exerted on the particle by the fluids and G G G G G G G G G j lub the lubrication correction. Also present is the force exerted by the electric field on the particle, which is assumed to have a charge Q j homogeneously distributed throughout its volume.
The interaction potential is given by the Hertzian model with the direction between particles' centersn n n jk = (x x x j − x x x k )/||x x x j − x x x k ||, and δ the overlap between the two particles, δ = max{2r p − ||r r r j − r r r k ||, 0}. We take the stiffness to be k H = 1.0(m 0 ∆t −2 ∆x −1/2 ), although we see no influence of its value in any of our considered system due to the small overlaps and colloid volume fractions involved. The lubrication force G G G G G G G G G j lub models an observed repulsive interaction between two particles approaching each other in a fluid medium [21]. It has its origin in strong pressure gradients generated by the flow of fluid out of the gap between the particles as they approach each other. Here we follow the procedure of Ladd et al. [21], including lubrication effects directly as a force of the form with the particle velocities v v v j = dx x x j /dt, and the sum over k runs over all particles such that k = j and ||x x x j − x x x k || − 2r p < 2 3 ∆x. These two forces provide a simplified description of what are in reality highly complex colloid/colloid interactions. A more involved potential than Eq. (36) would be needed to model the interactions at small distances, where van der Waals and other intermolecular interactions are relevant. The derivation of the correct interaction potential at these scales is not trivial. For example, the commonly used DLVO potential is not directly applicable, as the electrostatic interaction between Debye layers is expected to be already resolved at large distances, although not when the separation between the colloids becomes small, ||r j − r k || < ∆x. Moreover, a more general description is needed for colloids at fluid interfaces, a topic of active research [39][40][41][42][43]. Finding a correction for electrostatic interactions, such as Eq. (37), is suggested as a future direction of research.
Finally, only the force resulting from the interaction between fluids and colloids rests to be specified, G G G G G G G G G j fl . This is resolved by projecting the particles onto the lattice, marking as solid the sites that are covered by the particle. The interaction is then a consequence of the boundary conditions applied to fluid advecting towards particle solid nodes. Furthermore, as the particle moves, two other processes need to be specified: the conversion of fluid sites to solid, and the conversion of solid sites to fluid. Fig. 1 shows an illustration of these processes.
Fluid advecting onto particle solid sites is bouncedback half way between the solid and the fluid sites, by setting where c c c d = −c c c d * . The second term corresponds to a correction that takes into account the particle movement [21,22,44]. The force density transferred to the particle from every bounce-back collision is then g g g g g g g g g j bb (r r r i ) = When creating solid sites, as the particle moves and occupies a new lattice site, the momentum of the fluid being covered must be transferred to the particle, in order to ensure global conservation of momentum. Therefore an additional term has to be added to the particle, Moreover, in the last case of a particle uncovering a solid site, new fluid is created from an equilibrium distribution with the velocity of the particle at that site, such that The newly added densityρ σ corresponds to the average of the neighboring sites, where the sum runs over the sites' neighbors. Correspondingly, the final contribution to the total force on the particle is The total force from fluid interactions is then obtained by summing all individual contributions at every timestep, such that G G G j fl = ( j g g g g g g g g g j bb (r r r i ) + g g g j c (r r r i ) + g g g j d (r r r i ))∆x 3 , with the sum running through all relevant lattice sites of each process.
A final correction due to the presence of colloids in binary fluids involves the computation of Shan-Chen forces. As can be seen in Eq. (17), the density of the neighbors of a given site need to be accessed to determine the pseudopotentials, but right next to the particle these might be inside the particle. When solid sites are just considered to have no fluid the result is an artificial increase in the density next to the particle [22]. To solve this problem we use the methodology proposed by Jansen et al [22], and set the interface sites inside the particle to have a virtual fluid concentration equal to the average of the fluid neighboring sites,ρ σ . No advection or collision steps are performed on the virtual fluid sites. Nonetheless, these are considered for Shan-Chen forces, balancing the forces at the fluid side and thus preventing the formation of a high density layer [22].
For the boundary conditions of the ion fluxes with the colloids we follow an equivalent procedure as proposed by Kuron et al. [24]. As in general an electrical double layer is formed around the colloids, where ion concentrations are considerably larger than in the bulk, the displacement of charges due to the creation and removal of solid sites leads to sudden electric field variations, which result in large fluctuations of the particles' velocities. In order to reduce these discretization effects, the flux of ions is modified to take into account the volume fraction covered by the particle at each site. The overall procedure was inspired by Noble-Torczynski hydrodynamic boundary conditions, which generalize the previously described bounce-back procedure for sites that are partially filled by a solid [45]. The solid fraction ψ s (r r r i ) is determined at every lattice site, with ψ s = 0 for lattice site volumes that do not intersect with any particle, and ψ s = 1 for volumes completely inside a particle, that is, solid sites (see Fig 1). For the cases of partial overlap, the solid fraction is exactly given by the intersection of a sphere with a cube of side length ∆x at the position of the site r r r i . As there is no simple analytical solution for this geometrical problem, ψ is approximated by subdividing each lattice site h-times, and counting the number ν of cubic central points that are inside the spherical colloid. It then follows that ψ s = ν/(h + 1) 3 . The ionic fluxes (Eq. (16)) are modified to take into account the solid fraction field ψ s , such that at equilibrium the amount of ions at a site is proportional to the fluid fraction, ψ f = 1−ψ s . The fluxes then take the form The creation and deletion of solid sites has to also take into account the ions present in the fluid. At the creation of a solid site, ions are isotropically displaced to the fluid neighboring sites. Although this might seem artificial, it is always a relatively small amount of ions, as ψ s ∼ 1 (see Fig. 1). Moreover, it is the most straightforward method that is both local and directly conserves ions. Similarly, when removing a solid site, ions in the newly created fluid are simply set to zero.
It is important to notice that as the colloid moves, the number of solid sites it occupies is not constant; in order to keep a constant and homogeneous charge on the colloid, we redistribute Q i at every time-step so that the colloidal charge at each site q c (r r r i ) = (Q i (r r r i )/v p )ψ s (r r r i ), with the particle volume v p = (3/4)πr 3 p . Finally, the permittivity of the particles ε p (r r r) also has to be considered when solving Poisson's equation. For sites that are partially covered by the particles we take a simple interpolation using the solid fraction field between the permittivity of the fluids (already an interpolation depending on fluid concentration, as specified in Eq. (13)), and the permittivity of the particles, such that

A. Fluid-fluid interface
In this section the ion distribution next to a planar fluid interface is investigated. We validate the ion diffusive flux and corresponding hydrodynamic force given by solvation. Two different solvents are considered, which have different Gibbs transfer energies ∆µ ± = 0, as well as different permittivities, ε a = ε b . The results are compared with an analytical solution for ion distributions at interfaces derived by Onuki et al [46]. Bier et al. give a solution for the electrostatic potential of an interface at x = 0, with fluid a for x < 0 and b for x > 0 [47]: Here A = (1 + b ) cosh(κ i s) + (b + ) sinh(κ i s); b = ε a /ε b ; = n a /n b , with the subscript denoting the bulk values of the respective solvent. The screening lengths are given by κ a,b = (2βe 2 n a,b /ε a,b ) 1/2 and at the interface κ i = 2βe 2 n b /ε a . The potential difference between the two bulk phases φ D ≡ φ b − φ a is the Donnan potential. The parameter s shifts the potential as a way to capture interfacial effects [47], such as a smoothly varying ε(x), present in our case.
The same setup was previously considered by Capuani et al [17], where solvent/solvent interactions were modeled using a free-energy based on the concentration field c(r r r), instead of using Shan-Chen pseudopotentials as in this study. In this respect, the most significant difference between the two approaches is the equilibrium state of a demixed mixture. While in our case there is a small but significant fraction of the minority fluid at the bulk phases, the concentration field in the free-energy formalism is truncated such that at the bulk it takes the perfectly demixed values, |c| = 1 [17]. Nevertheless, the theoretical derivation of Onuki does not assume perfectly demixed phases [46]. In our case, agreement with the theoretical predictions was found only for the solvation potential normalized by the differences in concentration ∆ρ a , as in Eq. (24). When the expression in Ref. [17] is used, which assumes perfect species separation, the convergent values of φ PB (x = ±∞) are not recovered, as the Donnan potential in that case does not take the form (∆µ − − ∆µ + )/2e. We thus remark that the solvation potential in Eq. (24) is a more general case that does not assume perfect phase separation.

B. Dielectric droplet deformation
In the following test, a droplet of solvent b (oil) immersed in solvent a (water) is considered in an external electric field E E E e = (0, 0, E z ). In the absence of an electric field, the balance of surface tension and pressure determines the shape of the droplet, which minimizes the surface area: a sphere in 3D, a disk in 2D. When an external field is applied and the permittivities of the solvents are different, the droplet can elongate either in the direction of the applied field or perpendicular to it, depending on the permittivity contrast. The final shape is set by a balance between Laplace's pressure and the Kelvin force density. Notice that, as the system is ion-free, it allows us to independently test the dielectric force applied to the solvents, as in (29) the only non-vanishing term is the Kelvin force F F F E , Eq. (27).
As deformations are expected to be symmetric in the axis of the field, as also to optimize the computation time, only quasi-two-dimensional geometries are considered, by setting l x = 4∆x, while l y = l z = l. The deformation δ is defined as the normalized difference of half the droplet length in theŷ andẑ direction, δ ≡ (b − a)/(b + a), with a inŷ and b inẑ. For small deformations the equilibrium shape is an ellipse, and thus a and b correspond to the semi-minor and semi-major axes, respectively. We compare our results to the existing analytical solution in the limit of small deformations [48], with r d the undeformed droplet radius and γ the surface tension. The surface tension is obtained using the Young-Laplace equation, ∆p = γ/r d , where the pressure difference ∆p is taken to be between the center of the droplet and a point outside, computed from the equation of state of the Shan-Chen model, Eq. (20). The size of the droplet and thus δ is determined at sub-lattice resolution by fitting the local concentration profile c(y, z) next to the droplet interface by a hyperbolic tangent, and defining the interface position as the points where c = 0. As in this configuration all electric forces, except for dielectric, are zero, the degree of the deformation is determined by the relative magnitude of the electric forces to the interfacial tension stresses, characterized by the electrocapillary number Ca E = ε a r d E 2 z /γ. The predicted deformation can thus be rewritten as δ th = Ca E χ 2 /4.
Boundary conditions are taken to be periodic, and thus the simulation can be considered to be of a square lattice of identical droplets. We study the system size dependency by considering two normalized droplet sizes,r d ≡ r d /l ∈ {0.1, 0.25}, for each system size l ∈ {64∆x, 128∆x}, corresponding to volume fractions ν ∈ (0.03, 0.2). Initially a circular oil region is set to relax for 10 5 ∆t, and each simulation is then run from this equilibrated state for 2 × 10 5 ∆t. Both the strength of the applied field E z and the permittivity contrast χ are varied. The droplet is seen to deform in the direction of the applied electric field, as shown in Fig. 3a. The effect of the differences in permittivity is evident when looking at the electric field stream lines (Fig. 3b), which result in a dominant dielectrophoretic force at the poles of the drop in theẑ direction (Fig. 3c). This force is eventually balanced by the increase in surface tension due to the higher curvature of the deformed droplet.
An excellent agreement between the computed deformation and δ th is obtained forẼ z < 0.1 and r d ≥ 30∆x (see Fig. 4). (Here we have defined, for clarity, the dimensionless electric fieldẼ z = E/E 0 , with E 0 = m 0 e −1 ∆x∆t −2 ). In these ranges, the normalized defor- mation is seen to be independent ofẼ z and shows the expected quadratic dependence with χ. This dependence breaks forẼ z > 0.1 and high χ; we do not expect Eq. (47) to be valid for these deformations. Smaller droplets show errors within 15% even for smallẼ z , most probably as a consequence of the finite-size diffuse interface. The results are consistent with those obtained in [17], where the solvent's interactions are modeled through a free-energy model, in such a way that the interface width and the surface tensions are parameters of the system. In our case, both of these quantities are a result of the Shan-Chen model, solely determined by the interaction parameter G, and are measured from the resulting equilibrium state. Increasing the electric field further leads to large elongated droplets with pointed ends, as shown in Fig. 5. Similar shapes have been observed in experiments for drops under shear for low viscosity ratios between the two fluids [49]. When the field is turned off, the drops return to their equilibrium shape, and no breakup is observed even for the highestẼ z considered.

C. Neutral electrolyte droplets
In the previous case no ions were considered, the deformation of the droplet being driven by the difference in permittivity between the two fluids. In the following, ions are added to the fluids, and the deformation of a drop in an electric field is studied. Modifying the Gibbs transfer energies ∆µ ± allows us to produce ratios of bulk ion concentrations between the two fluids of up to 32   100. At this limit we expect the influence of the ions in the outer fluid to be negligible, and therefore to be close to the most common experimental configuration where the drop's fluid has a much higher conductivity than the medium. In order to isolate the effects of the ions from other contributions, we take both fluids to have equal permittivities, χ = 0, as well as equal densities and relaxation times, implying equal viscosities. We investigate the qualitative aspects of the deformation for variations in the average concentrationn = n(r r r)dr r r/l x l y l z , and E z .
As the total electric force applied on the droplet now has three contributions (cf. Eq. (29)), it is not evident that Ca E retains its relevance as it, for example, ignores the effects of increasing ion concentrations, which increase the electrostatic forces that lead to deformation.
We perform simulations of quasi-two-dimensional periodic domains of size {l x , l y , l z } = {4∆x, 128∆x, 384∆x}. The system is initialized with a circular droplet of fluid b of size r d = 32∆x, and a homogeneous ion concentration n = 10 −3 ∆x −3 . We first perform a relaxation phase during 2 × 10 5 ∆t with E z = 0, until diffusive fluxes vanish and the ion densities reach an equilibrium distribution. The Gibbs transfer energies are set to ∆µ ± = −4k B T , which results in a ratio of bulk ionic densities between the two species of n b /n a ∼ 80. For simplicity we set equal permittivities on both fluids, such that ε f r = 80, with ε f r = ε f /ε 0 . Comparing the measured Debye length λ D = ε f k B T /e 2n with experimental values known for water gives ∆x = 1nm, for ionic concentrations in the dilute limit 10 −3 mol/l.
Although we expect finite-size effects to be relevant, we observe that the total electric field E E E at the boundaries is much smaller than E z , and thus the qualitative aspects are not expected to be significantly altered by the domain size.
The electric field polarizes the drop as ions with different charges move in opposite directions. As a conse-quence, the electrostatic force, acting on opposite directions at each pole, deform the droplet, as shown in Fig. 6. For low E z the drop elongates and reaches stable ellipsoidal shapes. Beyond a critical E z the deformed drop becomes unstable, and after acquiring a dumbbell shape it breaks up in a pinch-off manner [50]. The morphology of the breakage varies with E z : slow elongations lead to round drop ends and thick connecting necks, with no satellite drops formed at breakup. Higher velocities, for higher applied fields, lead to greatly deformed ends perpendicular to the flow direction, a consequence of strong fluid drag and dielectric forces pushing in opposite directions. Capillary forces eventually break the slender neck, and two or more satellite drops are formed. Overall the qualitative aspects of the deformation and breakup coincide with previous experimental observations [51], and numerical solutions of leaky-dielectric models in various limits [50,[52][53][54]. Here we have shown the capability of the model to study drop breakup in electric fields; further studies could quantify the effects of ion distributions at different concentration limits, providing information unreachable by previous numerical methods which assume all charges to be at the interface [55]. Future quantitative studies should also take into account the known dependency of the degree of deformation and breakup with the grid resolution [56].

D. Charged droplets
We now look at the deformation of droplets with a nonzero electric charge. A quasi-two-dimensional periodic box of size {l x , l y , l z } = {4∆x, 384∆x, 384∆x} with a droplet of radius r d = 32∆x is considered. The same relaxation procedure as in the previous system leads to the droplet having a negative net charge, when setting different Gibbs transfer energies for the two ion species ∆µ ± = ±4k B T . The concentration of ions is fixed at n = 10 −3 ∆x −3 , which sets the average charge in the droplet atq d ∼ −3 × 10 −3 e. The permittivities are set as in the previous systems, ε r = 80. We vary the electric field E z .
As expected, drops accelerate in the direction of the electrostatic force qE E E. Eventually they either reach a terminal velocity with an equilibrium deformed shape, or go through a breakup process. Initially the spherical shape widens perpendicular to the direction of flow, acquiring an oblate spheroidal shape [57]. Afterwards, the drop bends, acquiring a bell-like shape with lobules at the ends connected by an increasingly long and thin neck, as shown in Fig. 7. The ultimate breaking process is analogous to pinch-off, when surface tension forces break the thin neck. Depending on the strength of the field, at breakup the neck either withdraws to the end lobules, in which case only two drops are produced (Fig. 7a), or breaks up further into smaller drops, a situation seen for the highest E z considered (Fig. 7b).
Deformation and breakup are a result of hydrodynamic   forces due to the relative velocity of the drop with the ambient fluid [58]. Previous studies of macroscopic drops and/or bubbles often refer to this process as secondary atomization, secondary breakup or droplet disintegration [58,59]. We observe that the breakup process for a charged droplet in an electric field is similar to the bag-breakup process known from previous experiments and numerical simulations of drops or bubbles flowing in either a fluid or a gaseous medium. After initial deformation, the center of the drop gets pushed downstream forming, in three-dimensions, a bag attached to a toroid. In our two-dimensional case, the analogous situation is of the observed neck connecting two separate lobules. Eventually the bag breaks, and is observed to fragment into a number of smaller drops [58,59]. We observe fragmentation only for the highest electric fields (see Fig. 7), most probably as a consequence of the resolution limiting the size of smaller drops. In this last case, the bag is observed to create a plume at its center, most probably a precursor of a breakup process referred to as multimode breakup in experiments which, consistently, is seen to occur when increasing the forcing beyond bag breakup [58]. The similarities of charged and uncharged phenomenologies can be further supported by considering the observed invariance of the breakup process with two dimensionless numbers: the Weber number, defined as the ratio of inertia to surface tension, We = ρ a u 2 d r d /γ, with u d the terminal velocity of the drop; and the Ohnesorge number, taken as the ratio of the drop's viscous to surface tension forces, Oh = η a / √ ρ a r d γ. In our systems, Oh ∈ (0.1, 1.0), and We ∈ (10, 30). At these values, various experimental studies have found the breakup process to be at either the bag or multimode type [58,59].
Until now only quasi-two-dimensional systems have been considered, as the deforming droplets were always expected to be axisymmetric in the direction of the applied field. We wish to remark that our numerical implementation is performant enough to simulate threedimensional systems of the cases previously presented. We realized individual simulations of selected cases and always observed qualitatively equivalent behaviours. As an example, we show in Fig. 8 the break-up process of a charged droplet of size r d = 64∆x, with parameters as in Fig. 7b, in a periodic cube of side length l = 384∆x. In this geometry one can clearly see how the elongated-disk drop breaks initially at its center, forming a toroid. Negative charges, initially forming a spherical double layer, are quickly accumulated at the front of the droplet. The fluid flux also significantly advects the charges through the outer and inner (after breakup) sides of the drop. Overall, the breakup process is consistent with the quasitwo-dimensional drops.
In summary, we conclude that the effect of the double layer on the dynamics of breakup is not highly significant in this parameter range, as the drop breakup process is to a large extent analogous to a non-charged drop deformed under the effect of hydrodynamic forces. Nevertheless, we do not expect this to be the case for other salt concentration limits. In general we have shown that our numerical methodology is capable of capturing the qualitative aspects of this process in the right order of the relevant dimensionless numbers. Future studies could systematically vary the concentration of ions to quantify the influence of electric forces on drop deformation and breakup.

E. Colloid electrophoresis
In the following sections electrokinetic colloidal suspensions are studied. In order to validate the coupling of the colloidal particles to solute and solvents, we first consider the basic electrophoretic setup in which a single, charged particle moves in an electrolyte solution under the effect of an electric field. We compare our results with experimental measurements [60] and lattice Boltzmann simulations [61]. In the latter study [61], the link-flux method was used to solve the electrokinetic equations in a similar manner as in the present study although, importantly, only fixed particles were considered, while in our case the particles are free to move.
The velocity at which a colloidal particle with charge Q moves in a weak electric field is linear with E, with the proportionality constant refered to as the electrophoretic mobility µ = v/E. A colloidal particle is usually surrounded by counter-ions (ions with a total charge −Q), which create an electric field opposing the movement of the particle. The charges on the colloid and the cloud of counter-ions is referred to as double layer. The rest of the electrolyte ions in the solution, referred to as salt or co-ions, increase the drag on the particle and decrease the electric field, thus reducing the mobility µ of the particle relative to the case with no salt [62]. The theoretical derivation of µ as a function of the charge of the particle Q, the salt concentration and, in the case of solutions, the solid fraction η is not trivial, and has been done only for the limiting cases of small or large double layer thickness [62]. In the following we validate our simulations by comparing µ with previous numerical work and experiments, usually realized at the mili-or microscale.
A single colloid of radius r p = 4.0∆x is situated in a periodic box of size l. Due to the use of periodic boundary conditions the system can be considered as a cubic colloidal crystal of solid fraction η = (4/3)πr 3 p /l 3 . We modify the solid fraction by varying l ∈ (24, 128)∆x, that is ν ∈ (10 −4 , 10 −2 ), and set the material parame- ters as in previous experiments and simulations [60,61]: colloid charge Q = 30e, and a Bjerrum length λ B = e 2 /(4πεk B T ) = 1.3∆x. The permittivity of the fluids and particles is set such that ε f r = 80, and ε p r = 10. A good agreement is observed with the previously measured dependency of the electrophoretic mobility with the solid fraction, µ(ν), as shown in Fig. 9. As in previous works, we have expressed our results using the dimensionless mobilityμ = 6πηl B µ/e. Consistent deviations with experimental values increase with the applied field, although they remain within a 5% error for E z ∈ (10 −4 , 10 −2 ). We expect the deviations at high forcing to be a consequence of the fluid terminal velocity being close to c s , namely Ma ∼ 1, where many of the model assumptions are no longer valid. On the other hand, for very low fields,Ẽ z = 10 −4 , the very low velocity of the colloid is significantly affected by discretization effects, not only due to flux of ions but also due to the hydrodynamic boundary conditions. These generate large fluctuations in µ (as can be seen by the error bars in Fig. 9), although the average remains consistent with previous studies for ν > 10 −2 .

F. A particle at a fluid interface
As final exemplary cases for the possibilities of the numerical method, particles at interfaces of electrolyte solutions are studied. We begin by considering the behavior of a single particle at a planar fluid interface, and afterwards show the possibility of simulating colloidal suspensions on the surface of droplets. These cases represent common scenarios of scientific and technological relevance. Nanoparticles at fluid interfaces are common in various in-development technological advancements, such as advanced coating processes for molecular electronics and optical devices [63,64], stabilization of emulsions [65] or better control and transport capabilities in micro-and nanofluidic devices [66].
Colloidal particles at fluid interfaces are highly stable, as the particle-oil and particle-water surface tensions are µ Exp. fluid [60] Exp. crystal [60] Sims [60] Sims [60] Sims [61] FIG. 9. Dimensionless electrophoretic mobilityμ as a function of the solid fraction ν, for a colloidal particle of size rp = 4.0∆x, and charge Q = 30e. Also shown is experimental data for suspensions of latex particles of diameter 64nm in a fluid (+) or crystalline (×) state, as detailed in [60]. Furthermore, data is shown for simulations using a lattice-Boltzmann/Molecular Dynamics combination as described in [60], for particles of total charge Q = 20e and radius 4∆x ( ), and Q = 30e and radius 6∆x ( ). Finally, data for simulations using an equivalent algorithm as this study, although with fixed particles, is also included ( ) [61].
in most cases much smaller than the oil-water one. Although the energy needed for detachment of the interface scales quadratically with the radius of the particle [67], at the nanoscale the energies needed for detachment are still 10-100 times larger than k B T [68,69]. Therefore we ignore, as a first approximation, temperature fluctuations.
A single particle of radius r p = 8∆x is placed at a fluid interface in a periodic box of size l = 64∆x. The radius of the particle was chosen to minimize the effects of the finite width of the interface. As the system is periodic, the setup is equivalent to an infinite crystalline colloidal suspension with solid fraction ν = 8 × 10 −3 .
The interface without the colloid is previously relaxed for 10 5 ∆t, such that the concentration of ions, initially homogeneously set to n(r r r) =n, reaches an equilibrium distribution given by the different Gibbs transfer energies ∆µ ± . These are fixed at ∆µ ± = 4.0k B T . Permittivities are set as in the single colloid electrophoresis case, ε f r = 80, and ε p r = 10. We varyn and the charge of the colloid Q, and measure the displacement of the colloid from its equilibrium position with no solutes, z * 0 , z * = |z c − z * 0 |. The displacement of the particle is seen to increase linearly withn, as shown in Fig. 10. At a critical salt concentration the ion's osmotic pressure is strong enough to detach the particle from the interface. Close to this limit, the linear behavior is lost, and the displacement saturates, probably an effect of the steep interface deformations interacting with the discretized particle shape.
The equilibrium position of a particle at a fluid interface is usually quantified by the contact angle θ c . This is the angle formed by the tangent of the interface and the tangent to the particle at the triple contact point. We observe θ c to be independent of the displacement of the particle. This is to be expected, as the wetting properties of the particle have not been modified, and even though the particle has displaced, the interface deforms in such a way that θ c is kept constant, θ c = 90 o in our case. In addition, we wish to remark that having access to θ c in experimental situations at the nanoscale can present a significant challenge. More commonly the apparent contact angle θ a is measured, formed by the tangent of the interface at the triple contact point with the horizontal. In that case we naturally observe a variation with salt concentration, as sin(θ a ) = z * /r p .
Surprisingly we see no strong effect of the charge of the particle in z * 0 , for all considered charges and salt concentrations. As shown in Fig. 10, charged particles generate a double layer, with most of the charge concentrated on the conducting side. Nevertheless, the resulting electrostatic force is observed to be at most an order of magnitude lower than both the ions osmotic pressure and the solvation forces. As the total concentration of ions around the particle is mostly independent of Q (the increase in negative ions next to the colloid compensates the loss of positive ions), the total electric force F F F e only slightly varies with Q, and therefore the equilibrium position is only marginally affected.
In summary, we have briefly shown the possibility of studying nanoparticle adsorbtion at fluid interfaces with resolved charge distributions, using the presented numerical methods. Future studies could explore the stability of the particle as a function of electrolyte ion concentration, particle charge, and the possibility of higher control using an external electric field. Research in this direction could help to improve our understanding of electrodipping forces at the nanoscale [70,71], a necessary step to obtain a general interaction potential for particles at fluid interfaces.

G. Colloid coated droplets
Finally we demonstrate the possibility of simulating the dynamics of colloidal suspensions adsorbed at curved fluid interfaces, such as the surface of a droplet. Colloidcoated droplets are present in many physical and chemical processes. In Pickering emulsions, the coating of droplets by colloids can stabilize emulsions without the need of surfactants. Emulsions can then be used as templates for fabrication of capsules composed of colloids: colloidosomes [65,72,73]. Colloidosomes can be functionalized as containers of different substances, such as drugs, general reagents, or even cells, with several potential applications in the food, biomedical and petroleum industries [73].
In this case we consider a droplet of radius r d = 32∆x covered by an homogeneous suspension of particles of size r p = 5∆x. The droplet parameters are identical as for the deforming charged droplets studied Section III B, as we wish to investigate the influence of the particles in the overall deformation of the droplet. For simplicity we consider all particles to be uncharged, Q j = 0. The total number of particles is N = 80, which gives a covering fraction of ν s = (1/4)N r 2 p /r 2 d ≈ 0.48. Permittivities are set as in the previous cases, ε f r = 80, and ε p r = 10. The effect of the covering particles on the deformation and breakup morphology of the droplet is dramatic. While without colloids the droplet underwent a process similar to bag breakup (see Fig. 8), in this case the breakup process is more similar to pinch-off: the drop deforms and squeezes through the particle coating, forming t = 6000Δt t = 2000Δt t = 10000Δt t = 14000Δt FIG. 11. Snapshots of a particle-coated charged droplet in an electric field. Parameters are identical as in Fig. 8, with rp = 5.0∆x. The drop interface c(r r r) = 0.0 (green) and colloids (white).
a long neck until it eventually breaks, as shown in Fig. 11. The results is the creation of a new, particle-free droplet, and the increase of ν s in the original droplet, such that it is now stable in the electric field. Including corrections to take into account the overlapping ion clouds at small distances would certainly reduce the equilibrium packing fractions, although we expect to observe the same qualitative dynamics. Due to the increased inertia of the coated droplets, their speed is considerably slower than the particle-free ones, so that drag forces do not significantly deform its shape. Overall we see phenomena worthy of future studies which might have practical significance, as the covering fraction of coated droplets could be increased by applying an external electric field and thus removing the excess fluid from the drop.

IV. CONCLUSIONS
A mesoscopic model for the simulation of electrokinetic phenomena of colloidal suspensions in fluid mixtures at the nanoscale was presented. The model follows a more microscopic description than the link-flux model originally proposed by Capuani et al., treating the kinetics of the ions as a response to individual forces instead of using a free-energy functional, and deriving the coupling between the different species using local pseudopotentials. Overall we have shown that binary mixtures of electrolytes can be successfully modeled using the Shan-Chen multicomponent pseudopotentials for both fluidfluid and fluid-ion interactions, providing a new, simple methodology for the simulation of such systems. Colloids are included via a Ladd coupling with the fluids, and a solid-fraction scheme of discretization that significantly reduces discretization errors. Results were shown to be consistent with the previous link-flux method, and agree with known theoretical solutions for ionic concentrations at fluid interfaces, droplet deformation and electrophoresis of a single colloid.
Several systems were shown to explore the posibilities of the model. Droplet deformation and breakup for neutral and charged droplets presents the same qualitative aspects as previous numerical and experimental studies at the microscale. The ability to integrate particles and fluid mixtures opens the possiblity to study at a new level of resolution the statics and dynamics of colloids at interfaces, both planar and curved. Our first investigations show that colloids at a planar electrolyte interface vary their equilibrium position when ions are added to one of the fluid, until the point were adsorbtion becomes unstable and the particle detaches. Furthermore, it is now possible to study (charged) colloidal suspensions adsorbed at the surface of a droplet; we have here shown a first exemplary case of a colloid-coated droplet breakup due to an external electric field. All these examples illustrate that the presented method can be used for studies of a variety of electrokinetic phenomena in parameter ranges so far unreachable.
This description of ion transport in electrolytes could be easily extended to include other effects. Additional forces acting on the ions can be readily added in Eq. (10), with the ion-solvent coupling force following directly from the friction coupling. More involved extensions beyond the Poisson-Boltzmann limit of no ion-ion interactions or solvent polarization, could include the use of modified Poisson-Nernst-Planck equations (MPNP) that take steric effects into account [74]. These could be included as additional terms in the ion flux, Eq. (9), and are possibly crucial in confined systems at high ion concentrations or strong external fields [74].
In summary, our work has shown that mesoscopic simulations can be used to study electrokinetic phenomena at the nanoscale, where obtaining experimental data becomes a challenge. The same algorithm could be applied to a variety of nanofluidic systems, such as droplet transport in inhomogeneous geometries, coalescence and generation of droplets, and nanomixers.