Nonlinear electrophoresis of ideally polarizable particles

Figliuzzi, B.; Chan, W. H. R.; Moran, J. L.; Buie, C. R.

doi:10.1063/1.4897262

We focus in this paper on the nonlinear electrophoresis of ideally polarizable particles. At high applied voltages, significant ionic exchange occurs between the electric double layer, which surrounds the particle, and the bulk solution. In addition, steric effects due to the finite size of ions drastically modify the electric potential distribution in the electric double layer. In this situation, the velocity field, the electric potential, and the ionic concentration in the immediate vicinity of the particle are described by a complicated set of coupled nonlinear partial differential equations. In the general case, these equations must be solved numerically. In this study, we rely on a numerical approach to determine the electric potential, the ionic concentration, and the velocity field in the bulk solution surrounding the particle. The numerical simulations rely on a pseudo-spectral method which was used successfully by Chu and Bazant [J. Colloid Interface Sci. 315(1), 319–329 (2007)] to determine the electric potential and the ionic concentration around an ideally polarizable metallic sphere. Our numerical simulations also incorporate the steric model developed by Kilic et al [Phys. Rev. E 75, 021502 (2007)] to account for crowding effects in the electric double layer, advective transport, and for the presence of a body force in the bulk electrolyte. The simulations demonstrate that surface conduction significantly decreases the electrophoretic mobility of polarizable particles at high zeta potential and at high applied electric field. Advective transport in the electric double layer and in the bulk solution is also shown to significantly impact surface conduction.

I. INTRODUCTION

Electrophoresis is defined as the motion of colloidal particles under the influence of an applied electric field, and provides an efficient way to manipulate charged particles in microfluidic devices. Applications include lab-on-a-chip technologies^1,2 or DNA transport and separation.³ The mathematical description of electrophoresis of particles with thin electric double layer builds upon the classical theory developed by Smoluchowski⁴ in 1903. Smoluchowski noted that the net force exerted by an externally applied electric field on a system constituted by a suspended spherical particle and the counter-ion cloud around it is null since the total charge of the system vanishes. However, the co-existence of a charged region in the electric double layer and of externally applied field leads to the apparition of an electro-osmotic flow, which is in turn responsible for the motion of the particle in the electrolyte solution. Smoluchowski demonstrated that the velocity u_EP of a uniformly charged, dielectric particle of spherical shape relative to the electrolyte is proportional to the applied electric field E

\begin{equation}\mathbf {u}_{EP} = \mu _{EP}\mathbf {E},\end{equation}

u_{E P} = μ_{E P} E,

(1)

where μ_EP denotes the electrophoretic mobility of the particle. Later, Morrison⁵ and Anderson⁶ demonstrated that Smoluchowski's result is also valid for single non-spherical particles, and Sellier⁷ generalized it to assemblages of non-spherical particles with equal zeta potentials. Smoluchowski established his formula for particles whose zeta potential ζ remains small when compared to the thermal voltage φ_T. However, in 2012, Schnitzer and Yariv⁸ demonstrated that Smoluchowski's formula is valid even at moderate zeta potentials, ζ/φ_T ∼ O(1), and breaks down when ζ/φ_T ∼ O(ln R/λ_D).

With the development of microfluidics, large applied electric fields are commonly encountered in lab-on-a-chip technologies.^1,2 As a consequence, the determination of the mobility μ_EP at high voltages has been a central focus in colloid science for several decades.^9,10 Classicaly, dielectric particles have been considered, for which surface charges remain immobile. Bikerman¹¹ showed in 1940 that the surface conductivity of the electric double layer cannot be neglected in comparison with the conductivity of the bulk solution at large zeta potential. In this situation, a net exchange of ions occurs between the electric double layer and the bulk solution, which results in local ion depletion and accumulation in the immediate vicinity of the particle. Overbeek,¹² Dukhin and co-workers,^13–20 and O'Brien and co-workers^21–23 made significant progress toward the understanding of nonlinear electrophoresis of particles with fixed surface charges. Dukhin and co-workers^13–20 notably recognized that the electrophoretic mobility generally depends on the applied electric field and described the key role played by concentration polarization in electrokinetics phenomena at high voltages. In particular, they demonstrated that concentration polarization results in a fluid motion driven by the concentration gradient. This phenomenon, referred to as diffusio-osmosis, was independently rediscovered by Prieve et al. in 1984.²⁴ The combination of electro-osmosis and diffusio-osmosis leading to particle motion, referred to as diffusio-phoresis, has subsequently been investigated by Zaltzman and Rubinstein.^25,26 In 2010, Rica and Bazant²⁷ derived analytical formulas describing the electrophoretic mobility of colloids under direct current. When the electrolyte is subject to a DC electric field, concentration gradients arise in the bulk solution. As a consequence, the flow in the immediate vicinity of the colloidal particle has both electro- and diffusio-osmotic contributions. Rica and Bazant's analysis relies on the thin-double layer approximation and neglects both surface conduction and advective transport in the bulk solution. Diffusio-osmotic flow also plays a key role in the electrophoresis of cation-selective particles. For these particles, diffusio-osmotic flow can occur even in absence of surface conduction. These particles were recently studied by Khair,²⁸ who presented numerical simulations for AC and suddenly applied electric fields.

Most of the studies on nonlinear electrophoresis of particles with fixed surface charge rely on weak field approximations which enable a linearisation of the model equations. Recently, Schnitzer and Yariv have developed analytical models aimed at describing the strongly nonlinear regime^8,29,30 for particles with thin electric double layers. At asymptotically large applied voltage, they demonstrated that advective transport in the bulk electrolyte results in a uniform salt concentration. They also showed that the asymptotic matching between the current density emerging from the electric double layer and the current density in the electro-neutral bulk solution is incompatible with the asymptotic matching between respective salt fluxes. To resolve this apparent contradiction, they demonstrated that a diffuse boundary layer emerges in the overlap region of the electric double layer and bulk electrolyte.⁸

Nonlinear electrokinetics experienced a renewed interest in the mid 1990s, when various nonlinear electrokinetic phenomena were discovered and explored, including AC-electro-osmosis³¹ and induced-charge electrokinetic phenomena around conducting colloids.^32–34 Ideally polarizable particles are an archetypal problem in induced charge electrokinetics. These particles consist of an ideal conductor which enables a recombination of their surface charges when an electric field is applied. As a result, the zeta-potential varies on the surface of the particle and an induced electric field appears in the bulk solution. Electro-osmotic flows around a perfect conductor have been widely studied, notably by Squires and Bazant,^33,34 who suggested several applications in microfluidic devices, and Chu and Bazant,³⁵ who conducted numerical simulations to determine the electric potential and the ionic concentration around a conductor. Yariv³⁶ in 2008 and Bazant et al.³⁷ in 2009 have subsequently studied induced-charge electrophoresis. Relying on the Gouy-Chapman model to describe the electric double layer, Yariv obtained an asymptotical expression relating the electrophoretic mobility to the applied electric field,³⁶ which states that the electrophoretic mobility tends to zero at large applied electric fields. This result contradicts experimental evidence and the decay of the mobility with the applied electric field in Yariv's formula can be attributed to the inability of the Gouy-Chapman model to properly describe the physics of the electric double layer at large voltages.^38–40 More recently, Schnitzer and Yariv have derived analytical macroscale models accounting for surface conduction mechanisms in induced-charge electro-osmosis,⁴¹ in one of the first attempts to incorporate surface conduction in the description of induced-charge electrokinetic phenomena.

Classical approaches to determine the electrophoretic mobility of a charged particle rely on the determination of the electric potential, the ionic concentration, and the velocity field in the solution surrounding the particle. For large particles with thin electric double layer, the physical domain is classically divided into two regions, namely, the bulk and the electric double layer. Both regions are described by asymptotic solutions which are matched using appropriate boundary conditions. Several models have been developed to describe the electric double layer. The Gouy-Chapman model, developed independently in 1910 by Gouy⁴² and in 1913 by Chapman,⁴³ is usually the starting point of all attempts to describe nonlinear effects in the electric double layer. This model relies on dilute solution theory, where ions are described as non-interacting point charges. However, at high zeta potential, the Gouy-Chapman model predicts diverging ionic concentrations in the electric double layer. A first model aimed at correcting the non-physical divergence of the Gouy-Chapman description was proposed by Stern in 1924, who introduced the concept of compact layer.⁴⁴ As reviewed by Bazant et al.,³⁷ the first model accounting for steric effects due to the finite size of ions in the electric double layer was developed by Bikerman in 1942. Bikerman's model was rediscovered in 2007 by Kilic et al.,^39,40 who investigated the influence of steric effects due on the charge dynamics and conductivity of the electric double layer. In a 2009 paper,³⁷ Bazant et al. identified the importance of steric effects in induced-charge electrophoresis. Accounting for ionic steric effects, they found an asymptotic expression demonstrating that the electrophoretic mobility scales as the square root of the applied electric field at large voltages. This result fundamentally changes the prediction of Yariv,³⁶ which relies on dilute solution theory to describe the electric double layer at high voltage. The model of Kilic et al.^39,40 has been subsequently applied by Khair and Squires⁴⁵ to study ion steric effects on the electrophoresis of dieletric particles with uniform surface charge.

In this paper, we present a numerical model aimed at studying electrophoresis of ideally polarizable spherical particles at large applied electric fields. The novelty of this model is that it accounts for ionic steric effects, surface conduction, and advective transport in the bulk solution all together. The model equations are highly nonlinear and have to be handled numerically. The manuscript is organized as follows. In Sec. II we describe the mathematical model used to find the solution of the electrokinetics equations in the electric double layer. Specifically, following the model of Kilic et al.,^39,40 we rely on a modified Poisson-Boltzmann equation that takes into account steric effects due to the finite size of ions to obtain a modified potential distribution across the electric double layer. Additionally, we express the electrophoretic slip velocity at the bulk-electric double layer interface as a combination of electro-osmosis and diffusio-osmosis. In Sec. III, we describe the model used to find the corresponding solution in the bulk electrolyte, as well as the effective boundary conditions used to match both asymptotic solutions, discussing in detail the coupling between the electric potential, the ionic concentration, and the velocity fields in the solution around the particle. In particular, this treatment considers the contributions of diffusion, electromigration, and advection to the ionic fluxes in the bulk, and accounts for the presence of a body force. In Sec. IV, we discuss the implications of the model in the weakly nonlinear regime and present the results of numerical simulations in situations where crowding effects and surface conduction are significant. Conclusions are finally drawn in Sec. V.

II. ELECTRIC DOUBLE LAYER

A. Electric potential in the electric double layer

We rely on the model developed by Kilic et al.^39,40 to describe the electric potential distribution within the electric double layer (EDL). In this model, the electric field is determined by the local mean charge density, and the total free energy

$\mathcal {F}$

F

in the electric double layer can be expressed as the sum of the electrostatic energy

$\mathcal {U}$

U

and of an entropic contribution

$- T\mathcal {S}$

- T S

⁠. The electrostatic energy can be written as a function of the excess electric potential ϕ, defined as the electric potential difference ψ − ψ_b between the electric double layer and the bulk solution, and of the local ionic concentrations c₊ and c₋:

\begin{equation}\mathcal {U} = \int_V \bigg [{-}\frac{\epsilon _m}{2}|\nabla \phi |^2 + ze(c_+ - c_-)\phi \bigg] \mathbf {dr} .\end{equation}

U = \int_{V} [- \frac{ε_{m}}{2} {| \nabla ϕ |}^{2} + z e (c_{+} - c_{-}) ϕ] dr .

(2)

In this expression, z denotes the ion valence, ε_m the solvent permittivity, assumed to remain constant within the EDL, and e the elementary charge. In presence of crowding effects, the entropic contribution yields⁴⁰

\begin{eqnarray}\begin{array}{ccc} {-}T\mathcal {S} = \displaystyle\int _V \frac{kT}{\gamma }\bigg [\gamma c_+ \ln (\gamma c_+) + \gamma c_- \ln (\gamma c_-) \\[6pt]\quad\;+\, [1 - \gamma (c_+ - c_-)]\ln [1 - \gamma (c_+ + c_-)]\bigg ]\mathbf {dr}, \end{array}\end{eqnarray}

\begin{matrix} \begin{matrix} - T S = \int_{V} \frac{k T}{γ} [γ c_{+} \ln (γ c_{+}) + γ c_{-} \ln (γ c_{-}) \\ + [1 - γ (c_{+} - c_{-})] \ln [1 - γ (c_{+} + c_{-})]] dr, \end{matrix} \end{matrix}

(3)

where k denotes the Boltzmann constant and T ambient temperature. The dimensionless packing parameter γ introduced in the expression of

$- T\mathcal {S}$

- T S

accounts for ionic steric effects in the solution and can be related to the ionic diameter a and to the bulk solution concentration c_b,∞ far from the particle through the estimate³⁹

\begin{equation}2\gamma = 2a^3 c_{b, \infty } = \nu .\end{equation}

2 γ = 2 a^{3} c_{b, \infty} = ν .

(4)

Setting the variational derivatives

$\delta \mathcal {F}/ \delta \phi$

δ F / δ ϕ

and

$\delta \mathcal {F}/ \delta c_\pm$

δ F / δ c_{\pm}

equal to zero, we obtain the modified Poisson-Boltzmann equation

\begin{equation}\nabla ^2 \phi = \frac{2ze c_b}{\epsilon _m}\frac{\sinh \displaystyle\frac{ze\phi }{kT}}{1 + 2 \nu \sinh ^2 \displaystyle\frac{ze\phi }{2kT}}.\end{equation}

\nabla^{2} ϕ = \frac{2 z e c_{b}}{ε_{m}} \frac{\sinh \frac{z e ϕ}{k T}}{1 + 2 ν \sinh^{2} \frac{z e ϕ}{2 k T}} .

(5)

Here, c_b denotes the bulk concentration in the immediate vicinity of the electric double layer. As reviewed by Bazant,³⁷ Eq. (5) was derived for the first time by Bikerman³⁸ in 1942, and has since been independently rediscovered in several studies. For ν = 0, crowding effects are neglected and we recover the classical Poisson-Boltzmann equation. The present study focuses on large spherical particles whose radius R is orders of magnitude larger than the characteristic thickness of the EDL, which is classically given by the Debye length

\begin{equation}\lambda _D = \sqrt{\frac{\epsilon _m kT}{2z^2e^2 c_{b, \infty }} }.\end{equation}

λ_{D} = \sqrt{\frac{ε_{m} k T}{2 z^{2} e^{2} c_{b, \infty}}} .

(6)

For such particles, one can demonstrate through a rigorous matched asymptotic analysis⁴⁶ that the electric double layer attains a quasi-equilibrium in the direction normal to the surface, as sketched in Fig. 1. Gradients normal to the surface are indeed much stronger than gradients parallel to the surface. Therefore, the modified Poisson-Boltzmann equation becomes

\begin{equation}\frac{\partial ^2 \phi }{\partial y^2} = \frac{2ze c_b}{\epsilon _m}\frac{\sinh \displaystyle\frac{ze\phi }{kT}}{1 + 2 \nu \sinh ^2 \displaystyle\frac{ze\phi }{2kT}}.\end{equation}

\frac{\partial^{2} ϕ}{\partial y^{2}} = \frac{2 z e c_{b}}{ε_{m}} \frac{\sinh \frac{z e ϕ}{k T}}{1 + 2 ν \sinh^{2} \frac{z e ϕ}{2 k T}} .

(7)

At the surface of the particle, ϕ is equal to the zeta-potential ζ, defined as the potential drop between the particle surface and the bulk electrolyte.

FIG. 1.

View large Download slide

Schematic representation of the electrophoresis of an ideally polarizable particle of negative charge. In the thin double layer approximation, the calculation is performed by considering two asymptotic solutions in the bulk solution (left) and in the electric double layer (right).

At this point, we can employ dimensionless variables to facilitate the physical analysis. We define the Debye length λ_D as the characteristic length scale and the thermal voltage

\begin{equation}\varphi _T = \frac{kT}{ze}\end{equation}

φ_{T} = \frac{k T}{z e}

(8)

as the potential scale. With the aforementioned scales, the modified Poisson-Boltzmann equation reduces to

\begin{equation}\frac{\partial ^2 \phi }{\partial y^2} = \frac{ c_b \sinh \phi }{1 + 2 \nu \sinh ^2 \displaystyle\frac{\phi }{2}}.\end{equation}

\frac{\partial^{2} ϕ}{\partial y^{2}} = \frac{c_{b} \sinh ϕ}{1 + 2 ν \sinh^{2} \frac{ϕ}{2}} .

(9)

This equation can be integrated once to yield,³⁹ in dimensionless form,

\begin{equation}\frac{\partial \phi }{\partial y} = - \text{sgn} [\zeta ] \sqrt{ \frac{2 c_b}{\nu }\ln \bigg (1 + 2\nu \sinh ^2 \frac{\phi }{2} \bigg ) }.\end{equation}

\frac{\partial ϕ}{\partial y} = - sgn [ζ] \sqrt{\frac{2 c_{b}}{ν} \ln (1 + 2 ν \sinh^{2} \frac{ϕ}{2})} .

(10)

Equation (10) can in turn be integrated numerically to obtain the potential distribution across the EDL.

B. Electroosmotic flow in the electric double layer

Fluid motion in the electric double layer is the result of two distinct phenomena: first, the co-existence of a charged region in the EDL and of an externally applied electric field leads to electro-osmotic flow; second, the applied field polarizes the cloud of counter-ions surrounding the particle, which results in a fluid motion driven by the concentration gradient, a phenomenon usually referred to as diffusio-osmostic flow.^24,47

The low Reynolds number flow in the EDL is described by the Stokes equation. Within the lubrication approximation, the fluid motion across the EDL is asymptotically smaller than the flow along the EDL by a factor of O(λ_D/R), so that the Stokes equation in the e_y direction reads, at first order in O(λ_D/R),

\begin{equation}-\frac{\partial P}{\partial y} - \rho _E \frac{\partial \phi }{\partial y} = 0.\end{equation}

- \frac{\partial P}{\partial y} - ρ_{E} \frac{\partial ϕ}{\partial y} = 0 .

(11)

Combined with the Poisson-Boltzmann equation, relation (11) becomes

\begin{equation}\frac{\partial P}{\partial y} + 2zec_b \frac{\sinh \bigg ( \displaystyle\frac{ze \phi }{kT} \bigg )}{1 + 2 \nu \sinh ^2 \bigg ( \displaystyle\frac{ze \phi }{2kT} \bigg )} \frac{\partial \phi }{\partial y} = 0.\end{equation}

\frac{\partial P}{\partial y} + 2 z e c_{b} \frac{\sinh (\frac{z e ϕ}{k T})}{1 + 2 ν \sinh^{2} (\frac{z e ϕ}{2 k T})} \frac{\partial ϕ}{\partial y} = 0 .

(12)

Integration of (12) with respect to y yields the osmotic pressure profile across the EDL:

\begin{equation}P(y) = P_{\infty } + \frac{2kTc_b}{\nu }\ln \bigg [1 + 2 \nu \sinh ^2 \bigg ( \frac{ze \phi }{2kT} \bigg )\bigg ].\end{equation}

P (y) = P_{\infty} + \frac{2 k T c_{b}}{ν} \ln [1 + 2 ν \sinh^{2} (\frac{z e ϕ}{2 k T})] .

(13)

The projection of the Stokes equation on e_x reads

\begin{equation}-\frac{\partial P}{\partial x} - \rho _E \frac{\partial \psi }{\partial x} + \eta \frac{\partial ^2 u}{\partial y^2} = 0.\end{equation}

- \frac{\partial P}{\partial x} - ρ_{E} \frac{\partial ψ}{\partial x} + η \frac{\partial^{2} u}{\partial y^{2}} = 0 .

(14)

Using expression (13) for the pressure found previously, we obtain

\begin{eqnarray}\begin{split}\eta \frac{\partial ^2 u}{\partial y^2} = \frac{2kT}{\nu }\ln \bigg [1 + 2 \nu \sinh ^2 \bigg ( \frac{ze \phi }{2kT} \bigg )\bigg ] \frac{\partial c_b}{\partial x}&& \\-\, \epsilon _m \frac{\partial ^2 \phi }{\partial y^2} \frac{\partial \psi _b }{\partial x}.&&\end{split}\end{eqnarray}

\begin{matrix} \begin{matrix} η \frac{\partial^{2} u}{\partial y^{2}} = \frac{2 k T}{ν} \ln [1 + 2 ν \sinh^{2} (\frac{z e ϕ}{2 k T})] \frac{\partial c_{b}}{\partial x} \\ - ε_{m} \frac{\partial^{2} ϕ}{\partial y^{2}} \frac{\partial ψ_{b}}{\partial x} . \end{matrix} \end{matrix}

(15)

Interestingly, as noted by Prieve et al.,²⁴ the velocity only depends upon the electric field and the concentration in the bulk solution. To obtain the slip velocity of the particle, we can integrate relation (15) twice to obtain

\begin{equation}u_s = \frac{\epsilon _m \zeta }{\eta }\frac{\partial \psi _b }{\partial x} - \frac{k^2 T^2 \epsilon _m }{\eta z^2 e^2 } \kappa (\zeta ) \frac{1}{c_b}\frac{\partial c_b}{\partial x},\end{equation}

u_{s} = \frac{ε_{m} ζ}{η} \frac{\partial ψ_{b}}{\partial x} - \frac{k^{2} T^{2} ε_{m}}{η z^{2} e^{2}} κ (ζ) \frac{1}{c_{b}} \frac{\partial c_{b}}{\partial x},

(16)

where the dimensionless coefficient κ(ζ) is given by

\begin{equation}\kappa (\zeta ) = \int _0^\zeta \int _0^1 \frac{\phi ^{\prime }}{2} \sqrt{ \frac{\ln \bigg (1 + 2\nu \sinh ^2 \displaystyle\frac{\phi ^{\prime } \phi ^{\prime \prime }}{2}\bigg )}{\ln \bigg (1 + 2\nu \sinh ^2 \displaystyle\frac{\phi ^{\prime } }{2}\bigg )} } \mathrm{d}\phi ^{\prime } \mathrm{d}\phi ^{\prime \prime },\end{equation}

κ (ζ) = \int_{0}^{ζ} \int_{0}^{1} \frac{ϕ^{'}}{2} \sqrt{\frac{\ln (1 + 2 ν \sinh^{2} \frac{ϕ^{'} ϕ^{''}}{2})}{\ln (1 + 2 ν \sinh^{2} \frac{ϕ^{'}}{2})}} d ϕ^{'} d ϕ^{''},

(17)

and where ζ denotes the dimensionless zeta potential at the surface of the particle. Equation (16) describes the slip velocity as a combination of electro-osmotic and diffusio-osmotic slips, and is an equivalent form of the Dukhin-Deryaguin slip formula¹⁷ for “first kind” electro-osmosis at a thin quasi-equilibrium double layer, which accounts for steric effects due to the finite size of ions in the electric double layer. Rica and Bazant obtained an equivalent formula for the Gouy-Chapman model in their study of electrodiffusiophoresis.²⁷ A formulation for coefficient (17) was first obtained by Kilic⁴⁸ as a general function of osmotic pressure. In this study, coefficient (17) has been derived for the particular case of the steric model. Again, we can consider a dimensionless version of Eq. (16) by introducing the characteristic velocity

\begin{equation}U^* = \frac{\epsilon _m \varphi _T^2}{\eta R}.\end{equation}

U^{*} = \frac{ε_{m} φ_{T}^{2}}{η R} .

(18)

With this scaling, we find that the dimensionless slip velocity u_s at the surface of the particle is

\begin{equation}u_s = \zeta \frac{\partial \psi _b }{\partial x} - \frac{\kappa (\zeta )}{c_b}\frac{\partial c_b}{\partial x}.\end{equation}

u_{s} = ζ \frac{\partial ψ_{b}}{\partial x} - \frac{κ (ζ)}{c_{b}} \frac{\partial c_{b}}{\partial x} .

(19)

III. BULK EQUATIONS

A. Electric potential and ionic concentration

Following most prior work on electrophoresis of charged spherical particles with thin electric double layer, we assume that the bulk solution is electroneutral, so that

\begin{equation}c_\pm = c_b.\end{equation}

c_{\pm} = c_{b} .

(20)

We also assume the diffusivities of the cations and anions to be the same. The ion fluxes in the bulk solution are given by the Nernst-Planck relation

\begin{equation}\mathbf {j^\pm } = -D\nabla c_b \mp \frac{zeD}{kT}c_b \nabla \Psi _b + c_b \mathbf {u},\end{equation}

j^{\pm} = - D \nabla c_{b} \mp \frac{z e D}{k T} c_{b} \nabla Ψ_{b} + c_{b} u,

(21)

which describes the flux of each ionic species as a combination of diffusion, electromigration, and advection. It is convenient to nondimensionalize the governing equations of the model to facilitate its analysis. Scaling the ionic flux by the quantity

\begin{equation}J = \frac{Dc_{b, \infty }}{R},\end{equation}

J = \frac{D c_{b, \infty}}{R},

(22)

we obtain the dimensionless Nernst-Planck relation

\begin{equation}\mathbf {j^\pm } = -\nabla c_b \mp c_b \nabla \Psi _b + \text{m}_b c_b \mathbf {u},\end{equation}

j^{\pm} = - \nabla c_{b} \mp c_{b} \nabla Ψ_{b} + m_{b} c_{b} u,

(23)

where m_b is a dimensionless ionic drag coefficient defined for the bulk solution by

\begin{equation}\text{m}_b = \frac{U^*R}{D} = \frac{\epsilon_m \varphi _T^2}{\eta D}.\end{equation}

m_{b} = \frac{U^{*} R}{D} = \frac{ε_{m} φ_{T}^{2}}{η D} .

(24)

The conservation of counter-ions and co-ions are expressed by the relations

\begin{equation}\left\lbrace \begin{array}{@{}l@{\quad }l@{}}\text{div} (\mathbf {j^+ }) = -\nabla ^2 c_b - \text{div}(c_b \nabla \Psi _b) + \text{m}_b \text{div}(c_b \mathbf {u}) = 0,\\\\\text{div} (\mathbf {j^- }) = -\nabla ^2 c_b + \text{div}(c_b \nabla \Psi _b) + \text{m}_b \text{div}(c_b \mathbf {u}) = 0, \end{array}\right.\end{equation}

\{\begin{matrix} div (j^{+}) = - \nabla^{2} c_{b} - div (c_{b} \nabla Ψ_{b}) + m_{b} div (c_{b} u) = 0, \\ div (j^{-}) = - \nabla^{2} c_{b} + div (c_{b} \nabla Ψ_{b}) + m_{b} div (c_{b} u) = 0, \end{matrix})

(25)

which can be combined to yield, in dimensionless form,

\begin{equation}\left\lbrace \begin{array}{@{}l@{\quad }l@{}}\nabla ^2 c_b - \text{m}_b \text{div}(c_b \mathbf {u})= 0 ,\\\\\text{div}(c_b\nabla \Psi _b ) = 0. \end{array}\right.\end{equation}

\{\begin{matrix} \nabla^{2} c_{b} - m_{b} div (c_{b} u) = 0, \\ div (c_{b} \nabla Ψ_{b}) = 0 . \end{matrix})

(26)

These equations describe the concentration and the electric potential fields in the bulk solution and are classically referred to as the Poisson-Nernst-Planck equations.³⁵

The inner and the outer asymptotic solutions are matched by appropriate boundary conditions expressing the conservation of ions across the interface. Chu and Bazant⁴⁹ developed a systematic strategy to obtain effective boundary conditions at the diffuse interface between the electric double layer and the bulk solution, demonstrating that excess quantities have to be considered to establish the boundary conditions of conservative relations. Here, the excess ion fluxes in the electric double layer are defined as the difference between ion fluxes in the diffuse layer and ion fluxes in the bulk solution:

\begin{equation}\mathbf {J_{exc}^{\pm }} = -\frac{D}{kT}(c_{\pm } - c_b)\mathbf {\nabla } \mu _b + c_{\pm }\mathbf {u} - c_b \mathbf {u_s},\end{equation}

J_{exc}^{\pm} = - \frac{D}{k T} (c_{\pm} - c_{b}) \nabla μ_{b} + c_{\pm} u - c_{b} u_{s},

(27)

where μ_b denotes the electrochemical potential in the bulk solution, u the fluid velocity in the EDL, and u_s the slip velocity. Note that relation (27) describes ionic fluxes in the EDL as a combination of diffusion, electromigration, and advection. The excess surface fluxes can be calculated by integrating (27) between y = 0 and +∞ in the electric double layer, to yield, in outer spherical coordinates,

\begin{eqnarray}\begin{split}\mathbf {J_S^{\pm }} = \int _{0}^{+\infty } \mathbf {J_{exc}^{\pm }}(y)dy = -\frac{D}{kT}c_b \Gamma _{\pm } \mathbf {\nabla }_S \mu _b&& \\+\, c_b \Gamma _{\pm }\mathbf {u_s} + \int _{0}^{+\infty } c_{\pm }(\mathbf {u} - \mathbf {u_s})\mathrm{d}y.&& \end{split}\end{eqnarray}

\begin{matrix} \begin{matrix} J_{S}^{\pm} = \int_{0}^{+ \infty} J_{exc}^{\pm} (y) d y = - \frac{D}{k T} c_{b} Γ_{\pm} \nabla_{S} μ_{b} \\ + c_{b} Γ_{\pm} u_{s} + \int_{0}^{+ \infty} c_{\pm} (u - u_{s}) d y . \end{matrix} \end{matrix}

(28)

In this relation, Γ_± denotes the surface adsorption coefficient,³⁹ defined by

\begin{equation}\Gamma _{\pm } = \frac{1}{c_b}\int _{0}^{+\infty }(c_{\pm }(y) - c_b)\mathrm{d}y,\end{equation}

Γ_{\pm} = \frac{1}{c_{b}} \int_{0}^{+ \infty} (c_{\pm} (y) - c_{b}) d y,

(29)

and ∇_S denotes the surface gradient operator in the outer spherical coordinates system. The equation for ionic transport normal to the boundary is

\begin{equation}j^{\pm }_{n} = -\frac{D}{kT}c_b\frac{\partial \mu _b}{\partial r}.\end{equation}

j_{n}^{\pm} = - \frac{D}{k T} c_{b} \frac{\partial μ_{b}}{\partial r} .

(30)

Finally, the conservation of ionic species at the interface of the EDL reads

\begin{equation}j^{\pm }_{n} + \text{div}_S\mathbf {J_S^{\pm }}= 0,\end{equation}

j_{n}^{\pm} + {div}_{S} J_{S}^{\pm} = 0,

(31)

so that

\begin{equation}\left\lbrace \begin{array}{@{}l@{\quad }l@{}}j^{+}_{n} + j^{-}_{n} = -\text{div}_S \big(\mathbf {J_S^{+}} + \mathbf {J_S^{-}}\big),\\\\j^{+}_{n} - j^{-}_{n} = -\text{div} _S \big(\mathbf {J_S^{+}} - \mathbf {J_S^{-}}\big), \end{array}\right.\end{equation}

\{\begin{matrix} j_{n}^{+} + j_{n}^{-} = - {div}_{S} (J_{S}^{+} + J_{S}^{-}), \\ j_{n}^{+} - j_{n}^{-} = - {div}_{S} (J_{S}^{+} - J_{S}^{-}), \end{matrix})

(32)

where div_S denotes the surface divergence operator. In outer spherical coordinates, these relations become

\begin{equation}\left\lbrace \begin{array}{@{}l@{\quad }l@{}}\dfrac{\partial c_b}{\partial n} = -\dfrac{1}{2}\text{div}_S \bigg [w\nabla _S \ln c_b + q\dfrac{ze}{kT} \nabla _S \Psi _b - \dfrac{w\mathbf {u_s}}{D} - \\\dfrac{1}{D}\int _{0}^{+\infty } (c_+ + c_-)(\mathbf {u} - \mathbf {u_s})\mathrm{d}y \bigg ],\\\\\dfrac{zec_b}{kT}\dfrac{\partial \Psi _b}{\partial n} = -\dfrac{1}{2}\text{div}_S \bigg [q\nabla _S \ln c_b + w\dfrac{ze}{kT} \nabla _S \Psi _b \\- \dfrac{q\mathbf {u_s}}{D} - \dfrac{1}{D} \int _{0}^{+\infty } (c_+ - c_-)(\mathbf {u} - \mathbf {u_s})\mathrm{d}y \bigg ], \end{array}\right.\end{equation}

\{\begin{matrix} \frac{\partial c_{b}}{\partial n} = - \frac{1}{2} {div}_{S} [w \nabla_{S} \ln c_{b} + q \frac{z e}{k T} \nabla_{S} Ψ_{b} - \frac{w u_{s}}{D} - \\ \frac{1}{D} \int_{0}^{+ \infty} (c_{+} + c_{-}) (u - u_{s}) d y], \\ \frac{z e c_{b}}{k T} \frac{\partial Ψ_{b}}{\partial n} = - \frac{1}{2} {div}_{S} [q \nabla_{S} \ln c_{b} + w \frac{z e}{k T} \nabla_{S} Ψ_{b} \\ - \frac{q u_{s}}{D} - \frac{1}{D} \int_{0}^{+ \infty} (c_{+} - c_{-}) (u - u_{s}) d y], \end{matrix})

(33)

where q denotes the excess charge stored in the electric double layer and w the excess ion concentration, defined by

\begin{equation}q = \int _{0}^{+\infty } \bigg [c_+(y) - c_-(y)\bigg ]dy = \int _{0}^{+\infty } \rho (y)\mathrm{d}y\end{equation}

q = \int_{0}^{+ \infty} [c_{+} (y) - c_{-} (y)] d y = \int_{0}^{+ \infty} ρ (y) d y

(34)

and

\begin{equation}w = \int _{0}^{+\infty } \bigg [ c_+(y) + c_-(y)- 2c_b \bigg ]\mathrm{d}y.\end{equation}

w = \int_{0}^{+ \infty} [c_{+} (y) + c_{-} (y) - 2 c_{b}] d y .

(35)

The excess ion concentration can be related to the Dukhin-Bikerman number, Du, through the relation

\begin{equation}\text{Du} = \frac{w}{Rc_b}.\end{equation}

Du = \frac{w}{R c_{b}} .

(36)

With the steric model developed by Kilic and Bazant^39,40 and described in Sec. II, the excess charge q stored in the electric double layer is given as a function of the zeta potential by

\begin{equation}q = - 2zec_{b, \infty } \lambda _D\text{sgn}(\zeta ) \sqrt{\frac{2c_b}{\nu }\ln \bigg [ 1 + 2\nu \sinh ^2 \bigg (\frac{ze\zeta }{2kT} \bigg )\bigg ]}.\end{equation}

q = - 2 z e c_{b, \infty} λ_{D} sgn (ζ) \sqrt{\frac{2 c_{b}}{ν} \ln [1 + 2 ν \sinh^{2} (\frac{z e ζ}{2 k T})]} .

(37)

The excess ion concentration is defined by the integral expression^39,40

\begin{equation}w = \int _0 ^{\zeta } \frac{2 c_{b, \infty } \lambda _D \sqrt{c_b} (1 - \nu )}{\sqrt{\dfrac{2}{\nu }\ln ( 1 + 2\nu \sinh ^2(\frac{u}{2}))}} \frac{(\cosh u - 1)\mathrm{d}u}{1 + 2\nu \sinh ^2 (\frac{u}{2})}\end{equation}

w = \int_{0}^{ζ} \frac{2 c_{b, \infty} λ_{D} \sqrt{c_{b}} (1 - ν)}{\sqrt{\frac{2}{ν} \ln (1 + 2 ν \sinh^{2} (\frac{u}{2}))}} \frac{(\cosh u - 1) d u}{1 + 2 ν \sinh^{2} (\frac{u}{2})}

(38)

which can be evaluated numerically.

Scaling the zeta potential by the thermal voltage and the excess charge stored in the electric double layer by the quantity

\begin{equation}q^* = zec_{b, \infty } \lambda _D ,\end{equation}

q^{*} = z e c_{b, \infty} λ_{D},

(39)

we obtain the dimensionless expression

\begin{equation}q = -2\text{sgn}(\zeta )\sqrt{\frac{2c_b}{\nu }\ln \bigg [ 1 + 2\nu \sinh ^2(\zeta )\bigg ]} .\end{equation}

q = - 2 sgn (ζ) \sqrt{\frac{2 c_{b}}{ν} \ln [1 + 2 ν \sinh^{2} (ζ)]} .

(40)

Similarly, scaling the zeta potential by the thermal voltage and the excess ionic concentration in the electric double layer by the quantity

\begin{equation}w^* = c_{b, \infty } \lambda _D,\end{equation}

w^{*} = c_{b, \infty} λ_{D},

(41)

we find

\begin{equation}w = \int _0 ^{\zeta } 2 \sqrt{c_b} \frac{\cosh u - 1}{1 + 2\nu \sinh ^2 (\frac{u}{2})} \frac{(1 - \nu ) du}{\sqrt{\frac{2}{\nu }\ln ( 1 + 2\nu \sinh ^2(\frac{u}{2}))}}.\end{equation}

w = \int_{0}^{ζ} 2 \sqrt{c_{b}} \frac{\cosh u - 1}{1 + 2 ν \sinh^{2} (\frac{u}{2})} \frac{(1 - ν) d u}{\sqrt{\frac{2}{ν} \ln (1 + 2 ν \sinh^{2} (\frac{u}{2}))}} .

(42)

The dimensionless boundary conditions (33) become then

\begin{equation}\left\lbrace \begin{array}{@{}l@{\quad }l@{}}\dfrac{\partial c_b}{\partial n} = -\dfrac{\epsilon }{2}\text{div}_S \bigg [w\nabla _S \ln c_b + q\nabla _S \Psi _b - \\\text{m} \bigg ( w \mathbf {u_s} + \int _{0}^{+\infty } (c_+ + c_-)(\mathbf {u} - \mathbf {u_s})\mathrm{d}y \bigg ) \bigg ],\\\\c_b\dfrac{\partial \Psi _b}{\partial n} = -\dfrac{\epsilon }{2}\text{div}_S \bigg [q\nabla _S \ln c_b + w \nabla _S \Psi _b - \\\text{m}\bigg ( q\mathbf {u_s} + \int _{0}^{+\infty } (c_+ - c_-)(\mathbf {u} - \mathbf {u_s})\mathrm{d}y \bigg )\bigg ], \end{array}\right.\end{equation}

\{\begin{matrix} \frac{\partial c_{b}}{\partial n} = - \frac{ε}{2} {div}_{S} [w \nabla_{S} \ln c_{b} + q \nabla_{S} Ψ_{b} - \\ m (w u_{s} + \int_{0}^{+ \infty} (c_{+} + c_{-}) (u - u_{s}) d y)], \\ c_{b} \frac{\partial Ψ_{b}}{\partial n} = - \frac{ε}{2} {div}_{S} [q \nabla_{S} \ln c_{b} + w \nabla_{S} Ψ_{b} - \\ m (q u_{s} + \int_{0}^{+ \infty} (c_{+} - c_{-}) (u - u_{s}) d y)], \end{matrix})

(43)

where ε denotes the ratio between the Debye length and the radius of the particle and m is a dimensionless ionic drag coefficient defined for the EDL:

\begin{equation}\text{m} = \frac{\epsilon_m \varphi _T^2}{\eta D}.\end{equation}

m = \frac{ε_{m} φ_{T}^{2}}{η D} .

(44)

These boundary conditions state that fluxes of ions are transported across the electric double layer/bulk interface to balance the tangential surface flux gradients in the electric double layer.

B. Determination of the zeta potential

The zeta potential is defined as the electric potential difference across the electric double layer

\begin{equation}\zeta (\theta ) = \Psi _P - \psi _b (\theta ),\end{equation}

ζ (θ) = Ψ_{P} - ψ_{b} (θ),

(45)

where the potential Ψ_P of the particle is unknown. To evaluate the electric potential of the particle, as demonstrated by Yariv,³⁶ we have to rely on global charge conservation on the particle surface. The total charge Q of the particle remains indeed invariant during the formation of the electric double layer, so that

\begin{equation}Q = - \int _S qdA = - 2 \pi R^2 \int _0^{\pi } q(\theta ) \sin \theta \mathrm{d}\theta ,\end{equation}

Q = - \int_{S} q d A = - 2 π R^{2} \int_{0}^{π} q (θ) \sin θ d θ,

(46)

q being the local surface charge, related to the zeta potential through relation (37) when the steric model developed by Kilic et al.^39,40 is employed to describe the physics of the EDL. In dimensionless form, we obtain

\begin{equation}Q = \int _0^{\pi } \text{sign}(\zeta )\sqrt{\frac{2c_b}{\nu }\ln \bigg [ 1 + 2\nu \sinh ^2 \bigg (\frac{\zeta }{2} \bigg )\bigg ]} \sin \theta \mathrm{d}\theta ,\end{equation}

Q = \int_{0}^{π} sign (ζ) \sqrt{\frac{2 c_{b}}{ν} \ln [1 + 2 ν \sinh^{2} (\frac{ζ}{2})]} \sin θ d θ,

(47)

where the global charge Q of the particle has been scaled by the quantity

\begin{equation}Q^* = 4 \pi R^2 zec_{b, \infty } \lambda _D.\end{equation}

Q^{*} = 4 π R^{2} z e c_{b, \infty} λ_{D} .

(48)

C. Velocity field and electrophoretic velocity

As evidenced by expression (26), advective transport introduces a coupling between the ionic concentration and the velocity. Thus, the calculation of the electrophoretic mobility can only be achieved through a complete determination of the velocity field. In the bulk solution, the velocity is solution of the Stokes equation

\begin{equation}\eta \nabla ^2 \mathbf {u} - \nabla P + \epsilon _m \nabla ^2 \Psi _b \nabla \Psi _b = \mathbf {0},\end{equation}

η \nabla^{2} u - \nabla P + ε_{m} \nabla^{2} Ψ_{b} \nabla Ψ_{b} = 0,

(49)

where u denotes the fluid velocity, P the pressure field, ε_m the solvent permeability, and Ψ_b the electric potential in the bulk electrolyte. In this equation, the term ε_m∇²Ψ_b∇Ψ_b refers to the effects of the body force. The Stokes equation is considered along with the continuity equation

\begin{equation}\nabla \cdot \mathbf {u} = 0.\end{equation}

\nabla \cdot u = 0 .

(50)

The problem considered here is very similar to the one of a sphere in a Stokes flow, with the notable exceptions that slip boundary conditions have to be considered at the surface of the sphere and that an external force appears in the expression of the Stokes equation. The field-induced variations of the bulk concentration indeed result in corresponding variations of the solution conductivity, which in turn result in a body force acting on the bulk fluid.

The classical resolution of the problem relies on the vorticity, defined by ω = ∇ × u. In spherical coordinates, for the bi-dimensional flow considered here, the vorticity is directed along the basis vector e_ψ. Using the continuity equation, we find that

\begin{equation}\nabla \times \mathbf {\omega } = \nabla \times \nabla \times \mathbf {u } = - \nabla ^2 \mathbf {u }.\end{equation}

\nabla \times ω = \nabla \times \nabla \times u = - \nabla^{2} u .

(51)

The momentum conservation can then be reformulated in terms of the vorticity to be

\begin{equation}- \eta \nabla \times \mathbf {\omega } - \nabla P + \epsilon _m \nabla ^2 \Psi _b \nabla \Psi _b = \mathbf {0}.\end{equation}

- η \nabla \times ω - \nabla P + ε_{m} \nabla^{2} Ψ_{b} \nabla Ψ_{b} = 0 .

(52)

Taking the curl of this relation, we finally obtain

\begin{equation}\eta \nabla ^2 \mathbf {\omega } + \nabla \times \epsilon _m \nabla ^2 \Psi _b \nabla \Psi _b = \mathbf {0}.\end{equation}

η \nabla^{2} ω + \nabla \times ε_{m} \nabla^{2} Ψ_{b} \nabla Ψ_{b} = 0 .

(53)

At this point, we introduce the Stokes stream function Ψ defined in spherical coordinates by

\begin{equation}u_r = \frac{1}{r^2 \sin \theta }\frac{\partial \Psi }{\partial \theta }, \quad u_\theta = -\frac{1}{r \sin \theta }\frac{\partial \Psi }{\partial r}.\end{equation}

u_{r} = \frac{1}{r^{2} \sin θ} \frac{\partial Ψ}{\partial θ}, u_{θ} = - \frac{1}{r \sin θ} \frac{\partial Ψ}{\partial r} .

(54)

The vorticity can be expressed as a function of the Stokes stream function through the relation

\begin{equation}\mathbf {\omega } = -\frac{1}{r \sin \theta } \mathcal {L} \Psi ,\end{equation}

ω = - \frac{1}{r \sin θ} L Ψ,

(55)

where operator

$\mathcal {L}$

L

is defined by

\begin{equation}\mathcal {L} = \frac{\partial ^2 .}{\partial r ^2} + \frac{1}{r^2}\frac{\partial ^2 .}{\partial \theta ^2} - \frac{\cot \theta }{r^2}\frac{\partial .}{\partial \theta }.\end{equation}

L = \frac{\partial^{2} .}{\partial r^{2}} + \frac{1}{r^{2}} \frac{\partial^{2} .}{\partial θ^{2}} - \frac{\cot θ}{r^{2}} \frac{\partial .}{\partial θ} .

(56)

We can finally reformulate the momentum conservation expressed by Eq. (52) in terms of the Stokes stream function to yield

\begin{equation}\eta \nabla ^2 \bigg (\frac{1}{r \sin \theta } \mathcal {L} \bigg )\Psi - \nabla \times \epsilon _m \nabla ^2 \Psi _b \nabla \Psi _b = \mathbf {0}.\end{equation}

η \nabla^{2} (\frac{1}{r \sin θ}, L) Ψ - \nabla \times ε_{m} \nabla^{2} Ψ_{b} \nabla Ψ_{b} = 0 .

(57)

We complete the problem description by specifying boundary conditions in the reference frame anchored at the center of the particle. At infinity, the flow is a steady uniform stream at velocity −U, where U is the electrophoretic velocity of the particle, which remains unknown. At the surface of the particle, the velocity field must match the slip velocity given by relation (16):

\begin{equation}\frac{1}{r^2 \sin \theta }\frac{\partial \Psi }{\partial \theta }\bigg |_{r = R} = 0 , \quad -\frac{1}{r \sin \theta }\frac{\partial \Psi }{\partial r}\bigg |_{r = R} = u_s.\end{equation}

\frac{1}{r^{2} \sin θ} \frac{\partial Ψ}{\partial θ} |_{r = R} = 0, - \frac{1}{r \sin θ} \frac{\partial Ψ}{\partial r} |_{r = R} = u_{s} .

(58)

D. Electrophoretic velocity

To close the model, we still have to determine U. Since the global charge of the particle and its surrounding EDL is zero, there is no net force acting on the system. As shown by Stone and Samuel,⁵⁰ for a spherical particle with an arbitrary slip distribution u_s(r) on its surface, the translational velocity is equal and opposite to the averaged slip velocity on the surface:

\begin{equation}\mathbf {U} = - \frac{1}{4 \pi R^2}\int _S \mathbf {u_s} \mathrm{d}S.\end{equation}

U = - \frac{1}{4 π R^{2}} \int_{S} u_{s} d S .

(59)

A similar expression was obtained by Anderson⁶ in 1989 for the specific case of electrophoresis. However, in our case, due to the existence of a body force, Eq. (59) does not apply.

We denote by σ_E the Maxwell stress tensor and by σ_H the hydrodynamic stress tensor. The mechanical equilibrium of the system constituted by the particle and its surrounding EDL reads

\begin{equation}\int _{S} \mathbf {\sigma _E} \cdot \mathbf {n}\mathrm{d}S + \int _{S} \mathbf {\sigma _H} \cdot \mathbf {n}\mathrm{d}S= 0,\end{equation}

\int_{S} σ_{E} \cdot n d S + \int_{S} σ_{H} \cdot n d S = 0,

(60)

where S denotes the external surface of the EDL. In this expression, the Maxwell stress tensor is given by

\begin{equation}\mathbf {\sigma _E} = \epsilon _m [ \nabla \psi _b \nabla \psi _b - (\nabla \psi _b \cdot \nabla \psi _b) \mathrm{I}].\end{equation}

σ_{E} = ε_{m} [\nabla ψ_{b} \nabla ψ_{b} - (\nabla ψ_{b} \cdot \nabla ψ_{b}) I] .

(61)

To facilitate the determination of the hydrodynamic force, we can make use of the reciprocal theorem.⁵¹ To that end, let

$\mathbf {\tilde{u}}$

\tilde{u}

and

$\mathbf {\tilde{\sigma }_H}$

{\tilde{σ}}_{H}

be the velocity field and hydrodynamic stress tensor, respectively, belonging to a companion Stokes flow about a spherical object in translation at a unit velocity in direction e_z in an unbounded fluid at rest at infinity. For the companion flow, the body force vanishes. The velocity field

$\mathbf {\tilde{u}}$

\tilde{u}

is classically given by

\begin{equation}\mathbf {\tilde{u}} = \frac{1}{2} \cos \theta \bigg ( \frac{3R}{r} - \frac{R^3}{r^3}\bigg )\mathbf {e_r} - \frac{1}{4} \sin \theta \bigg ( \frac{3R}{r} + \frac{R^3}{r^3}\bigg )\mathbf {e_\theta }.\end{equation}

\tilde{u} = \frac{1}{2} \cos θ (\frac{3 R}{r} - \frac{R^{3}}{r^{3}}) e_{r} - \frac{1}{4} \sin θ (\frac{3 R}{r} + \frac{R^{3}}{r^{3}}) e_{θ} .

(62)

The reciprocal theorem states that

\begin{eqnarray}\begin{split}\int _{S} \tilde{u} \cdot \mathbf {\sigma _H} \cdot \mathbf {n}\mathrm{d}S + \int _{S_{\infty }} \tilde{u} \cdot \mathbf {\sigma _H} \cdot \mathbf {n}\mathrm{d}S + \int _{V} \mathbf {\tilde{u}} \cdot \mathbf {f}\mathrm{d}V&&\\= \int _{S} \mathbf {u} \cdot \mathbf {\tilde{\sigma }_H} \cdot \mathbf {n}\mathrm{d}S + \int _{S_{\infty }} \mathbf {u} \cdot \mathbf {\tilde{\sigma }_H} \cdot \mathbf {n}\mathrm{d}S,&& \end{split}\end{eqnarray}

\begin{matrix} \begin{matrix} \int_{S} \tilde{u} \cdot σ_{H} \cdot n d S + \int_{S_{\infty}} \tilde{u} \cdot σ_{H} \cdot n d S + \int_{V} \tilde{u} \cdot f d V \\ = \int_{S} u \cdot {\tilde{σ}}_{H} \cdot n d S + \int_{S_{\infty}} u \cdot {\tilde{σ}}_{H} \cdot n d S, \end{matrix} \end{matrix}

(63)

where S_∞ denotes a surface at infinity (r → ∞). On the left-hand side of (63), the surface integral at infinity yields −6πηRU when projected on e_z. The hydrodynamic stress tensor of the companion flow is uniform over the surface of the particle.⁵¹ Thus, at the surface of the particle, the traction yields

\begin{equation}\mathbf {\tilde{\sigma }_H} \cdot \mathbf {n}|_{r = R} = -\frac{3\eta }{2R}\mathbf {e_z}\end{equation}

{\tilde{σ}}_{H} {\cdot n |}_{r = R} = - \frac{3 η}{2 R} e_{z}

(64)

so that

\begin{equation}\int _{S} \mathbf {u} \cdot \mathbf {\tilde{\sigma }_H} \cdot \mathbf {n}\mathrm{d}S = -\frac{3\eta }{2R}\mathbf {e_z} \cdot \bigg ( \int _S \mathbf {u_s} \mathrm{d}S \bigg ).\end{equation}

\int_{S} u \cdot {\tilde{σ}}_{H} \cdot n d S = - \frac{3 η}{2 R} e_{z} \cdot (\int_{S}, u_{s}, d, S) .

(65)

Similarly, the velocity of the companion flow is uniform over the surface of the sphere. Therefore, we have

\begin{equation}\int _{S} \mathbf {\tilde{u} } \cdot \mathbf {\sigma _H} \cdot \mathbf {n}\mathrm{d}S = U\mathbf {e_z} \cdot \int _{S} \mathbf {\sigma _H} \cdot \mathbf {n}\mathrm{d}S.\end{equation}

\int_{S} \tilde{u} \cdot σ_{H} \cdot n d S = U e_{z} \cdot \int_{S} σ_{H} \cdot n d S .

(66)

Finally, the integrand for the surface integral at infinity scales as 1/r when r tends to infinity, so that the integral vanishes. As a consequence, the reciprocity theorem can be rearranged to yield

\begin{eqnarray}\begin{split}\mathbf {e_z} \cdot \int _{S} \mathbf {\sigma _H} \cdot \mathbf {n}\mathrm{d}S = - 6 \pi \eta U R -\frac{3\eta }{2R}\mathbf {e_z} \cdot \bigg ( \int _S \mathbf {u_s} \mathrm{d}S \bigg )&&\\-\, \int _{V} \mathbf {\tilde{u}} \cdot \mathbf {f}\mathrm{d}V . &&\end{split}\end{eqnarray}

\begin{matrix} \begin{matrix} e_{z} \cdot \int_{S} σ_{H} \cdot n d S = - 6 π η U R - \frac{3 η}{2 R} e_{z} \cdot (\int_{S}, u_{s}, d, S) \\ - \int_{V} \tilde{u} \cdot f d V . \end{matrix} \end{matrix}

(67)

Using this expression along with (60), we can compute the electrophoretic velocity U of the particle:

\begin{eqnarray}\begin{split} U= \frac{1}{6\pi \eta R}\bigg ( \mathbf {e_z} \cdot \int _{S} \mathbf {\sigma _E} \cdot \mathbf {n}\mathrm{d}S -\frac{3\eta }{2R}\mathbf {e_z} \cdot \int _S \mathbf {u_s} \mathrm{d}S &&\\-\, \mathbf {e_z} \cdot \int _{V} \mathbf {\tilde{u}} \cdot \mathbf {f}\mathrm{d}V \bigg ). &&\end{split}\end{eqnarray}

\begin{matrix} \begin{matrix} U = \frac{1}{6 π η R} (e_{z} \cdot \int_{S} σ_{E} \cdot n d S - \frac{3 η}{2 R} e_{z} \cdot \int_{S} u_{s} d S \\ - e_{z} \cdot \int_{V} \tilde{u} \cdot f d V) . \end{matrix} \end{matrix}

(68)

IV. RESULTS AND DISCUSSION

A. Numerical model

The complete determination of the electrophoretic mobility requires the calculation of the ionic concentration c_b, of the electric potential ψ_b and of the velocity u. Thus, nonlinear electrophoresis of ideally polarizable particles with thin electric double layer is described by a rather complicated set of coupled partial differential equations.

In the bulk solution, the conservation of each ionic species is described by the advection-diffusion equations (26), while the Stokes stream function is solution of (57). Effective boundary conditions are introduced to match the asymptotic solutions found in the bulk solution and in the electric double layer, which rely on the conservation of ions through the diffuse electric double layer (Eq. (43)) and on the no-slip boundary condition at the interface between the bulk and the electric double layer (Eq. (58)). The velocity at infinity is obtained by considering the conservation of momentum (Eq. (59)) for the charged particle along with its surrounding electric double layer. Finally, the electric potential of the particle is obtained by employing the global charge conservation relation (Eq. (47)). For particular situations, for example when surface conduction is negligible in the electric double layer as it is the case for low applied electric field, the above equations can be simplified and handled analytically. However, in general, these equations must be solved numerically.

In this study, we used a pseudo-spectral method to solve the model equations. Chu and Bazant³⁵ relied on the same numerical method to perform the numerical calculation of the ionic concentration and of the electric potential fields around a conductor. The method is described in detail in Boyd.⁵² Following their study, we employed a tensor product of a uniform spaced grid for the polar angle direction and a semi-infinite rational Chebyshev grid for the radial direction to define the computational grid. We took advantage of the axisymmetry of the problem to reduce all the model equations to two dimensions. Finally, we used standard Newton iterations to solve the nonlinear equations involved in the mathematical description of the model. All numerical simulations were conducted using the software Octave. Note that the integral term describing advective transport in the EDL in Eq. (43) is computed numerically using the routine quad implemented in Octave. All simulations were realized on a grid of 80 points (radial direction) by 70 points (polar angle direction).

The mathematical model is governed by five dimensionless physical parameters. The first one is the dimensionless charge of the particle, defined by relation (48). Using the Debye description of the electric double layer, we can relate the dimensionless particle charge to the uniform equilibrium zeta potential of the particle:

\begin{equation}Q \simeq 2\zeta _0 .\end{equation}

Q ≃ 2 ζ_{0} .

(69)

In this relation, ζ₀ is scaled by the thermal voltage (25 mV for a univalent solution). The second parameter is the dimensionless applied electric field E. The characteristic electric field is defined here as

\begin{equation}E^* = \frac{\varphi _T}{R},\end{equation}

E^{*} = \frac{φ_{T}}{R},

(70)

and consequently depends on the particle size. For a 1 μm particle, the characteristic electric field is E* = 250 V/cm. The third parameter is the ratio ε between the Debye length λ_D and the particle radius R. The Debye length is essentially fixed by the bulk concentration c_{b, ∞}. For a highly concentrated electrolyte solution at 25 °C where c_{b, ∞} = 0.18 M, the Debye length is λ_D = 0.72 nm. Thus, for a 1 μm particle, we find ε ≃ 1.0 × 10⁻³. For an electrolyte solution at 25 °C where c_{b, ∞} = 1.0 × 10⁻³ M, the Debye length is λ_D = 9.6 nm, which results in ε ≃ 1.0 × 10⁻² for a 1 μm particle. The fourth parameters are the dimensionless ionic drag coefficients m and m_b given by (44), which arise in the equations describing advective-diffusive transport in the EDL and in the bulk solution, respectively. For an ionic diffusivity D = 10⁻⁹ m²/s, m and m_b can be estimated to be around 0.45. Finally, the last parameter is the packing parameter ν describing crowding effects, which is defined by relation (4) and depends on the ionic concentration in the solution. ν varies between 1.0 × 10⁻⁵ for diluted solutions and 0.1 for highly concentrated solutions.

B. Weakly nonlinear regime

Before discussing the numerical simulations, we briefly focus on the case of weak applied electric fields to compare this study to the analytical expressions found by Yariv³⁶ for the electrophoretic velocity of ideally polarizable particles. At low applied electric field, the excess charge q and the excess ionic concentration w remain of order O(1), so that the right-hand sides of Eq. (43) remain of order O(ε). It is consequently legitimate to consider the limit ε → 0, where surface conduction effects can be neglected. According to Eq. (26), the concentration is uniform in the bulk solution. The electric potential is then a harmonic function and can be uniquely determined with the boundary conditions³⁶

\begin{equation}\Psi _b(r, \theta ) = -E\bigg (r + \frac{R^3}{2r^2}\bigg ) \cos \theta .\end{equation}

Ψ_{b} (r, θ) = - E (r + \frac{R^{3}}{2 r^{2}}) \cos θ .

(71)

As a consequence, the zeta potential is

\begin{equation}\zeta (\theta ) = \Psi _P + \frac{3}{2}E\cos \theta .\end{equation}

ζ (θ) = Ψ_{P} + \frac{3}{2} E \cos θ .

(72)

When crowding effects are neglected, we recover the Gouy-Chapman model from charge conservation as shown in Eq. (47)

\begin{equation}Q = \int _0^{\pi } 2 \sinh \frac{\zeta }{2} \sin \theta \mathrm{d}\theta .\end{equation}

Q = \int_{0}^{π} 2 \sinh \frac{ζ}{2} \sin θ d θ .

(73)

Equation (73) can be used in conjunction with Eq. (72) to obtain an expression for the particle potential Ψ_P, which can be correlated with the charge Q of the particle, in a similar form to that obtained by Yariv (2008)³⁶

\begin{equation}\Psi _P = 2\sinh ^{-1} \left( \frac{3QE}{16 \sinh \frac{3E}{4}} \right) .\end{equation}

Ψ_{P} = 2 \sinh^{- 1} (\frac{3 Q E}{16 \sinh \frac{3 E}{4}}) .

(74)

For the particular form (71) of the electric potential, Ψ_P is simply proportional to the electrophoretic mobility of the particle. The numerical simulations performed with our model demonstrate convergence towards the results of Yariv when ν = 0 for different values of Q as evidenced in Figs. 2 and 3.

FIG. 2.

View large Download slide

Dimensionless electrophoretic mobility (a) and velocity (b) of a positively charged particle as a function of the dimensionless applied electric field for distinct values of the packing parameter ν. The dimensionless charge of the particle is Q = 1. Surface conduction is neglected in these simulations (ε = 0). The asymptotic approximation (80) is plotted in dashed line for E > 6. The analytical model (74) of Yariv is plotted in dashed line and perfectly match the curve obtained from numerical simulations for ν = 0.

FIG. 3.

View large Download slide

Dimensionless electrophoretic mobility (left) and velocity (right) of a positively charged particle as a function of the dimensionless applied electric field for distinct values of the packing parameter ν. The dimensionless charge of the particle is Q = 5. Surface conduction is neglected in these simulations (ε = 0). The asymptotic approximation (80) is plotted in dashed line for E > 6. The analytical model (74) of Yariv is plotted in dashed line and perfectly match the curve obtained from numerical simulations for ν = 0. In this simulation, the influence of steric effects can be observed even at small applied electric field due to the higher particle charge.

C. Crowding effects

Figs. 2 and 3 demonstrate that the mobility of a charged polarizable spherical particle is highly dependent on the packing parameter ν. This trend is markedly significant at large values of E.

At low values of E and Q, steric effects are relatively insignificant. Hence, the mobility values are similar for all values of ν considered, and are well predicted by the model developed by Yariv³⁶ describing the weakly nonlinear regime. This similarity breaks down for increasing Q as shown in Figure 3, even at low applied electric fields. Physically, since the zeta potential increases at higher values of Q, the extent to which the solvated ions are packed around the particle influences its distribution around the particle, as shown in Fig. 4.

FIG. 4.

View large Download slide

Dimensionless surface charge for distinct values of the packing parameter ν. The dimensionless applied electric field is E = 4 (a) and E = 7 (b). In these simulations, Q = 1 and ε = 0. The particle polarization can be observed in both figures. The surface charge diverges with increasing applied electric field when steric effects are neglected (ν = 0).

At high values of E, ζ increases significantly near the poles of the particle, and steric effects influence the electrophoretic mobility. Starting from global charge conservation (47), we can derive an asymptotic expression for the mobility in absence of surface conduction (ε → 0). In this case, the bulk concentration remains constant, and the zeta potential is given by

\begin{equation}\zeta (\theta ) = \frac{3}{2}E\cos \theta + \Psi _P.\end{equation}

ζ (θ) = \frac{3}{2} E \cos θ + Ψ_{P} .

(75)

As a consequence, the integral relation (47) expressing global charge conservation on the particle can be reformulated as

\begin{equation}Q = \frac{2}{3E}\int _{\zeta (\pi )}^{\zeta (0)} \text{sgn} (\zeta ) \sqrt{\frac{2}{\nu } \ln \bigg [1 + 2\nu \sinh ^2 \bigg (\frac{\zeta }{2} \bigg )}\bigg ] \mathrm{d}\zeta .\end{equation}

Q = \frac{2}{3 E} \int_{ζ (π)}^{ζ (0)} sgn (ζ) \sqrt{\frac{2}{ν} \ln [1 + 2 ν \sinh^{2} (\frac{ζ}{2})}] d ζ .

(76)

The integrand being uneven with respect to ζ, we find, in the case of a positively charged particle,

\begin{equation}\frac{3EQ}{2} = \int _{\frac{3}{2}E - \Psi _P }^{\frac{3}{2}E + \Psi _P } \sqrt{\frac{2 }{\nu } \ln \bigg [1 + 2\nu \sinh ^2 \bigg (\frac{\zeta }{2} \bigg )}\bigg ] \mathrm{d}\zeta .\end{equation}

\frac{3 E Q}{2} = \int_{\frac{3}{2} E - Ψ_{P}}^{\frac{3}{2} E + Ψ_{P}} \sqrt{\frac{2}{ν} \ln [1 + 2 ν \sinh^{2} (\frac{ζ}{2})}] d ζ .

(77)

At this point, we assume that the electric field is high enough such that

\begin{equation}\frac{3}{2}E - \Psi _P \gg 1 .\end{equation}

\frac{3}{2} E - Ψ_{P} ≫ 1 .

(78)

In this case, the argument in the logarithm in Eq. (77) is large, and the numerical value of the integrand can be approximated as

\begin{equation}\sqrt{\frac{2}{\nu } \ln \bigg [1 + 2\nu \sinh ^2 \bigg (\frac{\zeta }{2} \bigg )\bigg ]} \simeq \sqrt{\frac{2}{\nu }\bigg (\ln {\frac{\nu }{2}} + \zeta \bigg )}\end{equation}

\sqrt{\frac{2}{ν} \ln [1 + 2 ν \sinh^{2} (\frac{ζ}{2})]} ≃ \sqrt{\frac{2}{ν} (\ln \frac{ν}{2} + ζ)}

(79)

for reasonably large ν, to yield

\begin{equation}\Psi _P \simeq \frac{3QE}{4} \sqrt{\frac{\nu }{2\ln {\frac{\nu }{2}}+3E}} .\end{equation}

Ψ_{P} ≃ \frac{3 Q E}{4} \sqrt{\frac{ν}{2 \ln \frac{ν}{2} + 3 E}} .

(80)

At very high values of E, the limit of Eq. (80) finally becomes

\begin{equation}\Psi _P \simeq \frac{Q \sqrt{3\nu E}}{4} \propto \sqrt{E}.\end{equation}

Ψ_{P} ≃ \frac{Q \sqrt{3 ν E}}{4} \propto \sqrt{E} .

(81)

This asymptotic formula was first found and interpreted by Bazant et al.³⁷ in 2009. A physical account of relation (81) is that without considering the packing of solvated ions around the particle, the surface charge of the particle diverges with increasing electric fields. It follows from this that, as shown in Fig. 4, the zeta potential distribution around the particle is relatively symmetric, which leads to a low electrophoretic mobility. On the other hand, for large values of ν, the distribution of solvated ions around the particle is affected by steric effects due to the finite size of ions in the electric double layer, which preserves the asymmetry of the zeta potential distribution. The numerical simulations show good agreement with the asymptotic expression (80), as evidenced in Figs 2 and 3.

D. Influence of surface conduction

The numerical simulations also suggest that electrophoretic mobility decreases with increased surface conduction, as shown in Fig. 5 for the case of particles with increasing global charge. For the simulations in Fig. 5, the dimensionless electric field was held fixed at a value E = 0.5. At low values of Q and hence ζ, the excess charge q and excess ion concentration w in the electric double layer remain small. Hence, according to Eq. (43) describing ionic species conservation across the bulk/EDL interface, surface conduction is relatively insignificant. Since the particle does not significantly distort the concentration field in its vicinity, concentration polarization remains negligible, as shown in Fig. 6. Accordingly, Fig. 7 shows that the contribution of the diffusio-osmotic flow to the slip velocity is practically zero when the charge of the studied particle remains small.

FIG. 5.

View large Download slide

Dimensionless electrophoretic mobility of a positively charged particle as a function of the particle charge for distinct values of the ratio ε of the Debye length to the radius of the particle. We observe that electrophoretic mobility decreases with increased surface conduction, as shown by Khair and Squires⁴⁵ for the case of fixed surface charge particles.

FIG. 6.

View large Download slide

Concentration polarization at the immediate vicinity of the charged particle for ε = 0.01 (a) and ε = 0.02 (b). The simulation parameters are here E = 0.5, ν = 0.1.

FIG. 7.

View large Download slide

Diffusio-osmotic flow (a) and slip velocity (b) at the surface of a charged particle for ε = 0.02. The simulation parameters are here E = 0.5, ν = 0.1. For highly charged particles, significant diffusio-osmotic flow can be observed (a) which reduces the total slip at the surface of the particle (b).

At higher values of Q, the excess charge q and the excess ionic concentration w increase, so that surface conduction becomes significant in the EDL. As a consequence, the concentration field becomes polarized in the immediate vicinity of the particle. This physical effect is shown in Fig. 6, and is enhanced at high values of the parameter ε which relates the thickness of the EDL to the radius of the particle. As shown in Fig. 7, concentration polarization results in significant diffusio-osmotic flow, which in turns reduces the slip velocity at the surface of the particle. Accordingly, the net electrophoretic motion of the particle is significantly reduced.

A similar analysis can be conducted when larger electric fields are applied, as shown in Fig. 8. In this case, the zeta potential reaches high values on the surface of the particle. This results in significant surface conduction, as evidenced by the values of the Dukhin number plotted in Fig. 9, which in turn enhances concentration polarization and diffusio-osmotic flow. However, in this case, the diffusio-osmotic flow contribution to the slip velocity remains much smaller than the contribution of electro-osmotic flow, as shown in Fig. 9. Additionally, according to Eq. (26), concentration gradients enhance advective transport in the bulk solution, which significantly reduces ionic concentration fluctuations in the electric double layer. This result is in accordance with the asymptotic analysis conducted by Schnitzer and Yariv³⁰ for dielectric particles, which demonstrated that the concentration remains constant at first order in the bulk solution at high Peclet number. The ionic concentration field around the charged particle is shown in Fig. 10 when advective transport is neglected and when the dimensionless ionic drag coefficient takes a typical value m_b = 0.45.

FIG. 8.

View large Download slide

Dimensionless electrophoretic mobility of a positively charged particle as a function of the dimensionless applied electric field for distinct values of the ratio ε. The dimensionless ionic drag coefficients m_b and m are alternatively set to 0 and 0.45 to show the effects of advective transport in the bulk solution and in the EDL, respectively. The packing parameter ν = 0.015 corresponds to a bulk concentration c_{b, ∞} = 0.1M for an ionic diameter a = 0.5 nm. The Debye length is thus λ_D = 0.96 nm. For ε = 0.01, the particle radius is R = 96 nm, so that the characteristic electric field is E* = 2.5 × 10³ V/cm.

FIG. 9.

View large Download slide

Diffusio-osmotic flow (a) and slip velocity (b) at the surface of a moderately charged particle for distinct values of the parameter ε. The dimensionless ionic drag coefficient m_b is alternatively set to 0 and 0.45 to show the effects of advective transport in the bulk solution. The corresponding Dukhin number w is plotted on figure (c). The packing parameter is ν = 0.015, which corresponds to a bulk concentration c_{b, ∞} = 0.1 M for an ionic size a = 0.5 nm. The Debye length is thus λ_D = 0.96 nm. For ε = 0.01, the radius of the particle is R = 96 nm, so that the applied electric field is E* = 1.0 × 10⁴ V/cm.

FIG. 10.

View large Download slide

Ionic concentration field in the bulk solution for m_b = 0 (left) m_b = 0.45 (right), for ε = 0.01. The packing parameter is ν = 0.015, which corresponds to a bulk concentration c_{b, ∞} = 0.1M for an ionic diameter a = 0.5 nm. The Debye length is thus λ_D = 0.96 nm. For ε = 0.01, the radius of the particle is R = 96 nm, so that the applied electric field is E* = 1.0 × 10⁴ V/cm. Concentration gradients enhance advective transport in the bulk solution, which significantly reduces ionic concentration fluctuations in the electric double layer.

V. CONCLUSIONS

In this article, we have presented a numerical model aimed at describing nonlinear electrophoresis of ideally polarizable particles. We used the modified Poisson-Boltzmann equation obtained by Kilic et al.^39,40 to account for steric effects due to the finite size of ions in the electric double layer, and we modelled the ion fluxes in the bulk solution as a combination of diffusion, electromigration, and advection. We incorporated the model developed by Chu and Bazant³⁵ to match the double layer and the bulk solution asymptotically in a manner that accounts for surface conduction in the electric double layer.

The numerical simulations have demonstrated that this model, when applied to an ideally polarizable spherical particle, yields a non-zero electrophoretic mobility at high electric fields that scales approximately with the square root of the external electric field, in good agreement with the asymptotic formula derived by Bazant et al.³⁷ We have also quantified the influence of surface conduction on electrophoretic mobility (Fig. 5). The mechanism involved here is the emergence of concentration gradients in the bulk solution, which results in significant diffusio-osmotic flow around the particle. Hence, the diffusio-osmotic flow reduces slip velocity locally (Figs. 6 and 7). We also investigated the key role played by advective transport in the EDL, which was shown to enhance surface conduction. At high applied electric fields, surface conduction due to advective transport is comparable in magnitude to surface conduction due to chemical gradients (Fig. 8). Finally, we show that at high applied electric fields, advective transport in the bulk solution reduces concentration polarization, which in turn results in an increase of the electrophoretic mobility (Figs. 8–10).

The nonlinear relation between the applied electric field and the electrophoretic mobility could potentially open the way for new separation techniques and deposition methods for ideally polarizable particles. In the case of dielectric particles, Chimenti⁵³ and later Dukhin et al.⁵⁴ proposed a scheme for aqueous electrophoretic deposition in asymmetric AC electric fields, which relies on the nonlinear regime at high applied voltages. This scheme was successfully tested by Neirinck et al.⁵⁵ in 2009, and could easily be adapted to the case of ideally polarizable particles.

In this study, we have assumed the ionic diffusivities of the cations and anions to be the same. A possible extension of this work could be to account for the influence of the diffusivity ratio between cations and anions. Another assumption is that the solution permittivity remains constant across the electric double layer. However, several studies have recently suggested that the dielectric constant of the solvent decreases in regions with high ionic concentration due to the polarization of the solvent molecules caused by ions present.^56,57 This phenomenon has been termed excess ion polarization and can significantly influence the concentration profile in the EDL, as demonstrated by Hatlo et al.⁵⁸ Recently, Zhao and Zhai⁵⁹ have proposed a model describing the influence of dielectric decrement on several electrokinetic phenomena, including electro-osmosis and electrophoresis. In parallel, Bazant, Storey and co-workers have recently developed a model aimed at describing the physics of the electric double layer when the mean-field approximation breaks down. This model accounts for subtle electrokinetics effects occurring at large applied voltages.^60,61 A natural extension of our study could be the development of a numerical model accounting for dielectric decrement in the electric double layer surrounding the particle and for the break-down of the mean-field approximation at large applied electric fields.

ACKNOWLEDGMENTS

The authors wish to acknowledge one of the reviewers who noted an error in the original manuscript and derived the correct formula for Eq. (68), which incorporates the effects of the body force on the electrophoretic velocity. Thus Sec. III D largely resulted from the reviewers comments. The authors would also like to thank Professor M. Z. Bazant for his suggestions regarding the bibliography. The work presented in this paper was partially supported by the Shapiro Postdoctoral Fellowship (B.F. and J.L.M.) and an NSF Career Award (C.R.B.) under Grant No. 1150615.

REFERENCES

1.

T. M.

Squires

and

S. R.

Quake

, “

Microfluidics: Fluid physics at the nanoliter scale

,”

Rev. Mod. Phys.

77

,

977

–

1026

(

2005

).

https://doi.org/10.1103/RevModPhys.77.977

Google Scholar

Crossref

2.

H. A.

Stone

,

A. D.

Stroock

, and

A.

Ajdari

, “

Engineering flows in small devices: Microfluidics toward a lab-on-a-chip

,”

Annu. Rev. Fluid Mech.

36

,

381

–

411

(

2004

).

https://doi.org/10.1146/annurev.fluid.36.050802.122124

Google Scholar

Crossref

3.

J.-L.

Viovy

, “

Electrophoresis of DNA and other polyelectrolytes: Physical mechanisms

,”

Rev. Mod. Phys.

72

(

3

),

813

(

2000

).

https://doi.org/10.1103/RevModPhys.72.813

Google Scholar

Crossref

4.

M.

Von Smoluchowski

, “

Contribution à la théorie de l'endosmose électrique et de quelques phénomènes corrélatifs

,”

Bull. Acad. Sci. Cracovie

8

,

182

–

200

(

1903

).

https://doi.org/10.1051/jphystap:019040030091201

Google Scholar

Crossref

5.

F. A.

Morrison

Jr., “

Electrophoresis of a particle of arbitrary shape

,”

J. Colloid Interface Sci.

34

(

2

),

210

–

214

(

1970

).

https://doi.org/10.1016/0021-9797(70)90171-2

Google Scholar

Crossref

6.

J. L.

Anderson

, “

Colloid transport by interfacial forces

,”

Annu. Rev. Fluid Mech.

21

,

61

–

99

(

1989

).

https://doi.org/10.1146/annurev.fl.21.010189.000425

Google Scholar

Crossref

7.

A.

Sellier

, “

Sur l'électrophorèse d'un ensemble de particules portant la même densité uniforme de charges

,”

C.R. Acad. Sci., Ser. IIb: Mec., Phys., Chim., Astron.

327

(

5

),

443

–

448

(

1999

).

https://doi.org/10.1016/S1287-4620(99)80097-0

Google Scholar

Crossref

8.

O.

Schnitzer

and

E.

Yariv

, “

Macroscale description of electrokinetic flows at large zeta potentials: Nonlinear surface conduction

,”

Phys. Rev. E

86

(

2

),

021503

(

2012

).

https://doi.org/10.1103/PhysRevE.86.021503

Google Scholar

Crossref

9.

R. J.

Hunter

,

Foundations of Colloid Science (POD)

(

Oxford University Press

,

2000

).

Google Scholar

10.

H.

Lyklema

,

Introduction to Colloid Science

(

Elsevier

,

1995

).

Google Scholar

11.

J. J.

Bikerman

, “

Electrokinetic equations and surface conductance. a survey of the diffuse double layer theory of colloidal solutions

,”

Trans. Faraday Soc.

35

,

154

–

160

(

1940

).

https://doi.org/10.1039/tf9403500154

Google Scholar

Crossref

12.

J. Th. G.

Overbeek

, “

Theory of electrophoresis - the relaxation effect

,”

Kolloid-Beih.

54

,

287

–

364

(

1943

).

Google Scholar

Crossref

13.

S. S.

Dukhin

,“

Electrophoresis at large peclet numbers

,”

Adv. Colloid Interface Sci.

36

,

219

–

248

(

1991

).

https://doi.org/10.1016/0001-8686(91)80034-H

Google Scholar

Crossref

14.

B. V.

Derjaguin

,

S. S.

Dukhin

, and

V. N.

Shilov

, “

Kinetic aspects of electrochemistry of disperse systems. Part I. Introduction

,”

Adv. Colloid Interface Sci.

13

,

141

–

152

(

1980

).

https://doi.org/10.1016/0001-8686(80)87004-7

Google Scholar

Crossref

15.

S. S.

Dukhin

and

V. N.

Shilov

, “

Kinetic aspects of electrochemistry of disperse systems. Part II. Induced dipole moment and the non-equilibrium double layer of a colloid particle

,”

Adv. Colloid Interface Sci.

13

,

153

–

195

(

1980

).

https://doi.org/10.1016/0001-8686(80)87005-9

Google Scholar

Crossref

16.

S. S.

Dukhin

, “

Electrokinetic phenomena of the second kind and their applications

,”

Adv. Colloid Interface Sci.

35

,

173

–

196

(

1991

).

https://doi.org/10.1016/0001-8686(91)80022-C

Google Scholar

Crossref

17.

S. S.

Dukhin

and

B. V.

Derjaguin

,

Electrokinetic Phenomena

(

John Wiley & Sons

,

1974

).

Google Scholar

18.

S.

Barany

, “

Electrophoresis in strong electric fields

,”

Adv. Colloid Interface Sci.

147–148

,

36

–

43

(

2009

).

https://doi.org/10.1016/j.cis.2008.10.006

Google Scholar

Crossref

PubMed

19.

V.

Shilov

,

S.

Barany

,

C.

Grosse

, and

O.

Shramko

, “

Field-induced disturbance of the double layer electro-neutrality and non-linear electrophoresis

,”

Adv. Colloid Interface Sci.

104

,

159

–

173

(

2003

).

https://doi.org/10.1016/S0001-8686(03)00040-X

Google Scholar

Crossref

PubMed

20.

N. A.

Mishchuk

and

S. S.

Dukhin

, “

Electrophoresis of solid particles at large peclet numbers

,”

Electrophoresis

23

,

2012

–

2022

(

2002

).

https://doi.org/10.1002/1522-2683(200207)23:13<2012::AID-ELPS2012>3.0.CO;2-Y

Google Scholar

Crossref

PubMed

21.

R. W.

O'Brien

and

L. R.

White

, “

Electrophoretic mobility of a spherical colloidal particle

,”

J. Chem. Soc., Faraday Trans. 2

74

,

1607

–

1626

(

1978

).

https://doi.org/10.1039/f29787401607

Google Scholar

Crossref

22.

R. W.

O'Brien

and

R. J.

Hunter

, “

The electrophoretic mobility of large colloidal particles

,”

Can. J. Chem.

59

(

13

),

1878

–

1887

(

1981

).

https://doi.org/10.1139/v81-280

Google Scholar

Crossref

23.

R. W.

O'Brien

and

D. N.

Ward

, “

The electrophoresis of a spheroid with a thin double layer

,”

J. Colloid Interface Sci.

121

(

2

),

402

–

413

(

1988

).

https://doi.org/10.1016/0021-9797(88)90443-2

Google Scholar

Crossref

24.

D. C.

Prieve

,

J. L.

Anderson

,

J. P.

Ebel

, and

M. E.

Lowell

, “

Motion of a particle generated by chemical gradients. Part 2. Electrolytes

,”

J. Fluid Mech.

148

,

247

–

269

(

1984

).

https://doi.org/10.1017/S0022112084002330

Google Scholar

Crossref

25.

I.

Rubinstein

and

B.

Zaltzman

, “

Electro-osmotic slip of the second kind and instability in concentration polarization at electrodialysis membranes

,”

Math. Models Methods Appl. Sci.

11

(

02

),

263

–

300

(

2001

).

https://doi.org/10.1142/S0218202501000866

Google Scholar

Crossref

26.

B.

Zaltzman

and

I.

Rubinstein

, “

Electro-osmotic slip and electroconvective instability

,”

J. Fluid Mech.

579

,

173

–

226

(

2007

).

https://doi.org/10.1017/S0022112007004880

Google Scholar

Crossref

27.

R. A.

Rica

and

M. Z.

Bazant

, “

Electrodiffusiophoresis: Particle motion in electrolytes under direct current

,”

Phys. Fluids

22

,

112109

(

2010

).

https://doi.org/10.1063/1.3496976

Google Scholar

Crossref

28.

A. S.

Khair

, “

Transient phoretic migration of a permselective colloidal particle

,”

J. Colloid Interface Sci.

381

(

1

),

183

–

188

(

2012

).

https://doi.org/10.1016/j.jcis.2012.05.038

Google Scholar

Crossref

PubMed

29.

O.

Schnitzer

,

R.

Zeyde

,

I.

Yavneh

, and

E.

Yariv

, “

Weakly nonlinear electrophoresis of a highly charged colloidal particle

,”

Phys. Fluids

25

,

052004

(

2013

).

https://doi.org/10.1063/1.4804672

Google Scholar

Crossref

30.

O.

Schnitzer

and

E.

Yariv

, “

Strong-field electrophoresis

,”

J. Fluid Mech.

701

,

333

–

351

(

2012

).

https://doi.org/10.1017/jfm.2012.161

Google Scholar

Crossref

31.

A.

Ramos

,

H.

Morgan

,

N. G.

Green

, and

A.

Castellanos

, “

AC electric-field-induced fluid flow in microelectrodes

,”

J. Colloid Interface Sci.

217

(

2

),

420

–

422

(

1999

).

https://doi.org/10.1006/jcis.1999.6346

Google Scholar

Crossref

PubMed

32.

V. A.

Murtsovkin

, “

Nonlinear flows near polarized disperse particles

,”

Colloid J. Russ. Acad. Sci.

58

(

3

),

341

–

349

(

1996

).

Google Scholar

33.

T. M.

Squires

and

M. Z.

Bazant

, “

Induced-charge electro-osmosis

,”

J. Fluid Mech.

509

(

1

),

217

–

252

(

2004

).

https://doi.org/10.1017/S0022112004009309

Google Scholar

Crossref

34.

M. Z.

Bazant

and

T. M.

Squires

, “

Induced-charge electrokinetic phenomena: Theory and microfluidic applications

,”

Phys. Rev. Lett.

92

(

6

),

066101

(

2004

).

https://doi.org/10.1103/PhysRevLett.92.066101

Google Scholar

Crossref

PubMed

35.

K. T.

Chu

and

M. Z.

Bazant

, “

Nonlinear electrochemical relaxation around conductors

,”

Phys. Rev. E

74

,

011501

(

2006

).

https://doi.org/10.1103/PhysRevE.74.011501

Google Scholar

Crossref

36.

E.

Yariv

, “

Nonlinear electrophoresis of ideally polarizable particles

,”

Europhys. Lett.

82

,

54004

(

2008

).

https://doi.org/10.1209/0295-5075/82/54004

Google Scholar

Crossref

37.

M. Z.

Bazant

,

M. S.

Kilic

,

B. D.

Storey

, and

A.

Ajdari

, “

Towards and understanding of induced-charge electrokinetics at large applied voltages in concentrated solutions

,”

Adv. Colloid Interface Sci.

152

,

48

–

88

(

2009

).

https://doi.org/10.1016/j.cis.2009.10.001

Google Scholar

Crossref

PubMed

38.

J. J.

Bikerman

, “

Structure and capacity of electrical double layer

,”

Philos. Mag.

33

(

220

),

384

–

397

(

1942

).

https://doi.org/10.1080/14786444208520813

Google Scholar

Crossref

39.

M. S.

Kilic

,

M. Z.

Bazant

, and

A.

Ajdari

, “

Steric effects in the dynamics of electrolytes at large applied voltages. I. Double-layer charging

,”

Phys. Rev. E

75

,

021502

(

2007

).

https://doi.org/10.1103/PhysRevE.75.021502

Google Scholar

Crossref

40.

M. S.

Kilic

and

M. Z.

Bazant

, “

Steric effects in the dynamics of electrolytes at large applied voltages. II. Modified poisson-nernst-planck equations

,”

Phys. Rev. E

75

,

021503

(

2007

).

https://doi.org/10.1103/PhysRevE.75.021503

Google Scholar

Crossref

41.

O.

Schnitzer

and

E.

Yariv

, “

Induced-charge electro-osmosis beyond weak fields

,”

Phys. Rev. E

86

(

6

),

061506

(

2012

).

https://doi.org/10.1103/PhysRevE.86.061506

Google Scholar

Crossref

42.

M.

Gouy

, “

Sur la constitution de la charge électrique à la surface d'un électrolyte

,”

J. Phys. Théorique et Appliquée

9

(

1

),

457

–

467

(

1910

).

https://doi.org/10.1051/jphystap:019100090045700

Google Scholar

Crossref

43.

D. L.

Chapman

, “

A contribution to the theory of electrocapillarity

,”

Philos. Mag.

25

(

148

),

475

–

481

(

1913

).

https://doi.org/10.1080/14786440408634187

Google Scholar

Crossref

44.

O.

Stern

, “

The theory of the electrolytic double-layer

,”

Z. Elektrochem.

30

,

508

–

516

(

1924

).

Google Scholar

45.

A. S.

Khair

and

T. M.

Squires

, “

Ion steric effects on electrophoresis of a colloidal particle

,”

J. Fluid Mech.

640

,

343

(

2009

).

https://doi.org/10.1017/S0022112009991728

Google Scholar

Crossref

46.

E.

Yariv

, “

An asymptotic derivation of the thin-debye-layer limit for electrokinetic phenomena

,”

Chem. Eng. Commun.

197

(

1

),

3

–

17

(

2009

).

https://doi.org/10.1080/00986440903076590

Google Scholar

Crossref

47.

A. S.

Khair

and

T. M.

Squires

, “

Fundamental aspects of concentration polarization arising from nonuniform electrokinetic transport

,”

Phys. Fluids

20

(

8

),

087102

–

087102

(

2008

).

https://doi.org/10.1063/1.2963507

Google Scholar

Crossref

48.

M. S.

Kilic

, “

Induced-charge electrokinetics at large voltages

,” Ph.D. thesis (

Massachusetts Institute of Technology

,

2008

).

Google Scholar

49.

K. T.

Chu

and

M. Z.

Bazant

, “

Surface conservation laws at microscopically diffuse interfaces

,”

J. Colloid Interface Sci.

315

(

1

),

319

–

329

(

2007

).

https://doi.org/10.1016/j.jcis.2007.06.024

Google Scholar

Crossref

PubMed

50.

H. A.

Stone

and

A. D. T.

Samuel

, “

Propulsion of microorganisms by surface distortions

,”

Phys. Rev. Lett.

77

(

19

),

4102

–

4104

(

1996

).

https://doi.org/10.1103/PhysRevLett.77.4102

Google Scholar

Crossref

PubMed

51.

L. G.

Leal

,

Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes

(

Cambridge University Press

,

2007

).

Google Scholar

Crossref

52.

J. P.

Boyd

,

Chebyshev and Fourier Spectral Methods

(

Courier Dover Publications

,

2013

).

Google Scholar

53.

R. J. L.

Chimenti

, “

Electrophoretic separation

,” U.S. patent 5,106,468 (21 April

1992

).

Google Scholar

54.

A. S.

Dukhin

and

S. S.

Dukhin

, “

Aperiodic capillary electrophoresis method using an alternating current electric field for separation of macromolecules

,”

Electrophoresis

26

,

2149

–

2153

(

2005

).

https://doi.org/10.1002/elps.200410408

Google Scholar

Crossref

PubMed

55.

B.

Neirinck

,

J.

Fransaer

,

O.

Van der Biest

, and

J.

Vleugels

, “

Aqueous electrophoretic deposition in asymmetric ac electric fields

,”

Electrochem. Commun.

11

,

57

–

60

(

2009

).

https://doi.org/10.1016/j.elecom.2008.10.028

Google Scholar

Crossref

56.

D.

Ben-Yaakov

,

D.

Andelman

,

D.

Harries

, and

R.

Podgornik

, “

Beyond standard Poisson–Boltzmann theory: Ion-specific interactions in aqueous solutions

,”

J. Phys.: Condens. Matter

21

(

42

),

424106

(

2009

).

https://doi.org/10.1088/0953-8984/21/42/424106

Google Scholar

Crossref

PubMed

57.

D.

Ben-Yaakov

,

D.

Andelman

, and

R.

Podgornik

, “

Dielectric decrement as a source of ion-specific effects

,”

J. Chem. Phys.

134

(

7

),

074705

–

074705

(

2011

).

https://doi.org/10.1063/1.3549915

Google Scholar

Crossref

PubMed

58.

M. M.

Hatlo

,

R. H. H. G.

Van Roij

, and

L.

Lue

, “

The electric double layer at high surface potentials: The influence of excess ion polarizability

,”

Europhys. Lett.

97

(

2

),

28010

(

2012

).

https://doi.org/10.1209/0295-5075/97/28010

Google Scholar

Crossref

59.

H.

Zhao

and

S.

Zhai

, “

The influence of dielectric decrement on electrokinetics

,”

J. Fluid Mech.

724

,

69

–

94

(

2013

).

https://doi.org/10.1017/jfm.2013.152

Google Scholar

Crossref

PubMed

60.

M. Z.

Bazant

,

B. D.

Storey

, and

A. A.

Kornyshev

, “

Double layer in ionic liquids: Overscreening versus crowding

,”

Phys. Rev. Lett.

106

(

4

),

046102

(

2011

).

https://doi.org/10.1103/PhysRevLett.106.046102

Google Scholar

Crossref

PubMed

61.

B. D.

Storey

and

M. Z.

Bazant

, “

Effects of electrostatic correlations on electrokinetic phenomena

,”

Phys. Rev. E

86

(

5

),

056303

(

2012

).

https://doi.org/10.1103/PhysRevE.86.056303

Google Scholar

Crossref

2014

AIP Publishing LLC

Nonlinear electrophoresis of ideally polarizable particles

I. INTRODUCTION

II. ELECTRIC DOUBLE LAYER

A. Electric potential in the electric double layer

B. Electroosmotic flow in the electric double layer

III. BULK EQUATIONS

A. Electric potential and ionic concentration

B. Determination of the zeta potential

C. Velocity field and electrophoretic velocity

D. Electrophoretic velocity

IV. RESULTS AND DISCUSSION

A. Numerical model

B. Weakly nonlinear regime

C. Crowding effects

D. Influence of surface conduction

V. CONCLUSIONS

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Nonlinear electrophoresis of ideally polarizable particles

I. INTRODUCTION

II. ELECTRIC DOUBLE LAYER

A. Electric potential in the electric double layer

B. Electroosmotic flow in the electric double layer

III. BULK EQUATIONS

A. Electric potential and ionic concentration

B. Determination of the zeta potential

C. Velocity field and electrophoretic velocity

D. Electrophoretic velocity

IV. RESULTS AND DISCUSSION

A. Numerical model

B. Weakly nonlinear regime

C. Crowding effects

D. Influence of surface conduction

V. CONCLUSIONS

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only