## Abstract

We compute the electrostatic potential at the surface, or zeta potential *ζ*, of a charged particle embedded in a colloidal suspension using a hybrid mesoscopic model. We show that, for weakly perturbing electric fields, the value of *ζ* obtained at steady state during electrophoresis is statistically indistinguishable from *ζ* in thermodynamic equilibrium. We quantify the effect of counter-ion concentration on *ζ*. We also evaluate the relevance of the lattice resolution for the calculation of *ζ* and discuss how to identify the effective electrostatic radius.

## 1. Introduction

A quantitative understanding of the electrostatic interactions between charged macro-ions in solution is fundamental to comprehending a plethora of physical phenomena spanning from biology to material science, for example, the macroscopic and rheological properties of colloidal suspensions [1,2]. The details of such interactions are difficult to capture, as the effective interactions are determined by the interplay between the different components dissolved in solution (macro-ions, counter-ions, salt ions) and the solvent dielectric response [3]. In addition, when the system is driven out of equilibrium, for example, by an external electric field causing electrophoretic flow, hydrodynamic interactions between the solvent and solute species must also be included and can alter the equilibrium electrostatic interactions.

The simplest yet very useful model to describe electrolyte solutions is the Poisson–Boltzmann (PB) approach [1,4], which builds on a continuous description of the electrolyte, characterized in terms of the anion and cation local densities. Within this framework, it is possible to add charged macroscopic objects, with fixed surface charges, by accordingly changing the electrostatic boundary conditions of the PB equations. Notwithstanding that PB is a mean-field theory that defines ions as point-like and therefore neglects excluded-volume effects and spatial correlation due to their finite size, it has been successfully adopted to describe various physical systems, for example, to derive the electrostatic part of the Derjaguin–Landau–Verwey–Overbeek (DLVO) interparticle potential [1], which explains the interactions of weakly charged particles in solution at low volume fraction.

The zeta potential, *ζ*, defined as the electrostatic potential at the colloid surface, or slipping plane, plays a central role for charged colloidal dispersions, as it indirectly provides an estimate of the magnitude of the electric field between particles that results from the combined effect of ionized charged groups sitting on the particle, counter-ions released into solution and salt ions. Theoretically, Ohshima *et al.* [5] obtained an exact analytic expression relating the surface charge density to the zeta potential for an infinitely dilute spherical colloidal suspension by carefully approximating the nonlinear PB equation. This expression is of limited use, since the infinite dilution regime is difficult to achieve. Moreover, experimentally, *ζ* is normally derived out of equilibrium from electrophoretic mobility measurements. Since, *a priori*, the values of *ζ* in equilibrium and at steady state will differ, an accurate description of hydrodynamics must be included in order to derive the relationship between *ζ* and the colloid mobility.

Analytical results for the electrophoretic flow at finite dilution cannot, in general, be obtained. Computationally, different models have been used to describe electrophoresis [6–10]. Lobaskin *et al.* [8] and Dünweg *et al.* [9], using a lattice Boltzmann (LB) [11] solver coupled to molecular dynamics of ions, calculated the particle mobility for low to zero salt concentration, explicitly accounting for finite ion size, but without providing results for *ζ*. Kim *et al.* [10] show agreement with Ohshima’s results, but only for weakly charged particles at low salt concentration. Cell models [12] are also employed to describe electrokinetics and to derive mobility versus *ζ* curves; however, they rely on somewhat arbitrary boundary conditions for both electrostatics and hydrodynamics.

In this paper, we compute *ζ* using a mesoscopic model that couples a Navier–Stokes solver to the convection and diffusion of ions on a lattice. We will first investigate the difference between *ζ* values at thermodynamic equilibrium and at electrophoretic steady state. We then accurately analyse the dependence of *ζ* on colloid volume fraction, *Φ*, and salt and counter-ion concentrations, and compare them with predictions from PB at infinite dilution derived by Ohshima *et al.* [5]. Finally, we will address how to identify the electrostatic radius of a colloid on a lattice, as has been done previously for the hydrodynamic radius [13,14]. In the following section, we briefly introduce our model, the reference equations and the results that can be obtained using the PB equation for infinite dilution and at finite volume fractions. We present in §3 results for *ζ* as a function of the volume fraction, salt concentration and lattice refinement. Finally, we draw our main conclusions in §4.

## 2. Model

### (a) Electrophoretic flow

In order to derive a correct relation between the mobility and the zeta potential at steady state during electrophoresis, our model includes hydrodynamic interactions and electric forces due to local charge densities,
2.1
2.2
2.3
Here *D*_{k} and *z*_{k}, *k*=+,−, are the diffusivities and valences of positive and negative ions, *ρ*,** v**,

*p*

_{id}and

*η*correspond to the solvent density, velocity, ideal pressure and shear viscosity,

*e*is the electron charge,

*β*

^{−1}=

*k*

_{B}

*T*is the Boltzmann factor,

*φ*is the electrostatic potential and

*ϵ*is the solvent permittivity. Equation (2.1) expresses ion mass conservation during diffusion and advection, coupling ion dynamics to solvent motion. This is, in turn, described by the Navier–Stokes equation for a viscous fluid, which couples solvent dynamics to electrostatic forces due to local charge density (equation (2.2)). Finally, the Poisson equation enforces the electrostatic coupling between the charged species and the embedded solid particles. We solve these electrokinetics equations by combining an LB approach for the momentum dynamics with a numerical solver for the discretized convection–diffusion dynamics (equation (2.1)) based on link fluxes [15]. Each lattice site is labelled as solid or fluid, and we consider here a single spherical object of radius

*r*=

*a*embedded in a cubic lattice of volume

*L*

^{3}with periodic boundary conditions. Suspended particles are therefore resolved and interact with the neighbouring fluid through bounce-back [13,14]. The electrostatic potential is computed using a successive over-relaxation (SOR) scheme [16]. This method has been successfully applied to analyse different dynamical processes involving suspensions of charged objects [17–20].

### (b) Poisson–Boltzmann electrostatics

At equilibrium, using a mean-field approach that neglects excluded-volume effects and correlation between ions, the PB equation reads [4]
2.4
where a symmetric electrolyte *z*_{+}=*z*_{−}=*z* has been used for simplicity and *ρ*_{0}, the uniform macroscopic counter-ion and co-ion number concentration far from the colloid *ρ*_{+}=*ρ*_{−}=*ρ*_{0}, is assumed to be equal to that of an electrolyte reservoir with dissolved salt *c*_{salt}=2*ρ*_{0}.

An analytical solution of equation (2.4) for a single particle is not available, as the PB equation can be solved analytically only for a few symmetrical configurations. Ohshima *et al.* [5] solved equation (2.4) around a spherical colloid of charge *Q* using a perturbative approach and derived a quasi-exact analytical expression for *ζ*,
2.5
where is the inverse of the Debye length, *l*_{B}=*e*^{2}/4*πϵk*_{B}*T* is the Bjerrum length and is the adimensional zeta potential.

Useful global information can be obtained using the Debye–Hückel approximation, i.e. linearizing equation (2.4). When ,
2.6
from which one can derive the electrostatic potential around a charged spherical colloid as
2.7
In the linear regime, a charged particle in the presence of salt develops a screened, Yukawa-like electrostatic potential, with a characteristic length scale of *λ*_{D} and strength given by the zeta potential, *ζ*.

Equation (2.4) is strictly valid when the amounts of positive and negative charges dissolved in solution are equal, i.e. when the counter-ion concentration is negligible. Experimentally, this can be obtained by separating the colloidal particles from the bulk electrolyte solution with a semi-permeable membrane. For the experimental set-ups with non-zero concentration of counter-ions, equation (2.6) becomes
2.8
where *ρ*_{s0} and *ρ*_{c0} are the reference salt and counter-ion concentrations in the reservoir. Equation (2.8) stresses the fact that, especially at low salt concentration *c*_{s} (low *ρ*_{s0}) and high volume fraction *ϕ*=4*πa*^{3}/*L*^{3} (high *ρ*_{c0}), a deviation from the theoretical Debye–Hückel results can be observed.

Our aim is to measure *ζ* for different experimental regimes, analyse its behaviour beyond the linear approximation (when ) and understand how to overcome the intrinsic inaccuracies in the location of the electrostatic radius of any lattice model. In the LB approach, we will refer to *ζ* as the average electrostatic potential between lattice point pairs (boundary links) *l* that link a colloid and a fluid node calculated halfway between the two nodes using the linear approximation,
2.9
and we will analyse when such an assumption is representative of the colloid *ζ*.

## 3. Results

We have performed simulations for a single spherical particle of radius *a* embedded in cubic boxes with periodic boundary conditions, also varying the concentration of added salt. In order to avoid finite-size effects, for a given salt concentration we use a lattice length *L* at least three times the corresponding Debye length, *λ*_{D}. We have distributed the colloid charge *Q* uniformly on the lattice sites that the colloid fills, and choose *Q* to correspond to the predicted value of Oshima, according to equation (2.5). We will analyse both a magnitude of *Q* that lies in the linear regime, , and one for which the colloid is within the nonlinear regime, . We finally set the external electric field to *E*=0.01, well within the linear regime *E*≪*ζ*/*λ*_{D}.

The values for *ζ* obtained from electrophoresis are in principle different from those computed in equilibrium. This is analogous to different measurements for the size of a solvated particle, at equilibrium or at hydrodynamic steady state (hydrodynamic radius). We have therefore computed at thermodynamic equilibrium, i.e. with equilibrated charges and fluid at rest, and dynamically, when a steady state due to the external electric field is reached. In all simulations, our results show that the differences between the two values of *ζ* are smaller than the statistical precision. This similarity is expected because of the small magnitude of the perturbing electric field and colloid Péclet number, *Pe*=*av*_{coll}/*D*_{±}, where *v*_{coll} is the velocity of the colloid at steady state. For *Pe*≥1, the deformation of the salt cloud around the colloid induces significant deviations in the electrophoretic mobility [19]; we can hence envisage a corresponding departure of the dynamic *ζ* from its equilibrium counterpart. Typically, as for simulations here, *Pe*≪1 and particle mobility does not depend on *Pe*, hence the measured *ζ* values at equilibrium and steady state are statistically indistinguishable.

As the only theoretical result available, , is valid for a single particle in equilibrium at infinite dilution, we will compare the simulation results with equation (2.5) to assess the role of the counter-ion concentration at finite volume fraction. In addition, in order to analyse the extent to which the midpoint stands as a proper definition for the electrostatic radius, we run a pair of simulations for each colloid charge, volume fraction and salt concentration used, changing the resolution of the colloid in the lattice from low, using colloid radius *a*=4.5, to high, a *a*=9.0, also doubling the corresponding box sizes (e.g. from *L*=30 at low resolution to *L*=60 at high resolution for the highest volume fraction).

Figure 1 displays for a weakly charged colloid as a function of the system size and salt concentration, where the left (right) panel corresponds to a low (high) colloid resolution. The top panels show that the deviation of from Oshima’s prediction for varying *Φ* appears indirectly only through the corresponding change in the co- and counter-ion concentrations because the deviation decreases when the salt concentration is increased. Therefore, according to equation (2.8), we can gain understanding by analysing the dependence of on the ratio between the concentration of counter-ions released by the particle, *c*_{cions} (which decreases as *Φ* increases), and the concentration of salt dissolved in the system, *c*_{salt} (which is independent of *Φ*). The bottom panels confirm an agreement within 10 per cent between and for . In addition, comparing the left and right panels, we conclude that in the linear regime the lattice resolution does not play a significant role, as for a given salt concentration does not differ significantly and the use of a more computationally expensive refined lattice only reduces the error bars. In addition, we note from the bottom-right panel that discretization effects become more relevant when increasing the salt concentration. This sensitivity is associated to the reduction of *λ*_{D} and the corresponding loss of resolution in the electrostatic potential around the colloid.

Figure 2, organized analogously to figure 1, shows results for a strongly charged colloid, , for which the linearized Debye–Hückel approximation does not hold. We observe that the dependence of on *Φ* enters again indirectly through the relative changes in *c*_{cions}, and that the convergence towards can only be expected for high *c*_{salt}. However, as opposed to the linearized regime, adding salt can lead to a 20 per cent overestimation of for the highest salt concentration analysed (*ka*=4.5), especially when using a low-resolution lattice (left panels). In this strong coupling regime, the electrostatic potential decays faster than *λ*_{D} [1]. This strong nonlinear behaviour invalidates the midpoint as the natural choice for the electrostatic radius, which develops a dependence on the system parameters even when the effects of finite *Φ* and *c*_{cions} are negligible.

Figure 3 displays the measured electrostatic potential at the midpoint between a solid and fluid node, , together with the two extreme situations where the electrostatic potential is averaged over the solid and the fluid nodes corresponding to the boundary links. One can observe, as expected, that always lies in between the other two estimates of . Comparing the left and right panels, we conclude that the increase in resolution does not significantly affect , while decreasing the inaccuracy in the estimate of the limiting values for . Generically, the better resolved the colloidal particle (and the corresponding decay of the electrostatic potential, *φ*), the less spreading is seen between the different electrostatic potential estimates. When we increase *c*_{cions}, Oshima’s prediction no longer holds and equation (2.5) cannot be used to identify the electrostatic radius, as can be appreciated in the rightmost values of in the bottom panels. However, the weak dependence of *φ* in the boundary fluid nodes suggests that, when is no longer valid, we can still identify an electrostatic radius slightly larger than the one associated to . Analogously to the hydrodynamic radius, an ad hoc calibration of the electrostatic radius at all salt concentrations must be done to perform quantitative studies of the electrokinetics of colloidal suspensions with LB.

## 4. Conclusions

In this paper, we have analysed how to determine the electrostatic potential at the surface of a particle, or zeta potential *ζ*, for a colloidal suspension. To this end, we have made use of a hybrid mesoscopic model that couples a discrete lattice formulation of Boltzmann’s kinetic equation for the solvent to a discrete solution of the convection–diffusion equation for the charged ion species dissolved in the fluid, which implies treating counter- and salt ions as scalar fields at the PB level. We have found that, for weakly perturbing electric fields, it is not possible to distinguish between the equilibrium and dynamic zeta potentials. We expect that differences will arise when the deformation of the charge layer around the colloid becomes significant, a scenario that can be achieved, for example, when the Péclet number of the dissolved ions is not negligible.

In particular, the comparison with Oshima’s expression for at infinite dilution has allowed us to carry out quantitative checks to validate the code performance. We have seen that the counter-ion concentration has a significant effect on , leading to a decrease in its magnitude in both the linear and nonlinear regimes away from Oshima’s result. We have shown that the ratio between counter-ion and salt concentrations controls the departure from Oshima’s prediction, and that its prediction works reasonably well for . We have also assessed the relevance of the lattice resolution and have quantified its effects. We have seen that the mean of over the boundary nodes that determine the colloidal shape is, in general, a good estimate of the electrostatic radius. However, for highly charged colloids, a more refined choice of the particle radius will be in general needed. The effective radius will be slightly larger than predicted from owing to the nonlinear decay of the electrostatic potential around the particle. This effective electrostatic radius, which needs to be calibrated as a function of the salt concentration and particle radius, will in general differ from the particle hydrodynamic radius and requires a separate analysis. In most situations, we expect that an equilibrium calibration correcting from the finite counter-ion concentration will provide a quantitative estimate for the electrokinetics of colloidal suspensions.

## Acknowledgements

G.G. acknowledges support from the IEF Marie Curie scheme. I.P. acknowledges financial support from Dirección General de Investigación (Spain) and DURSI under projects FIS 2008-04386 and 2009SGR-634, respectively.

## Footnotes

One contribution of 25 to a Theme Issue ‘Discrete simulation of fluid dynamics: applications’.

- This journal is © 2011 The Royal Society