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(

## 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.^{3} The mathematical description of electrophoresis of particles with thin electric double layer builds upon the classical theory developed by Smoluchowski^{4} 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*

where μ_{EP} denotes the electrophoretic mobility of the particle. Later, Morrison^{5} and Anderson^{6} demonstrated that Smoluchowski's result is also valid for single non-spherical particles, and Sellier^{7} 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^{8} 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^{11} 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,^{12} 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.^{24} 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^{27} 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,^{28} 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.^{8}

Nonlinear electrokinetics experienced a renewed interest in the mid 1990s, when various nonlinear electrokinetic phenomena were discovered and explored, including AC-electro-osmosis^{31} 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,^{35} who conducted numerical simulations to determine the electric potential and the ionic concentration around a conductor. Yariv^{36} in 2008 and Bazant *et al.*^{37} 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,^{36} 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,^{41} 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^{42} and in 1913 by Chapman,^{43} 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.^{44} As reviewed by Bazant *et al.*,^{37} 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,^{37} 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,^{36} 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^{45} 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

_{b}between the electric double layer and the bulk solution, and of the local ionic concentrations

*c*

_{+}and

*c*

_{−}:

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^{40}

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

*a*and to the bulk solution concentration

*c*

_{b,∞}far from the particle through the estimate

^{39}

Setting the variational derivatives

Here, *c*_{b} denotes the bulk concentration in the immediate vicinity of the electric double layer. As reviewed by Bazant,^{37} Eq. (5) was derived for the first time by Bikerman^{38} 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

For such particles, one can demonstrate through a rigorous matched asymptotic analysis^{46} 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

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.

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

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

This equation can be integrated once to yield,^{39} in dimensionless form,

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*),

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

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

The projection of the Stokes equation on **e**_{x} reads

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

Interestingly, as noted by Prieve *et al.*,^{24} 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

where the dimensionless coefficient κ(ζ) is given by

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^{17} 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.^{27} A formulation for coefficient (17) was first obtained by Kilic^{48} 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

With this scaling, we find that the dimensionless slip velocity *u*_{s} at the surface of the particle is

## 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

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

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

we obtain the dimensionless Nernst-Planck relation

where m_{b} is a dimensionless ionic drag coefficient defined for the bulk solution by

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

which can be combined to yield, in dimensionless form,

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.^{35}

The inner and the outer asymptotic solutions are matched by appropriate boundary conditions expressing the conservation of ions across the interface. Chu and Bazant^{49} 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:

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,

In this relation, Γ_{±} denotes the surface adsorption coefficient,^{39} defined by

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

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

so that

where div_{S} denotes the surface divergence operator. In outer spherical coordinates, these relations become

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

and

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

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

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

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

we obtain the dimensionless expression

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

we find

The dimensionless boundary conditions (33) become then

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:

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

where the potential Ψ_{P} of the particle is unknown. To evaluate the electric potential of the particle, as demonstrated by Yariv,^{36} 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

*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

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

### 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

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}∇^{2}Ψ_{b}∇Ψ_{b} refers to the effects of the body force. The Stokes equation is considered along with the continuity equation

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

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

Taking the curl of this relation, we finally obtain

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

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

where operator

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

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):

### 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,^{50} 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:

A similar expression was obtained by Anderson^{6} 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

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

To facilitate the determination of the hydrodynamic force, we can make use of the reciprocal theorem.^{51} To that end, let

**e**

_{z}in an unbounded fluid at rest at infinity. For the companion flow, the body force vanishes. The velocity field

The reciprocal theorem states that

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.^{51} Thus, at the surface of the particle, the traction yields

so that

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

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

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

## 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^{35} 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.^{52} 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:

In this relation, ζ_{0} 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

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^{−3}. For an electrolyte solution at 25 °C where *c*_{b, ∞} = 1.0 × 10^{−3} M, the Debye length is λ_{D} = 9.6 nm, which results in ε ≃ 1.0 × 10^{−2} 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^{−9} m^{2}/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^{−5} 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^{36} 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^{36}

As a consequence, the zeta potential is

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

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)^{36}

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.

### 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^{36} 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.

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

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

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

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

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

for reasonably large ν, to yield

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

This asymptotic formula was first found and interpreted by Bazant *et al.*^{37} 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.

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^{30} 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.

## 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^{35} 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.*^{37} 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^{53} and later Dukhin *et al.*^{54} 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.*^{55} 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.*^{58} Recently, Zhao and Zhai^{59} 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.