2. Problem formulation

We consider a dielectric sphere with negative surface charge density $q_s$ , suspended in a binary electrolyte with ionic diffusivities $D_+$ and $D_-$ , and subjected to a DC electric field and a linear shear flow (figure 1 c). Our goal is to determine its electrophoretic velocity. To this end, we adopt the electrokinetic framework developed by Schnitzer and Yariv for large zeta potentials (Schnitzer & Yariv Reference Schnitzer and Yariv2012) (hereafter, the SY model), which has been successfully applied to electrophoresis (Schnitzer et al. Reference Schnitzer, Zeyde, Yavneh and Yariv2013; Schnitzer & Yariv Reference Schnitzer and Yariv2014) and CP electroosmotic flows (Calero et al. Reference Calero, Fernández-Mateo, Morgan, García-Sánchez and Ramos2021; Fernández-Mateo et al. Reference Fernández-Mateo, García-Sánchez, Calero, Morgan and Ramos2021). The SY model provides the governing equations for the electric potential $\phi$ , ionic concentration $c$ and fluid velocity $\boldsymbol {u}$ in the bulk electrolyte under the thin EDL approximation, appropriate for the experiments mentioned above where microscopic particles are dispersed in electrolytes with typical values of the Debye length in the range 10–100 nm.

We adopt a reference frame fixed at the centre of mass of the particle, and express the governing equations in non-dimensional form using the following characteristic scales: the particle radius $a$ as the length scale; the bulk ionic concentration $c_0$ as the concentration scale; and the thermal voltage for the electric potential, $V_{{th}}=k_{{B}} T/{ze}$ , where $T$ is the absolute temperature, $k_{ B}$ is the Boltzmann constant, $z$ is the ionic valence and $e$ is the elementary charge. For a $1:1$ electrolyte at room temperature, the thermal voltage is $V_{{th}}\approx 25\,\mathrm{mV}$ . Velocities are scaled by $\varepsilon V_{{th}}^2/(a\eta )$ , which in turn defines the characteristic time scale as $a^2\eta /(\varepsilon V_{{th}}^2)$ . Consequently, diffusion coefficients are scaled by $\varepsilon V^2_{{th}}/\eta$ . The dimensionless ambipolar diffusion coefficient is $D=2D_+D_-/(D_++D_-)$ . The non-dimensional governing equations are

(2.1) \begin{eqnarray} D{\nabla} ^2c-\boldsymbol {u}\boldsymbol{\cdot }\boldsymbol{\nabla }c&=&0, \end{eqnarray}

(2.2) \begin{eqnarray} \boldsymbol{\nabla }\boldsymbol{\cdot }(c\boldsymbol{\nabla }\phi )-\alpha \boldsymbol {u}\boldsymbol{\cdot }\boldsymbol{\nabla }c&=&0, \end{eqnarray}

(2.3) \begin{eqnarray} -\boldsymbol{\nabla }\!p+ {\nabla} ^2{\boldsymbol{u}}+{\nabla} ^2\phi \boldsymbol{\nabla }\phi =0,\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}&=&0. \end{eqnarray}

Equation (2.1) is the convection–diffusion equation for the ion concentration, while the equation for the electric potential (2.2) corresponds to the continuity of the electrical current with $\alpha \equiv(D_+^{-1}-D_-^{-1})/2$ . Equations (2.3) are the Stokes equations for the fluid velocity ${\boldsymbol{u}}$ and pressure $p$ (negligible Reynolds number) including a body force term of electrical origin.

The SY model also establishes the boundary conditions (BCs) at the particle surface. The model assumes that the particle zeta potential is much larger than the thermal voltage, i.e. $\zeta _0 \gg 25\,\mathrm{mV}$ , an assumption that is well satisfied in the experiments described above involving carboxylate microspheres. Under these conditions, co-ions are effectively excluded from the diffuse layer, while the counter-ion flux normal to the surface is balanced by the surface current within the EDL. As a result, the following BCs are obtained:

(2.4) \begin{equation} \left (\frac {\partial c}{\partial n}- c\frac {\partial \phi }{\partial n}\right ) = 0\quad \mathrm{and}\quad \left (\frac {\partial c}{\partial n}+ c\frac {\partial \phi }{\partial n}\right ) = -2{\textit{Du}}\boldsymbol{\nabla} _s^2\left (\phi +\ln c\right ), \quad \mathrm{at}\quad r=1, \end{equation}

where $\boldsymbol{\nabla} _s$ denotes the surface gradient operator tangent to the particle, and ${\textit{Du}}$ is the dimensionless Dukhin number, which quantifies the relative contribution of surface to bulk currents (Delgado et al. Reference Delgado, Gonzalez-Caballero, Hunter, Koopal and Lyklema2005).

Concerning the fluid velocity at the particle surface, the SY model prescribes the following slip BC, which coincides with the expression derived in (Prieve et al. Reference Prieve, Anderson, Ebel and Lowell1984) in the context of diffusiophoresis:

(2.5) \begin{eqnarray} {\boldsymbol {u}}^s&=&\zeta \boldsymbol{\nabla} _s\phi - 4\ln \left [\cosh (\zeta /4)\right ]\boldsymbol{\nabla} _s \ln c, \end{eqnarray}

where the non-dimensional zeta potential $\zeta$ is related to the non-dimensional surface charge density, $q_s$ , and ion concentration through the Gouy–Chapman equation $q_s=2\sqrt {c}\sinh {(\zeta /2)}$ (Hunter Reference Hunter1993). The BC for the fluid velocity at the particle surface is

(2.6) \begin{eqnarray} {\boldsymbol {u}}={\boldsymbol {u}}^s+\boldsymbol{\varOmega }\times \boldsymbol {r} \quad \mathrm{at}\quad r=1 ,\end{eqnarray}

where $\boldsymbol{\varOmega }$ is the angular velocity of the spherical particle.

In the far field ( $\boldsymbol{r}\rightarrow \infty$ ), the electric potential and ionic concentration approach their imposed values, $\phi =-\beta z$ and $c=1$ , with $\beta =E_0(a/V_{{th}})$ . Since the reference frame is attached to the particle, the fluid velocity approaches $\boldsymbol{u}=\Gamma x\hat {\boldsymbol{z}}-\boldsymbol{U}$ far from the particle, where $\Gamma =(\varepsilon V_{{th}}^2/a^2\eta )\dot {\gamma }$ and $\boldsymbol{U}$ is the electrophoretic velocity. The translational and angular velocities, $\boldsymbol{U}$ and $\boldsymbol{\varOmega }$ , are determined by enforcing zero net hydrodynamic force and torque, $\oint _{S} \tilde {\boldsymbol{T}}\boldsymbol{\cdot }{\rm d}\boldsymbol{S} = 0$ and $\oint _{S} (\boldsymbol {r}\times \tilde {\boldsymbol{T}})\boldsymbol{\cdot }{\rm d}\boldsymbol{S} = 0$ , where $\tilde {\boldsymbol{T}}$ is the hydrodynamic stress tensor and the integrals are evaluated over a closed surface enclosing the sphere and located sufficiently far from it.

Although the SY model was originally derived for a particle in a quiescent fluid, an imposed flow could, in principle, affect ion transport within the EDL and modify the BCs (2.4). Nevertheless, following the analysis in Yariv et al. (Reference Yariv, Schnitzer and Frankel2011) for the streaming potential generated by a falling sphere, the charge-advection contribution scales as $q_s u \delta$ , with $\delta$ the dimensionless Debye length. In the thin-EDL limit and under typical experimental conditions, this term is negligible compared with the electromigration current associated with surface conductance. Hence, the BCs (2.4) remain unchanged.

3. Problem analysis and perturbation scheme

We seek an analytical expression for the particle velocity in the limit of weak electric and shear forcing. To this end, the governing fields are expanded asymptotically in powers of the dimensionless electric-field strength, $\beta$ , and the dimensionless shear rate, $\Gamma$ ,

(3.1) \begin{eqnarray} c&=&1+\beta c_{01}+\Gamma \beta c_{11}+\ldots , \end{eqnarray}

(3.2) \begin{eqnarray} \phi &=&\beta \phi _{01}+\Gamma \beta \phi _{11}+\ldots , \end{eqnarray}

(3.3) \begin{eqnarray} \boldsymbol{u}&=&\beta \boldsymbol{u}_{01}+\Gamma \boldsymbol{u}_{10}+\Gamma \beta \boldsymbol{u}_{11}+\ldots , \end{eqnarray}

where the first subscript denotes the order of the expansion in $\Gamma$ , while the second denotes the order in $\beta$ . The fields $\phi _{10}$ and $c_{10}$ are omitted here since, in the thin-EDL limit and without an applied electric field, the shear flow does not perturb the electric potential or the ionic concentration, as discussed in the previous section (Yariv et al. Reference Yariv, Schnitzer and Frankel2011). For the same reason, terms proportional to $\Gamma ^2$ do not appear in our second-order expansion. Likewise, higher-order electric-field contributions do not affect lateral migration because the system remains axially symmetric around the sphere in the absence of shear flow. Consequently, our perturbation expansion up to second order retains only the velocity term proportional to $\Gamma \beta$ .

3.1. Solutions at first order in the electric field and zeroth order in the shear flow

The equations for $c_{01}$ and $\phi _{01}$ are, respectively,

(3.4) \begin{equation} D{\nabla} ^2c_{01} = 0\quad \mathrm{and}\quad {\nabla} ^2\phi _{01} = 0, \end{equation}

with BCs at the particle surface, $r=1$ ,

(3.5) \begin{equation} \frac {\partial \phi _{01}}{\partial r} - \frac {\partial c_{01}}{\partial r} = 0, \quad \frac {\partial c_{01}}{\partial r} = -{\textit{Du}}\boldsymbol{\nabla} _s^2(\phi _{01} + c_{01}). \end{equation}

BCs far from the particle are ( $r \to \infty$ ): $c_{01} \to 0, \quad \phi _{01} \to -r \cos \theta$ .

The solutions for $c_{01}$ and $\phi _{01}$ are

(3.6) \begin{equation} c_{01}=\frac {3 {\textit{Du}}}{2+4{\textit{Du}}}\frac {\cos \theta }{r^2},\qquad \phi _{01}=-r\cos \theta +\frac {{\textit{Du}}-1}{2+4{\textit{Du}}}\frac {\cos \theta }{r^2}, \end{equation}

which are coincident with the expressions in Schnitzer & Yariv (Reference Schnitzer and Yariv2012) and in Hunter (Reference Hunter1993) for a thin EDL and a highly charged surface.

The velocity field $\boldsymbol{u}_{01}$ is governed by the equations

(3.7) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_{01} = 0, \quad {\nabla} ^2\boldsymbol{u}_{01} = \boldsymbol{\nabla }\!p_{01}, \end{equation}

with BC at the particle surface: $ \boldsymbol{u}^s_{01} = \zeta _0\boldsymbol{\nabla} _s\phi _{01} - 4\ln (\cosh (\zeta _0/4) )\boldsymbol{\nabla} _s c_{01}$ .

Far from the particle, the velocity field approaches the negative electrophoretic velocity, $\boldsymbol{u}_{01} \to -\boldsymbol{U}_{01}$ , which is determined by the zero-net-force condition on the particle. The solution for $\boldsymbol{U}_{01}$ is

(3.8) \begin{equation} U_{01}=\frac {\zeta _0(1+{\textit{Du}})+4{\textit{Du}}\log \left (\cosh (\zeta _0/4)\right )}{1+2{\textit{Du}}}. \end{equation}

This expression is in agreement with the results given by Schnitzer et al. (Reference Schnitzer, Zeyde, Yavneh and Yariv2013).

3.2. Solutions at zeroth order in the electric field and first order in the shear flow

At this order, the velocity field $\boldsymbol{u}_{10}$ is the only field that remains to be solved, and it is governed by the Stokes equations

(3.9) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_{10} = 0, \quad {\nabla} ^2\boldsymbol{u}_{10} = \boldsymbol{\nabla }\!p_{10}. \end{equation}

Far from the particle, the velocity approaches the undisturbed shear flow: $ \boldsymbol{u}_{10} \to x\hat {\boldsymbol{z}}$ as $r \to \infty$ . The following no-slip BC applies at the particle surface: $ \boldsymbol{u}_{10} = \boldsymbol{\varOmega }\times \boldsymbol{r}$ , at $r=1$ , where the angular velocity $\boldsymbol{\varOmega }$ is determined by the hydrodynamic zero-net-torque condition on the particle. The solution for this velocity field is well known (Leal Reference Leal2007, ch. 8) and can be expressed as

(3.10) \begin{equation} \boldsymbol{u}_{10} = x\hat {\boldsymbol{z}} - \left [ \frac {x\hat {\boldsymbol{z}} + z\hat {\boldsymbol{x}}}{r^5} + \frac {5}{2r^3}\left (1-\frac {1}{r^2}\right )\frac {zx}{r^2}\boldsymbol{r} \right ]\!, \end{equation}

with the particle rotation given by $\boldsymbol{\varOmega } = -(1/2)\hat {\boldsymbol{y}}$ .

3.3. Solutions at first order in the electric field and first order in the shear flow

Equations governing $c_{11}$ and $\phi _{11}$ are

(3.11) \begin{eqnarray} D{\nabla} ^2c_{11} &=& \boldsymbol{u}_{10}\boldsymbol{\cdot }\boldsymbol{\nabla }c_{01}, \end{eqnarray}

(3.12) \begin{eqnarray} {\nabla} ^2\phi _{11} &=& \alpha \boldsymbol{u}_{10}\boldsymbol{\cdot }\boldsymbol{\nabla }c_{01}, \end{eqnarray}

with BCs at the particle surface, $r=1$ ,

(3.13) \begin{equation} \frac {\partial \phi _{11}}{\partial r} - \frac {\partial c_{11}}{\partial r} = 0, \quad \frac {\partial c_{11}}{\partial r} = -{\textit{Du}}\boldsymbol{\nabla} _s^2(\phi _{11} + c_{11}). \end{equation}

Far from the particle, $r \to \infty$ : $c_{11} \to 0$ and $\phi _{11} \to 0$ .

Using the results from the previous sections, the source term in (3.11) can be expressed in terms of the following associated Legendre polynomials of order 1, $P^1_1(\cos (\theta ))$ and $P^1_3(\cos (\theta ))$ :

(3.14) \begin{equation} \begin{aligned} \boldsymbol{u}_{10}\boldsymbol{\cdot }\boldsymbol{\nabla }c_{01} = & A\frac {\cos \varphi }{r^2} \left [\left ( -\frac {2}{5}-\frac {1}{r^3}+\frac {9}{10r^5} \right )P^1_1+\left ( \frac {2}{5}-\frac {2}{3r^3}+\frac {4}{15r^5} \right )P^1_3 \right]\!, \end{aligned} \end{equation}

with $A=({1}/{D})({3{\textit{Du}}}/{(2+4{\textit{Du}})})$ . Because the source term in (3.14) decays as $\mathcal{O}(r^{-2})$ , the solution of (3.11) cannot satisfy the far-field BC: the general solution for $c_{11}$ presented below does not vanish as $r\to \infty$

(3.15) \begin{align} c_{11} = -\left [f(r) + \frac {A_1}{r^2} + A_2 r\right ]\cos \varphi P_1^1(\cos \theta ) + \left [g(r) + \frac {B_1}{r^4} + B_2 r^3\right ]\cos \varphi P_3^1(\cos \theta ), \end{align}

where $f(r) = A (-({1}/{5}) + ({1}/{4r^3}) - ({1}/{20r^5}) )$ , $g(r) = A (-({1}/{30}) + ({1}/{9r^3}) + ({1}/{30r^5}) )$ and $A_1,A_2,B_1$ and $B_2$ are integration constants.

For (3.12), this issue can be circumvented by introducing the change of variables $\psi _{11}=\phi _{11}-\alpha Dc_{11}$ . Consequently, (3.12) reduces to ${\nabla} ^2\psi _{11}=0$ , which ensures that $\psi _{11}$ vanishes as $r\to \infty$ . BCs for $\psi _{11}$ on the particle surface are

(3.16) \begin{equation} \frac {\partial \psi _{11}}{\partial r}+(\alpha D-1) \frac {\partial c_{11}}{\partial r} = 0, \quad \frac {\partial c_{11}}{\partial r} = -{\textit{Du}}\boldsymbol{\nabla} _s^2(\psi _{11} + (\alpha D+1)c_{11}), \end{equation}

while it vanishes far from the particle ( $\psi _{11}\to 0$ as $r\to \infty$ ).

The breakdown of the regular perturbation scheme for $c_{11}$ is analogous to the Whitehead paradox for heat transfer from a sphere in a uniform flow, where no uniformly valid solution exists for the temperature field in the limit of ${\textit{Pe}} \to 0$ (Leal Reference Leal2007, ch. 9). Here, ${\textit{Pe}}$ denotes the Péclet number, which quantifies the ratio of convective to diffusive transport. In our problem, ${\textit{Pe}}=\Gamma /D$ when defined using the particle radius as the characteristic length scale. This paradox arises because the local Péclet number (which in our problem scales as ${\textit{Pe}}(r) = \Gamma r^2 / D$ ) cannot remain arbitrarily small throughout an unbounded domain. Consequently, a proper solution for $c_{11}$ and $\psi _{11}$ – including the unambiguous determination of the integration constants – must be sought within the framework of a singular perturbation analysis using the method of matched asymptotic expansions, as detailed in the next section.

Before proceeding with the asymptotic matching, we present the general solution for  $\psi _{11}$

(3.17) \begin{equation} {\psi _{11} = -\frac {C_1}{r^2}\cos \varphi P_1^1(\cos \theta ) + \frac {D_1}{r^4}\cos \varphi P_3^1(\cos \theta ),} \end{equation}

where the constants $C_1$ and $D_1$ are coupled to those in (3.15) via the BCs at $r=1$ (3.16)

(3.18) \begin{eqnarray} -2C_1 + (\alpha D-1)(f'(1) - 2A_1 + A_2) &=& 0, \end{eqnarray}

(3.19) \begin{eqnarray} -4D_1 +(\alpha D-1)( g'(1) - 4B_1 + 3B_2) &=& 0, \end{eqnarray}

(3.20) \begin{eqnarray} f'(1) - 2A_1 + A_2 &=& 2{\textit{Du}}(C_1 +(\alpha D+1) (f(1) + A_1 + A_2)), \end{eqnarray}

(3.21) \begin{eqnarray} g'(1) - 4B_1 + 3B_2 &=& 12{\textit{Du}}(D_1 + (\alpha D+1) (g(1) + B_1 + B_2)). \end{eqnarray}

4. Far-field analysis and asymptotic matching

As mentioned above, the convective term in (2.1) cannot be neglected at arbitrarily large length scales, as it becomes comparable to or dominant over diffusion in regions where $r \ge {\textit{Pe}}^{-1/2}$ . Consequently, we must seek an outer solution, denoted as $C(\boldsymbol{r})$ , valid in this far-field domain. This solution will subsequently be matched asymptotically with the near-field (inner region) solution, which holds where $r \ll {\textit{Pe}}^{-1/2}$ . To this end, we introduce the rescaled coordinate $\bar {r} = r \sqrt {{\textit{Pe}}}$ . This scaling yields a distinguished limit that balances the diffusive and convective terms in the equation for $C(\boldsymbol{r})$ as ${\textit{Pe}} \to 0$ (Leal Reference Leal2007, ch. 9.F).

Under the coordinate transformation $\bar {r} = r \sqrt {{\textit{Pe}}}$ , the conservation equation (2.1) at order $\beta$ becomes

(4.1) \begin{equation} {\nabla} ^2 C - \bar x\frac {\partial C}{\partial \bar z}=\left ( -{\textit{Pe}}^{3/2}\frac {5}{2}\frac {\bar z\bar x}{\bar r^5}\boldsymbol{\bar{{r}}} - {\textit{Pe}}^{5/2}\frac {5}{2}\frac {\bar z\bar x}{\bar r^7}\boldsymbol{\bar{{r}}} + {\textit{Pe}}^{5/2}\frac {\bar x\hat z +\bar z\hat x}{2\bar r^5}\right ) \cdot \bar {\nabla }C. \end{equation}

By neglecting higher-order terms in the limit ${\textit{Pe}} \to 0$ , we obtain the leading-order governing equation in the outer region

(4.2) \begin{equation} \bar {\nabla} ^2 C_0 -\bar x \frac {\partial C_0}{\partial \bar z} = 0. \end{equation}

The asymptotic matching condition requires that the inner limit of the outer solution matches the outer limit of the inner solution, namely, $\lim _{\bar {r} \to 0} C_0(\bar {r}) = \lim _{r \to \infty } c(r)$ .

A general solution to (4.2) that satisfies the far-field BC ( $\lim _{\bar {r} \to \infty }C(\bar {r}) = 0$ ) was derived by Elrick (Leal Reference Leal2007, ch. 9.F) in the context of the temperature distribution produced by a point heat source in a simple shear flow. This solution reads as follows:

(4.3) \begin{equation} T_0 (\bar {\boldsymbol{r}})= \frac {M}{2\sqrt {\pi }}\int _{0}^{\infty } \frac {\text{d}s}{(1+s^2/12)^{1/2}s^{3/2}}\exp \left [-\frac {(\bar z - \bar xs/2)^2}{4s(1+s^2/12)} - \frac {\bar y^2+\bar x^2}{4s}\right ]\!, \end{equation}

where $M$ is a constant and $s$ a dummy variable of integration. In our case, given the far-field behaviour of $c_{01}$ ((3.6)) the source for $C_0$ corresponds to a point dipole located at the origin and oriented along the $z$ -direction. By analogy with electrostatics, where a dipole field is related to a monopole field via a directional derivative, the solution to our problem can be obtained directly from (4.3) as $C_0(\bar {\boldsymbol{r}}) = -\partial T_0 / \partial \bar {z}$

(4.4) \begin{equation} C_0 (\bar {\boldsymbol{r}}) = \frac {M}{4\sqrt {\pi }}\int _{0}^{\infty } \frac {(\bar z - \bar xs/2)\text{d}s}{(1+s^2/12)^{3/2}s^{5/2}}\exp \left [-\frac {(\bar z - \bar xs/2)^2}{4s(1+s^2/12)} - \frac {\bar y^2+\bar x^2}{4s}\right ]\!. \end{equation}

In the limit $r \to 0$ , the dominant contribution to the integral in (4.4) arises from the region where $s \sim \bar r^2 \to 0$ . Consequently, the asymptotic behaviour of (4.4) near the origin can be determined analytically by expanding the integrand for small values of $s$ . It is straightforward to show that $C_0(\bar {\boldsymbol{r}}) \approx M ( \bar z/\bar r^3 - \bar x/4\bar r + \bar z^2\bar x/4\bar r^3)$ , where we have retained terms that, upon reverting the change of variables, yield contributions up to order $\textit{Pe}$ . Indeed, the evaluation yields

(4.5) \begin{align} C_0(\boldsymbol{r}) = M\left ({\textit{Pe}}^{-1}\frac {z}{r^3} - \frac {x}{4r} + \frac {z^2x}{4r^3}\right ) = \frac {M}{\textit{Pe}}\frac {\cos \theta }{r^2} - \frac {M}{4}\cos \varphi \left (-\frac {4}{5}P_1^1(\cos \theta ) + \frac {2}{15}P_3^1(\cos \theta )\right )\!. \end{align}

By comparing this result with the $c_{01}$ term of the inner solution, we obtain $M/\textit{Pe}=\textit{DA} = ({3\textit{Du}}/{(2+4\textit{Du})})$ . Asymptotic matching up to order $\textit{Pe}$ then requires that $\textit{Pe}\,c_{11}|_{r \gg 1} = (C_0 - \textit{DA}\cos \theta /\bar {r}^2 )|_{\bar {r} \ll 1}$ , which implies that $A_2 = B_2 = 0$ .