Ignacio Pagonabarraga
email: ipagonabarraga@ub.edu

Daan Benjamin Rotenberg

D Frenkel

Recent advances in the modelling and simulation of electrokinetic effects: bridging the gap between atomistic and macroscopic descriptions

Recent advances in the modelling and simulation of electrokinetic effects: bridging the gap between atomistic and macroscopic descriptions

Introduction

Charged materials are ubiquitous in nature. This abundance stems among other from the outstanding ability of water to dissolve charges. Its high electric permittivity favors ionic dissociation, and most water-soluble macromolecules and colloids are charged. Charge thus plays a key role in their interactions and significantly affects the collective properties of a large variety of materials. Electrostatics is also involved in the properties of biological systems at cellular and molecular scales 1 , where electrostatic correlations play an important role in their equilibrium properties 2 . From the dynamical point of view, the interaction between charges is only one aspect of the problem, and the coupling between solvent and charge flow, referred to as electrokinetic phenomena, must be properly taken into account.

Electrokinetic effects have long been studied in colloidal suspensions. The electrophoretic mobility of colloids results not only from the direct effect of the applied electric field on the colloid, but also from the force acting on the suspending fluid which is charged in the vicinity of the surface. This same force is also exploited (electroosmosis) to generate a fluid flow near a charged, solid surface, e.g. in microfluidic devices 3 . Recent experimental advances have made it possible to address the behaviour of charged systems in new regimes, such as the electrokinetic response of dilute colloidal suspensions with low salt concentration. In these conditions, the screening of the colloidal charge by the salt micro-ions is less efficient and as a result, the dynamics is more sensitive to the colloid volume fraction than for a suspension of uncharged or strongly screened colloids 4 . This scenario, most easily achieved for small colloids, is relevant for protein suspensions.

Electrokinetic effects also play a role in biological systems. For example, vesicles respond to applied electric fields 5 , giving rise to a variety of phenomena, including membrane flow, electroporation and fusion [START_REF] Haluska | Proc. Nat. Acad. Sci[END_REF] . Nonlinear induced electrokinetic flows influence the properties of membranes outside of equilibrium [START_REF] Lacoste | [END_REF] and the potential drop across membranes can destabilize them 8 . Electrokinetic effects have also been exploited in novel imaging techniques, such as ion-conductance microscopy 9 , which can operate under ambient conditions to obtain nanoscale-resolution images of living cells and can be used to probe their local mechanical properties 10 . Modern manipulation techniques make it possible to follow the motion of individual macromolecules in heterogeneous environments. Such studies have made it possible to analyse the electrophoretic response of microtubules and their orientation transverse to the electric field [START_REF] Van Den Heuvel | Proc. Nat. Acad. Sci[END_REF] , arising from a competition between bending and electrokinetic forces [START_REF] Van Den Heuvel | [END_REF] . Single molecule manipulation is of particular relevance in the case of DNA electrophoresis through nanopores 13 . As the ionic current is sensitive to the charge of the polymer segment in the pore, experimental studies of ionic currents during DNA translocation can in principle be used to discriminate between different types of nucleic acids 14 . There is also growing interest in the behaviour of dispersions of liquids of widely differing permittivity. Experiments on water droplets in oil have shown that charge accumulated at the interface between the two phases can stabilize water-in-oil dispersions [START_REF] Leunissen | Proc. Nat. Acad. Sci[END_REF] . Finally, understanding of the coupling of charge transport and flow in porous materials is of relevance for the study of dispersion of radio-active ions in the soil and for the understanding of seismo-electric phenomena (see e.g. [START_REF] Pride | [END_REF] ).

The theoretical description of the equilibrium properties of charged colloids has a long history, going back to the seminal work of Gouy, Debye and Hückel and, more specifically for colloids, the work of Derjaguin, Landau, Verwey and Overbeek. Here, we will not review the literature on the equilibrium theory (for a recent review that includes references to the historical literature, see ref. 17 ). Suffice it to say that determining the equilibrium properties suspensions of charged colloids, emulsions and porous materials is already difficult, because of the heterogenous nature of the systems (solid-liquid or liquid-liquid interfaces, internal structure of dissolved macromolecules...), and the range of length scales involved: Molecular size (microions and solvent molecules at the surface), diffuse layer over which the fluid is charged, radius of the colloidal particle or pore size, distance between colloids. Moreover, the dynamical phenomena associated with each of these length scales can occur over a wide range of time scales. This multiscale nature of the problem makes it difficult to treat all species on the same footing, e.g. using all-atom Molecular Dynamics simulations.

When modelling electro-kinetic effects, the challenges posed by the long-range nature of electrostatic interactions are compounded by the fact that hydrodynamic interactions are also long-ranged. Hence, both in the structure and the dynamics there is a coupling between very different length scales. Moreover, the dynamics affects the structure: hydrodynamic interactions among the macroions change the distribution of the micro-ions, which in turn modify the electrostatic interactions. The material heterogeneity at intermediate scales is intricate in most of the examples discussed. Therefore, approaches to model these systems must be flexible in dealing with bound-ary conditions and cannot be limited to simple, regular geometries. This constraint is especially severe when combined with long range electrostatic interactions. Several approaches have been proposed in recent years to arrive at an integrated description of electro-kinetic phenomena. In this paper, we first recall the basic concepts of electrokinetics (Section 2). We then analyze the (often complementary) strategies that aim at capturing the relevant dynamic processes at different length and time scales (Section 3). Finally, we discuss several illustrative examples of the use of these methods (Section 4).

Basic concepts

Electrokinetic phenomena refer to coupled solvent and charge flow in the vicinity of a charged surface. Electro-osmosis, for example, is the solvent flow arising in the presence of an electric field near a stationary charged solid surface. Conversely, a fluid flow near such a surface transports charges and hence gives rise to an electrostatic potential difference, referred to as the streaming potential. The corresponding electrokinetic effects in the case of charged solid particles in a stationary fluid are electrophoresis and generation of a sedimentation potential, respectively.

Any theoretical description of the phenomena listed above necessarily involves three ingredients: 1) a description of the surface charge distribution (e.g. on the surface of a colloidal particle or on the pore walls in a porous medium), 2) a description of the distribution and transport of the microscopic ions (counterions and dissolved salt) and 3) a description of the dynamics of the liquid phase (or phases, in the case of liquid/liquid interfaces).

It should be stressed that in heterogeneous media such as we discuss here, electroneutrality is only satisfied globally -the fluid carries a net charge in the vicinity of charged surfaces. A simplification of the problem results if we can use the "macroscopic" Maxwell equations to describe the relation between the electric field strength and the charge distribution in the liquid phase with relative dielectric constant, ε r . In that case, the electrostatic potential ψ satisfies the Poisson equation ∇ 2 ψ = -ρ el /(ε 0 ε r ), where ρ el = ∑ N k=1 ρ k z k e is the local charge density due to a N-component electrolyte mixture and e stands for the unit charge, while ε 0 corresponds to the dielectric constant of vacuum . Here ρ k denotes the number density of the ionic species with valency z k . In equilibrium at temperature T , the ionic densities follow, to a first approximation, a Boltzmann distribution ρ k = ρ 0 k e -β z k eψ , where ρ 0 k denotes the number density of species k in a (possibly hypothetical) osmotic reservoir in contact with the system of interest. β = 1/k B T , where k B denotes the Boltzmann constant. For the simple geometry where the solid particle or the fluid interface can be regarded as a planar surface, for small potential drops across the interface, this Poisson-Boltzmann theory can be linearized. In this so-called Debye-Hückel limit, the ionic densities vary as e -κx , where x denotes the distance from the surface and κ

-1 ≡ 4πl B ∑ k ρ k z 2 k -1/2 is
referred to as the the Debye screening length. The Bjerrum length l B = e 2 /4πε 0 ε r k B T is the distance at which the electrostatic interaction between two unit charges equals the thermal energy. In water at room temperature, l B = 7 Å and for a 0.1 mol/L solution of a symmetric 1-1 electrolyte κ -1 ∼ 1 nm. Note that the Debye-Hückel theory ignores all direct correlations between ionic species -in this respect, it is a mean-field theory.

Within this mean-field picture, the time evolution of ionic densities is usually described by the Nernst-Planck equation that accounts for the combined effect of advec-tion, diffusion and migration:

∂ ρ k ∂t + ∇ • (ρ k u) = ∇ • D k ∇ρ k + D k z k k B T ρ k ∇ψ . (1)
The key transport properties in the Nernst-Planck equation are the diffusivities D k of the various ionic species. The local fluid velocity u is assumed to satisfy the Navier-Stokes equation

∂ ∂t u + u • ∇u = ν∇ 2 u - ∇P ρ + f V ρ (2)
where ν = η/ρ is the kinematic viscosity of the embedding fluid of density ρ whilst the local force density f V includes electrostatic forces -ρ el ∇ψ. For example, an applied electric field can only induce flow in the diffuse layer close to the surface, where the fluid is charged. In contrast, in the bulk of the liquid, fluid flow is affected by convective (u • ∇u) and viscous (ν∇ 2 u) momentum transfer, plus any body forces (such as gravity) that act on all fluid elements, irrespective of their charge. The set of equations listed above (Poisson, Nernst-Planck, Navier-Stokes) contains the required ingredients to describe the coupled ionic and solvent fluxes. The numerical solution of these equations is however not an easy task as they are coupled, non-linear partial differential equations. Moreover, the boundary conditions are generally complex (irregular porous medium, moving colloidal particles...). Nevertheless, the details of the flow and of the ionic charge distribution are often not our primary interest: When modelling a colloidal suspension, the aim is to include the effective interactions between macroions, mediated by the solvent and the microions in the diffuse layer. Such a coarse-grained picture may be accurate enough if there is a clear length and time scale separation between the dynamics of the macroions on the one hand and that of the solvent and microions on the other. If the colloidal size, as measured by its radius R, is much larger than the ionic radius, a, it is reasonable to use a coarse-grained description in which the ionic atmosphere is modeled as a continuum.

Even in that case, if ions interact strongly, electrostatic correlations cannot always be neglected. The transition between two regimes (surface-ion or ion-ion dominated) can be understood in terms of the ratio between two characteristic length scales, namely the Bjerrum length l B described above and the Gouy-Chapman length l GC = k B T ε/σ s e (with σ s the surface charge density of the colloid) at which the electrostatic energy of a unit charge in the field of the charged surface is equal to the thermal energy. The so-called strong-coupling limit is then such that z 3 l B /l GC 1. In such cases, a more detailed model that resolves the ionic structure is required to account for ionic correlations.

In the absence of driving fields, on time scales in which colloids relax structurally it is safe to neglect the detailed coupling of the colloids to the solvent dynamics. Nonetheless, one cannot exclude a priori situations where hydrodynamics can have a significant effect on the long-time dynamics of the colloids. What level of description is appropriate to account for this coupling depends on the particular situation. As always, the aim is to use is to use a model that is as simple as possible...but not simpler. For example, the electrostatic interactions between macroions are screened by the diffuse layers and decay as e -κr /r, whilst the hydrodynamic interactions may be accounted for in dilute suspensions via the Oseen tensor decaying as 1/r. These rather extreme simplifications are not always appropriate; they rely, for electrostatics, on a mean-field, frozen picture of the ionic diffuse layer, and for hydrodynamics, on the flow induced by distant spheres in an unbounded medium. Such descriptions assume in particular that the relevant timescales are large compared to the electrostatic relaxation time of the ionic cloud, the Debye time τ D ∼ 1/D ± κ 2 and the viscous relaxation time τ ν ∼ R 2 /ν.

Even under these conditions, if one is interested in the effect that the structure of the ionic distribution will have on driven system, the interplay between electrolyte advection and diffusion, as quantified by the Péclet number, Pe, must be taken into account. For example, at high Pe one cannot assume the equilibrium electrolyte distribution when analyzing the dynamics of such a driven system. An appropriate treatment of the fluid and electrolyte density dynamics must be carried out on the same footing.

Simulation strategies

In principle, all-atom molecular dynamics (MD) simulations of charged heterogeneous fluid media provide a proper means to model both the electrostatic and the electrokinetics of electrolytes. Such studies have allowed, for example, a careful and systematic analysis of the relation between wall-induced ionic ordering and the effective boundary conditions of electrolytes at solid surfaces 18 . Such studies made it possible to determine the limit of validity of the simple mean field approach 19 . However the computational cost of atomistic simulations prevents their use to investigate the large length-scales and long time-scales relevant for colloidal systems or for flow in porous media. Thus a variety of coarse-grained approaches have been developed to meet this modelling challenge. As shown in 1 for the particular case of a colloidal suspension, these different coarse-graining strategies differ in their description of both the charged species and of the embedding solvent, and can be combined to yield a variety of approaches that are complementary to each other, both conceptually and computationally. In general, particle-based approaches treat electrostatic interactions assuming periodic boundary conditions, although in some cases more realistic treatment of the system geometry has been carried out.

Figure 1: The dynamics of charged fluids involves different length scales. Accordingly, different levels of coarse-graining are possible, illustrated here in the case of a suspension of charged colloids. Both the microions and the solvent can be described at three different levels: (1) explicitly as ions and molecules, or (2) explicitly as continuous fields such as ionic densities ρ ± (r) or solvent velocity u(r), or (3) implicitly lumped into effective interactions between colloids. Simple, approximate expressions for the latter include a screened Coulomb interaction for the Potential of Mean Force (PMF) and the Oseen tensor for Hydrodynamic Interactions (HI).

In the following, we discuss some of the coarse-grained approaches before illus-trating their use in specific applications. We stress that our aim is to illustrate rather than to be comprehensive. In particular, our overview covers methods that differ in their description of the solvent that can either be modelled as collection of particles or as a continuum or that, at an even coarser level of description, can only manifest itself via effective hydrodynamic interactions. Similarly, the micro-ions can either be modelled explicitly as particles or as ionic densities, or their presence can be accounted for implicitly through their effect on the rest of the system. The different levels of descriptions discussed in the present article are summarized in 1.

Explicit Solvent

Implicit

Implicit solvent, explicit ions

As most of the space is occupied by the solvent, and ions normally constitute rather dilute suspensions, a number of modelling techniques have focused on retaining information about the molecular details of the charged species of interest, whilst simplifying the description of the solvent and treating it as an effective medium. The solvent motion, even at this coarse-grained level, may affect significantly the dynamics of a charged suspension due to the hydrodynamic flows it carries. However, accounting for hydrodynamic interactions among charged species, especially when dealing with colloidal suspensions, is not obvious, and hence there exist different ways to account for the collective effect of the solvent dynamics on charged suspensions.

Absence of hydrodynamic interactions

Electrolyte properties have often been studied within the framework of the Primitive Model, according to which the uncharged, polarizable solvent contributes an effective electric permittivity, ε r , while ions are described as charged hard spheres. One can then analyze the distribution of charged ions in this dielectric background and the corresponding resulting thermodynamic properties 20 . Within the Primitive Model, the major effect of the solvent on the dynamic properties of the ions disregards the details of the flows they generate but lumps all hydrodynamic effects in an effective friction coefficient, related to the solvent viscosity η by Stokes' law F = -6πηav, with a the ionic radius, and a random force resulting from the thermal collisions between solute and solvent which satisfies fluctuation-dissipation. When the friction is large, the solute velocity relaxes instantaneously and the Langevin dynamics to describe ionic motion simplifies to simple Brownian dynamics (BD).

Inclusion of hydrodynamic interactions

Although the approach described above may be interesting to analyze the structural changes induced by external forces on an electrolyte solution, it fails to provide a correct hydrodynamic description, as the total momentum of the system (ions plus solvent) is not conserved. This shortcoming can be overcome by modifying the effective friction that the moving ions feel as a result of the flow they induce when they are accelerated. More precisely, the flow generated by a moving point ion is not resolved explicitly but only the resulting force acting on a second, distant ion is included in the form of an interaction Oseen tensor. In the spirit of the Primitive Model, the finite size of the ions may be included in the hydrodynamic interaction (HI) tensor, resulting in the Rotne-Prager form 21 . This significantly improves the description of the transport properties of bulk electrolytes, including multicomponent solutions of mono and divalent ions, and micellar solutions 22 . Brownian dynamics of charged suspensions that approximate hydrodynamics by assuming that the Oseen tensor is pairwise additive, have been used extensively to analyze the long-time diffusion of charged colloidal suspensions and to compare both with analytical predictions and with experimental results 23 .

This tensorial approach to HI has also been extended by Netz and coworkers to capture the effect of solid boundaries on the solvent flow 24 . This description, where charged ions are described individually while solvent atoms are regarded as a continuous medium is not fully consistent when the sizes of these two types of particles are of the same order of magnitude, and the associated coarse graining, whilst conceptually understood in equilibrium conditions, does not hold in principle out of equilibrium.

Due to the long-range nature of the HI, approaches based on the Oseen propagator scale very badly with system size, and are in principle valid only for unbounded systems. The approach of Netz and coworkers implements Blake's propagator, which accounts for a single planar wall, also suffers from poor scaling. Building upon earlier works for unbounded media 25,26 , recent studies have proposed to deal with confined geometries using ideas similar to Ewald summation and particle-particle particle mesh [START_REF] Hockney | Computer Simulation Using Particles[END_REF] methods in electrostatics to speed up the slow convergence induced by the long range nature of the HI [START_REF] Hernandez-Ortiz | [END_REF] . However, to our knowledge these methods have not been applied in combination with electrostatic interactions to address electrokinetic phenomena.

An analogous approach has also been developed to speed-up Stokesian dynamics 29 , a method that includes the dynamic effect of the solvent in the corresponding friction matrix which determines colloidal diffusion and response to driving fields. This approach has been used to analyze the relevance of hydrodynamic interactions in the long-time diffusion of charged suspensions 30 . The quantitative comparison with experimental results has clarified the relevance of many-body hydrodynamic interactions. In particular, it has shown the need to go beyond the usual assumption of pairwise additivity to understand the collective diffusion of charged suspensions without having to invoke hydrodynamic screening.

Implicit solvent, implicit ions

In order to reach longer scales, the standard procedure is to average out the degrees of freedom corresponding to the microions, whose dynamics is assumed to be much faster than that of the macroions, and consider only the Potential of Mean Force between the latter. In equilibrium, the prime example of such an approach is the celebrated DLVO theory [START_REF] Verwey | Theory of the Stability of Lyotropic Colloids[END_REF] . The DLVO approach has been extended by a number of researchers (for recent reviews see e.g. 17,[START_REF] Dahirel | [END_REF] ). In particular, it has been possible to go beyond the frozen double-layer picture using density-functional-like approaches 33 or numerical computation from Brownian dynamics simulations with explicit microions 34 . At this stage, these approaches have not yet been integrated with a hydrodynamic description of electrokinetic phenomena.

One possibility to arrive at such a description is to solve the linearized electrohydrodynamic equations to obtain the propagator of the electrolyte solution. In the Debye-Hückel limit Ajdari and Long 35 derived the corresponding Green function of the system as a response to applied external fields. This strategy has the advantage that one can use it to study the Brownian dynamics of a macroion suspension without including the dynamics of the salt and counterions explicitly. Such an approach has been followed by Löwen and coworkers 36 to study the long-time dynamics of charged colloidal suspensions. In particular, when applied to mixtures of colloids with opposite charges they showed that the phenomenon of "laning"' that is observed in the absence of hydrodynamic interactions, is destroyed by HI and more complex patterns emerge. One drawback of the above approach is that it cannot deal with the deformation of the double layer, which is assumed to be always in thermal equilibrium. Moreover, it assumes that the fluid is unbounded and cannot account properly for non-electrostatic external driving forces either.

It is in principle possible to solve the dynamics of the electrolyte in terms of their local number densities, ρ ± using lattice kinetic approaches, exploiting the analogy of the Nernst-Planck equation and the Fokker-Planck equation which can be studied in terms of a phase space distribution f k (r, v,t) with an appropriate kinetic model which conserves particle number. Lattice methods, that are quite efficient numerically 37 can be exploited to achieve large length and time scales in the evolution of these systems at the expense of neglecting dynamic correlations between ions.

Explicit solvent: Particle based approaches

Dealing with hydrodynamic interactions through an implicit solvent is computationally expensive because of the slow decay of fluid flow due to momentum conservation. During the last decade alternative approaches have been developed to treat solvent flow through local rules. One strategy is based on introducing appropriate dynamics for effective particles, which recover the correct thermodynamic equilibrium and hydrodynamic modes of the system, even if the particles do not describe the atomic or molecular character of the solvent. Such approaches benefit from the knowledge of molecular simulation, and provide an appealing framework to understand the dynamics of the system. Among these approaches, Dissipative Particle Dynamics (DPD), introduced by Hoogerbrugge and Koelman 38 and formulated consistently by Español and Warren 39 models the solvent as point particles that interact through conservative as well as pairwise dissipative and random forces. DPD avoids hard-core exclusion effects, thus speeding up the computations to the extent that it becomes possible to reach hydrodynamically relevant time scales. Non-Hamiltonian interactions are proposed to ensure that proper hydrodynamic behaviour is recovered 40 .

Groot 41 proposed an extension of DPD which accounts for electrolyte dynamics. In order to avoid the collapse of oppositely charged point ions, local charges are smeared out around each DPD particle. A variant of the particle-particle particle mesh (PPPM) which solves the far field in real space, introduced to deal with systems with heterogeneities in the electrostatic properties of the system (such as the local permitivitty) 42 and general geometries, is then used to treat separately the near field and far field in-teractions between such charge distributions. Standard Ewald sums techniques can be adapted to study systems with soft potentials 43 , although PPPM has been shown to be more flexible and to provide better control of the numerical accuracy than Ewald 44 . This method allows a standard approach to deal with electrostatics, and the total force of the smeared charge distribution is then applied to the point particle. As the charge distribution in the system is affected by hydrodynamic flow, this approach provides a natural coupling between electrostatics and fluid motion. However, smearing the point charge into a Gaussian distribution introduces a new length scale in the system, which is usually comparable to the Bjerrum length, l B , whereas one would prefer to have a smearing on a range small compared with l B . The ability of DPD to treat hydrodynamic transport in complex geometries has been exploited to study the dynamics of a variety of charged, soft materials, including the self-assembly of charged surfactants 45 , polyelectrolyte suspensions 46 and brushes 47 as well as the effect that polyelectrolytes have on membrane properties 48 and their transport properties across bilipidic membranes 49 . DPD can be regarded as a momentum-conserving thermostat together with soft conservative forces. The use of alternative momentum-conserving thermostats with simpler structure, such as the Lowe-Andersen thermostat 50 , offers the possibility to combine molecular methods enforcing appropriate momentum conservation.

A combined study of a DPD fluid with resolved charged ions which interact through explicit electrostatic interactions has been carried out to compare the performance with a Lattice-Boltzmann (LB) method (see below) where also ions are treated explicitly 51 . A systematic procedure to map parameters between the two mesoscopic approaches is put forward, showing that the computational cost associated with the electrostatic coupling is equivalent in both approaches. The different interactions between the fluid and the solid wall in the two methods allow to assess the relevance of boundary conditions in electrokinetics.

Explicit solvent: Kinetic methods

Pure solvent

As an alternative to particle-based approaches, kinetic models have been introduced to exploit locality in dealing with the solvent dynamics. This usually involves including charged species in a pre-existing kinetic approach for a hydrodynamic system. A powerful tool to simulate hydrodynamics, which was later extended to deal with charged species, is the so-called Lattice-Boltzmann (LB) method. The LB approach follows the time-evolution of the one-particle density distribution of solvent molecules moving on a lattice using a discrete set of velocities. Accordingly, the model is discrete in time and space and the distribution function, f (r, v,t), evolves in two steps. In the first (collision) step, f relaxes locally to a prescribed equilibrium distribution, which can be understood as a low-velocity expansion of a Maxwell-Boltzmann distribution. In the second (propagation) step, f is advected to neighbouring nodes. Viscosity arises from the collision step. If the lattice satisfies certain minimal symmetry requirements, this kinetic model recovers hydrodynamic behaviour at long time and length scales [START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF] . Since all kinetic rules are local, this approach can easily deal with complex boundary shapes; in particular it becomes feasible to analyze the dynamics of complex fluids in arbitrary confining geometries. Colloidal suspensions can also be simulated using the same approach and one can enforce the appropriate time and length scale separations that characterize the hydrodynamics of such suspensions 53 .

Explicit microions

It is possible to couple the LB scheme to the motion of particles with sizes smaller than the lattice spacing. One then needs to follow the dynamics of these particles, which follow Newton's equation of motion together with a dynamic coupling which accounts for fluid friction. The particle experiences a friction force whenever its velocity differs from the fluid velocity from the nearby nodes. Since total momentum is conserved because there are no external forces involved in this coupling, the opposite of this force is also distributed among the neighbouring fluid nodes. To reach thermal equilibrium, it is necessary to add a stochastic force which compensates for the energy lost due to this dissipative coupling. This force follows a Gaussian process with zero mean and amplitude proportional to the friction coefficient. This approach has been successfully applied to study the dynamics of polymer and colloidal suspensions 54,55 and has been extended subsequently to include small charged particles 56 . These small particles can also be joined to form a larger macroion, such as a colloid. In this case the colloidal charge is located at the colloid center of mass while the small spheres which define its surface ensure, through Lennard-Jones potential interactions with the dissolved ions that these do not enter inside the colloid. This approach has been used to study the electrophoresis of colloidal suspensions 57 and also of polyelectroytes, showing the relevance of HI in their dynamical response 58 . Using the same approach, colloidal electrophoresis in the strong charge regime has been addressed 59 . Resolving individual ions has the advantage that the effect of electrostatic correlations on the dynamics of colloids and polymers can be captured. However, large colloids are difficult to reach; hence spatial separation between ionic and solvent sizes is lacking. In fact, the ion size sets the smallest relevant length scale in the simulations, which are then typically restricted to the nanometer scale.

Implicit microions

One way to reach larger length scales is to regard the collective, rather than individual, behavior of ions. In doing so, one loses in principle the possibility to account for electrostatic correlations due to the finite ionic size, but now the relevant length scales associated to the ionic densities are given by the Bjerrum and double-layer length scales, which are usually larger than the ionic molecular size, depending on the solvent under consideration (although for water at room temperature the Bjerrum length is comparable to an atomic size).

An alternative approach to address some of the problems mentioned above keeps the LB description of the fluid but, rather than following the dynamics of the individual ions, considers a density of anions and cations at each lattice node. Unlike the components in a LB binary mixture, these ion densities have no intrinsic velocities, and cannot therefore be regarded as one body distribution functions. One can set up dynamic rules discrete in space and time to reproduce the convection diffusion equations which characterize their motion. The coupling between the charged species and the fluid flow is twofold. On the one hand, the charge densities are advected with the local fluid flow; on the other hand, the local charge density results in a force and accelerates locally the fluid by adding a forcing term in to the LB model. Finally, the Poisson equation needs to be solved at each time step to determine the potential gradients and determine the local electric field. Real space solvers of the Poisson equation are not restricted to periodic boundary conditions. They allow in particular to account for general geometries and locally varying dielectric permittivites.

Originally, Warren 60 exploited moment propagation, a recursive scheme to diffuse a passive scalar in a LB scheme 61 , to study electro-osmosis of an electrolyte in a pore. Horbach and Frenkel 62 studied the foundations of the model and applied it to analyze the effect of the double layer width on colloidal sedimentation.

However, this approach suffers from "leakage" at fluid/solid boundaries. An improved implementation was subsequently developed 63,64 , which considered the diffusive fluxes along links joining nodes -the basic dynamic quantity which determines the temporal evolution of the charge densities. This approach has been applied to analyze the sedimentation of charged spheres 63 and disks 65 as well as to the dispersion of charged solutes in porous media 66 . The need to resolve the the charge density profile of the double layer on the lattice length scale in these approaches sets a stringent condition on the scales that can be captured numerically. In practice, it tends to limit the approach to the study of processes where the double layer is wide. It should be noted that the method of 63 relies on the knowledge of an approximate free-energy functional that depends on the charge densities. However, the method as such can also be applied to other non-ideal fluids 67 . In this context, the method described in the previous paragraph corresponds to an ionic mixture that behaves as an ideal gas, hence recovering Poisson-Boltzmann theory in equilibrium, although not restricted to its linearized, Debye-Hückel limit.

It is straightforward to extend the hybrid description of the electrolyte to model suspensions of charged colloids where colloids, unlike the small co-and counter-ions are modelled individually 68,69 . For example, the local structure of the LB method allows to couple the fluid to solid/liquid interfaces with general symmetries, favoring easy parallelization with good scalability 70 . The link-flux description for the temporal evolution of the charged densities also allows for a simple implementation of moving boundary condition without "charge leakage".

Fluid mechanics approaches

Instead of using a kinetic model to provide a more microscopic basis to describe the collective dynamic behaviour of the fluid, it is possible to solve numerically the Navier-Stokes equation, Eq. 2 on a lattice. In this context, the Fluid Particle Dynamics (FPD) 71 method regards explicit colloids as a fluid of large viscosity, with a viscosity profile determined by a tagged field centered around each colloid, φ . This profile (taken for convenience as a hyperbolic tangent) moves as the center of mass (c.o.m.) displaces in time. The c.o.m. velocity is obtained from the fluid evolution itself, weighted with the corresponding colloid's tagged field. Forces between colloids are also distributed following φ and added as appropriate contributions to f V in Eq. 2. Within this scheme, the solid/fluid interface is diffuse, but the computational effort is reduced to that of solving the Navier-Stokes equation on a lattice and convecting the particles according to their local fluid velocities. This approach has been generalized to account for electrolyte solutions 72 . The colloidal charge is localized on the colloid surface by distributing it according to ∇φ . Additional density fields for the electrolyte species, ρ ± are introduced, analogously to the hybrid method described in previous sections. The local charge provides an additional contribution to the force f V acting on the fluid, and the Nernst-Planck equation ( 1) is solved numerically on the lattice. This approach has been applied to study colloidal electrophoresis and its use to separate colloids with opposite charges.

In order to avoid distortions in the colloidal shape, FPD needs to operate at large viscosity contrasts, which comes at high numerical cost. In order to avoid numerical inaccuracies which result from compromises between the viscosity contrast and time steps in practical FPD implementations, it is possible to keep the same approach but regard φ as a frozen field and exclude the fluid from its interior. One then does not need to treat the colloid as a higher viscosity fluid, thus overcoming some of the previous limitations 73 . This alternative Smoother Profile (SP) method 74 has been applied to study the electrophoresis of colloidal suspensions as a function of colloid concentration and double layer width 75 . These methods offer complementary approaches to the use of standard techniques from fluid mechanics. The latter, traditionally exploited to analyze the fluid dynamics at macroscopic scales, have been also extended to consider electrolyte flows in generic microfluidic geometries 76 . An appropriate meshing of the space is required to deal with complex geometries.

Applications

The different approaches described in the previous section have been used to analyze a variety of physical systems. Below, we briefly discuss some representative examples to indicate the usefulness of the coarse-grained methods described and also to show the wealth of physical situations where charge dynamics, and the coupling with the underlying solvent, become relevant and should be accounted for properly.

Electrokinetics in colloidal suspensions

Electrophoresis

The understanding of electrophoresis is a classical subject in colloidal science [START_REF] Russel | Colloidal Dispersions[END_REF] . O'Brien and White performed a theoretical study of the electrophoretic mobility of a sphere in the Debye-Hückel regime for an arbitrary double layer [START_REF] O'brien | [END_REF] . However, for double layers that are outside the Smoluchowski (narrow double-layer) or Henry (very diffuse double layer) limits [START_REF] Russel | Colloidal Dispersions[END_REF] , expressions for the speed of electrophoresis are lacking and numerical simulations are needed to test the approximate theoretical predictions. Dünweg and Holm 4 analyzed the electrokinetics of charged colloids in the presence of counterions alone and in the low-salt regime, using lattice-Boltzmann simulations. In their approach, they make use of the so-called "raspberry method" 55 , where the colloid is regarded as a spherical shell covered by small particles, with radius smaller than the lattice spacing, which move as a rigid object. Dissolved ions are also modelled as small spheres that interact with the colloidal charge (located at the center of mass of the colloid). Electrostatic interactions are computed using a PPPM approach and both the particles that constitute the colloid and the moving ions interact with the surrounding fluid through momentum conserving, dissipative forces. This approach properly accounts for ionic correlations and short-range interactions, but it becomes prohibitively expensive for colloidal sizes that are much larger than the ionic size.

Combining experimental and numerical evidence, Ref. 4 has analyzed electrophoresis for high salt and low salt dispersions. The study shows that it is possible to analyze in a unified way the response of a suspension, regardless of the number of dissolved ions, in terms of the effective length κ -1 = (4πl B n c ) -1/2 , where n c = Z/V denotes the counterion number density. This length reduces to the standard Debye length when salt is present and n c corresponds to the total ion number density. The simulations show that when the double layer of different colloids start to overlap the contribution from the counterions to the electrophoretic mobility should be accounted for. Moreover, in the absence of added salt, the electrophoretic mobility goes through a maximum before leveling off at its limiting high-Z behaviour.

The numerical scheme is validated by comparing with experiments as a function of colloid volume fraction, Φ, and salt concentration. At low volume fractions the "evaporation" of the counterion layer around the colloid leads to a sharp increase of the electrophoretic mobility. On increasing volume fraction, on the other hand, the colloid acquires an effective charge that is not very sensitive to the bare colloidal charge, Z. This broad analysis identifies the main dependence of the electrophoretic mobility on the relevant system parameters, characterized in terms of the reduced electrophoretic mobility, µ red ≡ µ/µ H (Z = 1), where µ H = Ze/6πηR is the electrophoretic mobility in the Henry limit for a charged sphere of radius R in a fluid of viscosity η. The proposed scaling , µ red (Φ, Z, l B /R) µ red (Φ, Z), where Z = Zl B /R, reproduces the experimental results in all the regimes analyzed. It is possible to identify the regime where the response is linear in the applied electric field using the Green-Kubo linearresponse expression for the electrophoretic mobility 79 .

The same numerical approach, based on the 'raspberry model", has also been used to analyze the effect of a symmetric electrolyte on the electrophoretic mobility of a highly charged colloid 59 . In this regime a part of the electrolyte is strongly attached to the colloid, although the study shows that hydrodynamics reduces this strong coupling with respect to the equilibrium electrolyte structure. The dynamic distortion of the charge distribution close to he colloid is also sensitive to the ion properties. For asymmetric electrolytes a reversal of the electrophoretic mobility is reported 59 , in accordance with experimental observations.

The electrophoretic mobility of polyelectrolytes has been explored in the same framework 58,80 . These charged polymers are modelled as a chain of charged beads joined with non-linear springs. The electrophoretic mobility has a dependence on the number of monomers in accordance with the Zimm model, which displays a maximum as a function of the polymer charge analogous to observations obtained for the electrophoretic mobility of a highly charged colloid. Alternatively, an implicit solvent model has been used to analyze DNA electrophoresis 81 , showing that on increasing the salt concentration, the counterions cease to be dragged along with the chain. This change in the dynamical behaviour of the counterions induces an anomaly in the counterion excess conductivity, which also affects the polyelectrolyte mobility. For charged DNA segment diffusion the effect of salt concentration is weaker 82 . A minimum in the diffusivity is found at intermediate screening lengths as a result of a balance between the number of counterions dragged by the polymer and the effective hydrodynamic radius of the neutralizing cloud.

Studies of the electrophoretic mobility of charged suspensions have also been carried out using the SP approach 74 . The computed mobilities compare well with the predictions of O'Brien and White in the regime where the double layer is of the same order as the colloid radius (Rκ = 1 and 1/2) even for high colloidal charge, (up to Z = 500) and large field strengths. These findings are surprising as this is the regime where the linearized theory might be expected to fail. Simulations containing up to a few hundred colloids and Φ ≤ 0.6 show that the electrophoretic mobility decays initially as Φ 1/3 , as predicted by Ohshima 83 . Deviations from the (low-density) Ohshima prediction are only observed when the double layers of different colloids start to overlap.

The Fluid Particle Dynamics (FPD) method 71 has also been used to study colloidal electrophoresis to study the effect of double-layer overlap 72 . An enhancement of the velocity fluctuations with the electrolyte double layer width is observed; a simple theoretical argument predicts an increase of the relative velocity fluctuations that scales with (κR) -1/2 , in agreement with the numerical observations. Simulations of suspensions of oppositely charged colloids make it possible to quantify the efficiency of electrophoretic techniques to induce particle separation.

Sedimentation

A complementary electrokinetic geometry corresponds to charged colloid sedimentation. As opposed to the previous situation, there is now a net force acting on the suspension. As a result, the flow fields generated by the sedimenting colloids are longer ranged. Theoretical predictions in this case are restricted to the narrow double layer limit 84 . Numerical simulations using a Lattice-Boltzmann approach and disregarding individual ions have verified that Booth's predictions are recovered for narrow double layers and have extended the study of the terminal sedimentation of charged spherical colloids for intermediate and wide double layers 63 . It was found that the reduction in the sedimentation velocity is maximal for a double layer with a size comparable to that of the colloid, indicating that in this regime the distortion of the double layer is the largest. The results showed that for sedimentation finite-size effects are more relevant than for electrophoresis, due to the external force acting on the system. The effect of colloidal shape has also been addressed analyzing the sedimentation of charged disks 65 . Particle asymmetry leads to asymmetric mobilities. Double layer distortion and polarization is now different for transverse or longitudinal sedimentation, and for large applied fields this asymmetry generates fluid vortices at the corners of the colloidal particle, a feature absent with spherical colloids. At high charges, in the nonlinear regime, the sedimentation velocity becomes less sensitive to the colloidal bare charge. In 2 we display the velocity field around a highly charged disk sedimenting transversally. One can appreciate that the disk radius determines the characteristic size over which the fluid flow is coupled to the particle motion. Although barely visible, in the rear corners of the disk, the nearly stagnant flow develops vortices due to the large accumulation of electrostatic isopotential lines.The authors have analyzed how to identify an effective disk size which accounts for this charge renormalization.

Electrical conductivity

The Primitive Model described in section 2 has been the basis of modern theory of transport in electrolyte solutions. In particular, it allowed to go beyond the Debye-Hückel-Onsager limiting laws. Brownian dynamics simulations and integral equations based on the Mean Spherical Approximation (MSA) demonstrated the necessity to include hydrodynamic interactions between ions to describe the conductivity up to molar concentrations 20 . Complex mixtures of ions, including micellar suspensions 22 have also been investigated at this level of description (explicit microions and HI tensor). It was found, for example, that HI tend to increase the self-diffusion coefficient of macroions 85 .

Dahirel et al. have also used a coarse-graining procedure to average out the degrees of freedom of the microions 86 . BD simulations with explicit microions are used to compute 1) the PMF between two macroions and 2) the effective diffusion coefficients for a single macroion infinitely diluted in a bath of microions. These two ingredients are then used in BD simulations of an effective one-component system of macroions, to compute the self-diffusion coefficient of macroions at finite macroions concentration. It is important to consider not only the static effect of the microions (PMF), but also their dynamical influence on the macroions, because the electrostatic friction considerably lowers the diffusion coefficient of the macroions, even when the concentration of the latter is small. This approach has recently been extended to include the effect of HI, and investigate the self-diffusion and mutual diffusion of macroions, their electrophoretic mobility and the conductivity of colloidal suspensions 87 .

Acoustophoresis

Another interesting coupling between solvent and ionic flows is the acoustophoretic effect, predicted by Debye in 1933 88 . It occurs when an ultrasonic wave is applied to an electrolyte solution: Local charge separation results in the generation of an electric field. While Debye's treatment was based on point particles, later studies at the MSA level allowed to estimate the effect of the finite ionic size and of the HI on the so-called ionic vibration potential (IVP) in electrolytes and colloidal vibration potential (CVP) in colloidal suspensions 89 . It is very useful to combine data from acoustophoresis and conductivity, because the relative contribution of the micro and macroions is not the same in both experiments. There is however room for improvement in the models, and numerical simulations with one of the mesoscopic methods described in the present paper will certainly provide valuable insight in this field. Such simulations might also provide the basis for a better understanding of the seismo-electric effect used for underground exploration 90 .

DNA translocation through nanopores

Recent experiments have analyzed the translocation of DNA through nano-holes in a membrane. The use of the lattice-Boltzmann method has allowed to identify the different contributions to the mobility of the DNA 91 . It was shown that such a description is sufficient to capture the low charge and high charge regimes of DNA, and to recover Manning condensation in equilibrium. The authors first regarded DNA as a uniformly charged rigid rod. Such an assumption allowed them to identify the role of the cylinder-to-pore cross section when translocation is driven by an applied electric field perpendicular to the perforated membrane. The dynamic approach underlying the lattice Boltzmann method allowed the authors to identify the different contributions to the electric current and to distinguish steady-state from transient effects, something that was problematic in earlier atomistic simulations 92 and in experiments 93 . At high salt concentration, there is a decrease in the ionic current upon translocation, which is progressively smaller as the overall DNA charge increases. However, the situation changes at low salt concentrations, when the double layer becomes of the same order as the free space between the DNA molecule and the edge of the translocation channel. As shown in 3.a the electric current decreases when the DNA enters the pore if the double layer is narrower than the pore size, while for a lower salt concentration an increase is observed, as displayed in fig. 3.b. In both cases the effect becomes more marked as the DNA charge increases. Such a dependence of the diffusive ionic current arises as a result of the competition between the decrease in available cross section, which leads to a corresponding decay in the current, and the local deformations in the electrostatic field around the pore which favors the transport of ions across the pore. In general, the net current depends on the thickness of the double layer, which tends to increase the electric current, and on the partial blockage of the hole by the DNA, which tends to decrease the net electric current. 91

with permission from ACS).

The flexibility of the approach allowed the authors to consider the effect of distributing the charge along the DNA molecule following a helical distribution in the present of a charged "probe" in the channel wall. In this case, the current through the pore follows the charge distribution, with a period equal to half the translocated DNA pitch. The amplitude of the signal can be modulated with the applied field, and its modulation is optimized for double layer widths smaller than the pitch and of the same order than the charged patch. The translocation of a flexible polymer chain under the action of an applied electric field has also been addressed 94,95 , assuming that the effect of the external force can be accounted for as an effective force localized at the hole. A complete analysis including flexibility and full electrokinetic coupling is still missing.

Solid-fluid interactions

The relevance of molecular structure and electrolyte specificity on the dynamical response to applied forces has been investigated using extensive molecular dynamics simulations. Using an implicit solvent approach, with fully resolved ions, it has been possible to study the electro-osmosis near an atomically resolved solid substrate 24 . The numerical results showed that there is a connection between the dynamics of ions in the presence of a solid wall and the traditional Helmholtz-Smoluchowski theory (HS) which predicts the liquid velocity profile, u(z) = u(0)

+ ε 0 ε r E η [ψ(z) -ψ(0)],
under the action of an applied electric field, E parallel to the solid interface, where η stands for the liquid viscosity and ψ is the electrostatic potential, which varies only within the double-layer. Moreover, deviations from the standard theory arise from surface heterogeneities, because charge correlations on a homogeneous solid surface does not invalidate HS. The analysis of the effective friction that a flowing electrolyte exerts on different types of solid substrates indicated that the interaction can be dominated by the hydrodynamic friction for moderate surface charge and/or low valency electrolytes, while for highly charged substrate and/or multivalent ions the direct electrostatic friction dominates 24 . These results generalize earlier ones obtained from Brownian dynamics simulations, in the absence of solvent, which had shown that the electrofriction between counterions and the ions constituting the charged surface have a pronounced effect in the mobility reduction of the former. The immobilization of a fraction of the counterions leads to a dynamic decrease of the effective charge of the solid substrate 96 .

The interplay between hydrodynamic, excluded-volume and electric friction induced by surface corrugations leads to effective dynamical charges smaller than the equilibrium ones 97 , in accordance with experimental observations which find discrepancies between the equilibrium charge measured on a solid surface through titration and the one derived from the zeta potential measured for example in electro-osmotic experiments. The dispersion interaction of the solvent and ions with the substrate has also an important impact on the electrokinetic response of the electrolyte 98 . In this respect, molecular dynamics simulations have also been useful to disentangle the effect of the solvent wettability and hence, the effective slip of the fluid at the solid surface. Slip induced by the substrate hydrophobicity induces a larger zeta potential when a pressure gradient is applied on an electrolyte solution confined in a charged slit. The variations of the effective zeta potential with the ionic strength differs qualitatively, depending on whether the substrate is hydrophilic or hydrophobic, and hence correlates with the effective slip length (which increases with the substrate hydrophobicity). While for hydrophilic substrates the zeta potential can be easily related to the substrate charge, this is not the case on hydrophobic walls, leading to increases in its value of one order of magnitude for moderate slip lengths 99 . As a result, slip can significantly enhance electrokinetic phenomena and that such a behaviour can be tailored both by modifying the substrate hydrophobicity and geometry 19 . The possibility to induce giant electroosmotic and electrophoretic effects, as well as other coupled phoretic phenomena 100 is potentially interesting for microfluidic applications, and experimental results indicate that such applications may be viable.

Further molecular dynamics studies have shed light on the anomalous electrokinetic phenomena induced by the molecular structure of the electrolyte layer close to the solid surface. It has been shown that if the degree of affinity between the ions and the solid substrate is asymmetric between anions and cations, the ionic structure is subsequently altered leading to changes in the effective ζ -potential which affects the observed electroosmotic flow 101 . In narrow channels, ion immobilization arising from the molecular interaction with the solid substrate has also been reported to lead under some circumstances to electro-osmotic flow reversal 102 .

Electrokinetics at liquid/liquid interfaces

Solid-liquid interfaces are not the only ones relevant for electrokinetics. Experimental evidence of charge separation at oil/water interfaces has recently been reported, including data obtained from scattering techniques that probe the interfacial ionic profiles on the microscopic scale 103 . This spontaneous charge separation is driven by the differential solvation properties of the ions, with e.g. a hydrophilic cation and a hydrophobic anion. This mechanism was also invoked to explain the metastability of crystals of µm-sized water droplets in oil: The charge separation at the interface results in an effective charge for the droplets, and a long-range repulsion across the low permittivity oil phase [START_REF] Leunissen | Proc. Nat. Acad. Sci[END_REF] . This separation was then exploited tune effective interaction between colloids 104 .

From a dynamical point of view, the electrophoretic mobility of droplets and bubbles has attracted much attention, in particular because the polarity of the interface deduced from these macroscopic experiments, i.e. a negative surface for air bubbles in water, seems in contradiction with microscopic surface-sensitive techniques and simulations 105 . Electric fields can be used to induce the motion of liquid droplets, e.g. for lab-on-chip manipulations 106 , and the coalescence of oppositely charged droplets 107 . In the latter case, large electric fields cause the droplets to bounce off each other, because of the dynamic coupling between electrostatic and capillary forces 108 .

Simulating such charged fluid-fluid interfaces is particularly challenging. It requires the description of the two solvents and several types of ions, and it should account for several features: the immiscibility of the solvents, the species-dependent affinity of the ions for both solvents, and electrostatics. The hydrodynamic behaviour of these systems is governed by the coupled capillary, electrostatic and solvation forces. The hybrid Lattice-Boltzmann / link-flux approach (see section 3.4.3) introduced for electrolyte solutions has been recently extended to tackle these complex systems 67 . The fluid is treated as a mixture characterized by a local composition φ = ρ o -ρ w ρ o +ρ w ∈ [-1, 1], with ρ o and ρ w the local oil and water densities. Phase separation results from the free energy of mixing, assumed to be of the Landau-Ginzburg form:

F mix = dr - 1 2 Bφ 2 + 1 4 Bφ 4 + 1 2 K(∇φ ) 2 . (3)
The first two terms correspond to the bulk phase behaviour, with minima for φ = ±1, while the last reflects the cost of sustaining interfaces. Ions contribute to the free energy of the system with an ideal, a solvation and an electrostatic term:

F ions = dr ∑ α=± ρ α (r) k B T (ln ρ α (r) -1) -µ α +V solv α (r) + z α e 2 ψ(r) (4)
where z ± = ±1 is the valency of the ions, µ α is a reference chemical potential and the electrostatic potential ψ is solution of the Poisson equation ∇

• [ε(r)∇ψ(r)] = -[ρ + (r) -ρ -(r)] e.
Both the solvation potentials V solv α (r) and permittivity ε(r) depend on the local composition φ (r) of the fluid. The derivatives of the free energy with respect to ρ ± and φ provide the local chemical potentials, whose gradients are the source of ionic and solvents fluxes (from the generalized Fick's law and the Cahn-Hilliard equation), and of thermodynamic forces, according to the Gibbs-Duhem equal-ity: f th V = -∇P = φ ∇µ φ + ρ + ∇µ ex + + ρ -∇µ ex -. These forces accelerate the fluid, which in turn contributes convective fluxes changing the local composition.

This free-energy based approach allows one to recover the analytical predictions obtained for a system in equilibrium in the linearized regime for a sharp, planar oil/water interface by Onuki 109 and van Roij et al. 110 . The difference in ionic solvation free energies ∆µ ± = µ o ±µ w ± results in the partitioning of the salt, with a salt concentration ratio ρ o s /ρ w s = e -β ∆µ av where ∆µ av = (∆µ + + ∆µ -)/2, and in a potential difference across the interface, the Donnan potential ψ D = ψ oψ w = (∆µ --∆µ + )/2e. The interface behaves qualitatively as a planar capacitor subject to a potential difference ψ D . Simulations show that the linearized Poisson-Boltzmann (LPB) theory provides a good description of the surface charge density and the ionic adsorption even in the moderately non-linear regime (β eψ D ≥ 1). This somewhat surprising result could be traced back to a compensation of errors. In this regime, the analytical LPB theory predicts smoother ionic profiles than the non-linear PB result for a sharp interface. However, the sharp interface model implies an unphysical jump in the ionic densities and overestimates the surface charge density with respect to the simulations with a smooth interface. As shown in 4 the ionic profiles obtained by re-exponentiating the LPB electrostatic potential is in good agreement with the simulation results, except in the very close vicinity of the interface. Figure 4: Equilibrium ionic profiles at a planar oil/water interface, in the case of hydrophobic anions and hydrophilic cations (for the details of the system parameters, see Ref. 67 ). ρ w corresponds to the ionic density in the bulk phases and ∆x stands for the lattice spacing. Simulation results (symbols) for cations and anions are compared to an analytical result for an infinitely thin interface (see text and Ref. 67 ). The insert shows the local composition φ and electrostatic potential ψ. This simulation method is however not limited to the linearized regime, nor to equilibrium. It can be considered as a numerical solver for the coupled Navier-Stokes / time-dependent density functional theory (DFT) equations, and opens the way for the simulation of the complex phenomena occurring at a charged liquid-liquid interface described above. Preliminary results indicate that it is possible to describe both the transient and stationary electro-osmotic flow generated by an electric field applied parallel to a planar charged oil/water interface.There is then no net external force, but the two fluids are accelerated in opposite directions, and the steady-state velocity profile results from the balance between electric forces and viscous dissipation at the interface.

Dispersion of charged tracers by an electro-osmotic flow

The diffusion of tracers in inhomogeneous flows results in their dispersion: A tracer pulse broadens under the combined effect of advection by the flow and diffusion in the direction perpendicular to the streamlines. Electro-osmotic flows (EOF) are inhomogeneous mainly in the vicinity of the surface, as can immediately be inferred from the Helmholtz-Smoluchowski equation. This feature explains the less pronounced dispersion of tracer pulses in EOF compared to Poiseuille flows generated by pressure gradients, and their usefulness for applications to separation science and analytical chemistry. The approach pioneered by Taylor and Aris 111,112 makes it possible to quantify electrolyte separation in electro-osmotic flows, as in the case of Taylor dispersion in microfluidic environments 113,114 . Such an understanding provides insight into the competing mechanisms that give rise to electrolyte dispersion in flows in channels 115 . Particularly interesting is the fact that tracers of different charge are dispersed differently by the same EOF, because they do not probe the same part of the flow. The dispersion of charged tracers by an EOF in a salt-free slit pore has been recently addressed using the hybrid LB approach presented in section 3.4.3, and the moment propagation method 61 to compute the Velocity Autocorrelation Function (VACF), Z k (t), of tracers with valency z k 66 . From the VACF, one computes the time-dependent diffusion coefficient in the direction perpendicular to the flow, D ⊥ k (t) = t 0 dt Z ⊥ k (t ), which describes how tracers probe the different parts of the slit pore. This function decays faster for co-ions than for neutral tracers, because the characteristic decay time scales as L 2 /D with L the width of the slab accessible to the tracer. The latter is smaller for co-ions than neutral species, since co-ions are repelled from the surfaces and confined to a smaller volume. Counterions also exhibit a faster decay at short times, since they are confined to the diffuse layer close to the surface, but a slower decay is observed at longer times. This corresponds to the motion from one diffuse layer to the one on the opposite surface, which requires to overcome an electrostatic barrier and is thus slower than the diffusion of neutral species over the same distance.

When the observation time is large compared to the characteristic time to explore the direction normal to the charged surfaces (and thus to the flow), the relevant quantities are the average velocity of each tracer in the direction of the flow and the dispersion coefficient. They are obtained from the VACF as vk = lim t→∞ Z k (t) and

D k = ∞ 0 dt Z k (t) -Z k (∞)
, respectively. The simulations, in agreement with approximate analytical results, indicate that the average velocity is larger (resp. smaller) for co-ions (resp. counterions) than for neutral tracers. This may seem paradoxical, since the direct effect of an electric field on the co-ions is to drive them in the direction opposite to the EOF, whose direction corresponds to the electric force on the excess of counterions. The explanation lies in the fact that for small ions in water, the electrophoretic velocity β Dz k E is small compared to the EOF whose typical velocity is u EOF = eE/2πηl B , with η the viscosity and l B the Bjerrum length. Then, the EOF velocity is larger in the center of slit, where most co-ions are found and smaller near the surfaces, where most counterions are located. This explains the behaviour of these species compared to neutral ones, which probe the whole slit with the same weight. The different localization also explains why the dispersion coefficient is larger (resp. smaller) for counterions (resp. co-ions) than for neutral species, since the flow is more inhomogeneous near the surfaces than in the center of the slit. Using this approach, it has also been possible to investigate quantitatively the dispersion coefficient, which 67 ). The effective diffusion coefficient of tracers with different charge (neutral •, cationic 2 and anionic tracers) is reported normalized by the molecular diffusion coefficient, as a function of the ratio between the electrolyte Debye length, κ -1 b , and the sphere radius, R. The inset shows the characteristic time to explore the pores for each tracer.

was found to vary as D k = D k 1 + f (αL, z k ) × Pe 2 , with Pe = u EOF L/D k the Péclet number, D k the molecular diffusion coefficient of species k, and f a function of the tracer charge z k , of the distance between surface L and of a screening length α -1 analoguous to the Debye length. Such variations can be seen in 5.a; the inset of this figure also reports the dispersion factor f (αL, z k ).

Moment propagation has also been used to investigate the diffusion of charged tracers through charged porous media 67 . 5.b reports the effective diffusion coefficient of charged tracers through the pores of a compact FCC lattice of negatively charged spheres. Results are given as a function of the ionic stength of the saturating fluid, expressed in terms of the ratio between the Debye length and the radius of the spheres. These pore-scale simulations allowed to relate the effective diffusion coefficient to the microscopic pathways followed by tracers of different charge. Cations are confined to the diffuse layers in the vicinity of the surfaces, and explore a small, connected fraction of the pores, while anions are mainly located near the center of large pores and must overcome an electrostatic barrier to go from one pore to the next. This interpretation is also confirmed by the characteristic time scale required for each tracer to explore its accessible porosity (see insert of 5b).

Outlook

In this article we have focused on the different modelling approaches to tackle dynamical phenomena in charged heterogeneous materials. The study of electrokinetics in these materials poses a number of challenges due to the different length and time scales involved The methods described are complementary in their goals, and provide a set of simulation tools to address a variety of phenomena, ranging from the fundamental understanding of boundary conditions and the role of ion correlations to the analysis of long-range pattern formation mediated by electric and flow coupling. These new methods, which capture the coupling of solvent flow with electrolyte diffusion and migration have proven a valuable tool to understand the interplay between these different phenomena. In the selected examples, we have shown that such a dynamic coupling is a key factor to understand some of the features observed experimentally, and that a proper handling of these phenomena is necessary to describe the dynamics in heterogeneous charged materials. The insight gained and the methods developed in recent years provide a useful starting point to explore a wide range of challenging problems. This includes for example the relevance of electrokinetics in cellular processes, such as the effects of charge transport in biological membranes and their role in the stability of such structures, or the electro-osmotic effects related to transport through ion channels.

Acknowledgments

Figure 2 :

 2 Figure 2: Flow field for a sedimenting disk of radius R with valency Z = 300 and aspect ratio 10 in an electrolyte solution with ionic strength κR = 1. ( Figure reproduced from Ref. 63 with permission from AIP).

Figure 3 :

 3 Figure 3: Diffusive ionic current, I diff normalized by the electrolyte current, I 0 as a charged rod of surface charge density σ and length l c passes through a membrane pore of diameter d larger than the rod section. The position of the rod is characterized by the distance x of its center of mass to the membrane surface. a) high ionic strength, κd = 1; b) low ionic strength, κd = 0.02. ( Figure reproduced from Ref. 91 with permission from ACS).

Figure 5 :

 5 Figure 5: (a) Dispersion coefficient of charged tracers by an electro-osmotic flow in a salt-free, negatively charged slit of width L, under and applied electric field E parallel to the surfaces. Results are given for anions (+), neutral tracers (2) and cations (•), for two surface charge densities, αL. The inset shows the dispersion factor f defined from the quadratic dependence of D k with E (see text), as a function of the surface charge density. Lines are results of an approximate analytical theory. ( Figure reproduced from Ref.66 with permission of IOP). (b) Diffusion of charged tracers in a modelporous medium consisting of a compact FCC lattice of negatively charged spheres ( see Ref.67 ). The effective diffusion coefficient of tracers with different charge (neutral •, cationic 2 and anionic tracers) is reported normalized by the molecular diffusion coefficient, as a function of the ratio between the electrolyte Debye length, κ -1 b , and the sphere radius, R. The inset shows the characteristic time to explore the pores for each tracer.

I

  .P. acknowledges financial support from MICINN (Spain) under project FIS2008-04386 and DURSI (2009SGR634). B.R. acknowledges financial support from AN-DRA (France), and from ANR under project SIMISOL (ANR-09-SYSC-012). D.F. acknowledges financial support from the Royal Society of London (Wolfson Merit Award) and from the ERC (Advanced Grant agreement 227758). I.P. and B.R. thank the Chemistry Department of the University of Cambridge for hospitality.

Table 1 :

 1 Summary of the levels of descriptions discussed in the present article. The meaning of the acronyms can be found in the corresponding sections.

	Solvent