2.2.1. The characteristic scales and the dimensionless parameters

For the present scenario, it will generally be insightful and convenient to work with dimensionless (unit-free) quantities, thus it is important to first consider the various characteristic scales and non-dimensional parameters that are relevant to the system under consideration. To this end, for any general variable, say $\xi '$, we will denote the corresponding characteristic value as $\xi _{c}$.

Table 1 lists all the characteristic scales chosen for the present analysis. One important point to note here is that based on the characteristic surface charge ($\sigma _c$, assumed to be a known quantity), we have defined a characteristic potential $\zeta _c$ – this potential indicates the order of magnitude of the ‘zeta potential’ on the particle surface. In other words, $\zeta _c$ gives an estimate of how much the potential drops across the diffuse part of the EDL. Thus both $\sigma _c$ and $\zeta _c$ quantify how strong the surface charge is on the particle. When a particle is uniformly charged, the choice of $\sigma _c$ is trivial; for non-uniformly charged particles, the surface average of $\sigma '(\theta,\phi )$ may be chosen as $\sigma _c$.

Table 1.

Characteristic scales chosen (Saville Reference Saville1977; Ajdari Reference Ajdari1996) to non-dimensionalize pertinent variables, equations and the boundary conditions.

Based on the above characteristic scales, we may define several non-dimensional numbers that are relevant to the present analysis: (i) the non-dimensional EDL thickness $\delta = \lambda _{D}/a$; (ii) the characteristic potential on the particle surface relative to the thermal potential (Ajdari Reference Ajdari1996; Ghosh et al. Reference Ghosh, Mukherjee and Chakraborty2021) $\bar {\zeta }_0 = e\zeta _{c} /kT$; (iii) the externally imposed field strength relative to the thermal potential (Saville Reference Saville1977) $\beta = eE_{0}a/kT$; (iv) the nominal Deborah number, which quantifies the extent of viscoelasticity (Bird et al. Reference Bird, Armstrong and Hassager1987; Ghosh et al. Reference Ghosh, Mukherjee and Chakraborty2021), $De = u_{c}\lambda _c/a$.

2.2.2. The key assumptions

The present analysis is predicated on a few assumptions that simplify the governing equations (discussed later), yet without sacrificing the essential physics. We discuss some of those main assumptions here.

First, we will assume the particle's surface charge density to be ‘weak’, which necessitates $\bar {\zeta }_0 \ll 1$. Second, we will confine ourselves to the thin EDL limit, which mandates that the characteristic EDL thickness $\lambda _{D} = \sqrt {\epsilon k_b T/2e^2 c'_0}$ is small compared to the particle radius (Ajdari Reference Ajdari1995), i.e. $\delta \ll 1$ – this condition is usually satisfied for most particles that undergo electrophoretic motion. Third, the nominal Deborah number $De$ is assumed to remain $O(1)$ or less. However, because of the weak surface charge ($\bar {\zeta }_0\ll 1$), the (dimensionless) velocity would scale as $O(\bar {\zeta }_0)$ and hence the effective Deborah number (Ghosh et al. Reference Ghosh, Mukherjee and Chakraborty2021) $De_{eff}$ will become $O(\bar {\zeta }_0 De)$ $\sim O(\bar {\zeta }_0) \ll 1$, which indicates that the flow is only weakly viscoelastic despite the nominal $De$ being $O(1)$. Fourth, the imposed electric field strength is assumed to be low to moderate, so that $\beta \sim O(1)$.

Fifth, the flow field is assumed to be quasi-steady – i.e. the angular and the translational velocities of the particle are governed by the instantaneous surface charge distribution. Since the true velocity in the surrounding fluid scales as $O(\bar {\zeta }_0 u_c)$, the angular velocity of the particle is expected to scale as $\varOmega _c \sim \bar {\zeta }_0 u_c/a$, thus the time scale associated with the change in the surface charge distribution becomes $\varOmega _c^{-1} \sim a/\bar {\zeta }_0 u_c$. Therefore, for the quasi-static assumption to be satisfied, one must ensure that this time scale is large compared to the characteristic time scale of the fluid itself, i.e. $\varOmega _c^{-1} \gg \lambda _c$, which ultimately requires $De_{eff} \ll 1$ – this is consistent with the third assumption stated above. As a result, the flow can adjust to any change in the surface charge (or potential) distribution almost instantaneously. Sixth, we ignore the presence of depletion layer around the particle, hence the entire body of the suspending fluid is assumed to obey the same constitutive relation – such a paradigm has been used extensively in the existing literature (Afonso, Alves & Pinho Reference Afonso, Alves and Pinho2009; Choudhary et al. Reference Choudhary, Li, Renganathan, Xuan and Pushpavanam2020; Li & Koch Reference Li and Koch2020; Mahapatra & Bandopadhyay Reference Mahapatra and Bandopadhyay2021). Seventh, the flow field is assumed to be viscosity dominated on account of low Reynolds number.

Finally, it is also important to note the limitations of the Oldroyd-B model used herein (Bird et al. Reference Bird, Armstrong and Hassager1987). Although this constitutive model can be applied to relatively higher strain rates (as compared to ordered fluids) and is successful in capturing stress relaxation and normal stress coefficients to a reasonable extent, it nevertheless fails to account for shear-thinning behaviour and non-zero second normal stress coefficients exhibited by several polymeric fluids (Natu & Ghosh Reference Natu and Ghosh2021). In addition, the applicability of the Oldroyd-B model at extremely high stretching rates, as observed inside the EDL, also becomes questionable. Despite this, the Oldroyd-B model is extremely useful in developing insights into the multi-dimensional flows of complex fluids, because it rightfully captures the essential physics of viscoelastic flows, as is evident from its widespread use in the existing literature (Phan-Thien Reference Phan-Thien1983; Bird et al. Reference Bird, Armstrong and Hassager1987; Tan & Masuoka Reference Tan and Masuoka2005; Aggarwal & Sarkar Reference Aggarwal and Sarkar2008; Mukherjee & Sarkar Reference Mukherjee and Sarkar2011; Turkoz et al. Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018).

3.4.1. The flow field

By substituting (3.5) into (3.1), it is straightforward to derive the $O(\bar {\zeta }_0)$ equations, which read ${\boldsymbol{\mathsf{\tau}}}_1 = 2\boldsymbol{\mathsf{D}}_1$, $-\boldsymbol{\nabla } p_{1} + \nabla ^2\boldsymbol{v}_1 = 0$ and $\boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{v}_1= 0$, subject to $\boldsymbol{v}_1(r = 1,\theta,\phi ) = \boldsymbol \varOmega _1 \times \hat {\boldsymbol{e}}_r + \boldsymbol{U}_1 +\boldsymbol{V}_{HS}^{(1)}$, where the expression for $\boldsymbol{V}_{HS}^{(1)}$ is specified in (3.8a). The force- and torque-free conditions at this order would translate into $\int _{S_{p}} (-p_1\boldsymbol{\mathsf{I}}+{\boldsymbol{\mathsf{\tau}}}_1) \boldsymbol{\cdot } \hat {\boldsymbol{e}}_{r}\,{\rm d}S = 0$ and $\int _{S_{p}} \hat {\boldsymbol{e}}_{r}\times ({\boldsymbol{\mathsf{\tau}}}_1\boldsymbol{\cdot }\hat {\boldsymbol{e}}_{r})\,{\rm d}S = 0$, from which $\boldsymbol{U}_1$ and $\boldsymbol \varOmega _1$ may be determined. The leading-order velocity may be represented using the Lamb's general solution (Pak & Lauga Reference Pak and Lauga2014) and reads

(3.9)\begin{align} \boldsymbol{v}_1 = \sum_{n=1}^{\infty}&\left[\boldsymbol{\nabla}\times\left(\boldsymbol{x} \varLambda_{-(n+1)}\right)+\boldsymbol{\nabla}\varPhi_{-(n+1)}-\frac{n-2}{2n(2n-1)}\, r^2\,\boldsymbol{\nabla}\mathcal{P}_{-(n+1)}\right.\nonumber\\ &\quad \left.{}+\frac{n+1}{n(2n-1)}\,\boldsymbol{x}\mathcal{P}_{-(n+1)}\right], \end{align}

where $\varLambda _{-(n+1)}$, $\varPhi _{-(n+1)}$ and $\mathcal{P}_{-(n+1)}$ are solid harmonics of order $-n-1$, and $r = |\boldsymbol{x}|$. For instance, $\mathcal{P}_{-(n+1)}$ has the expression

(3.10)\begin{equation} \mathcal{P}_{-(n+1)} = r^{{-}n-1}\sum_{l=0}^{n}P^{l}_n(\eta)(A_{ln}\cos l\phi+\tilde{A}_{ln}\sin l\phi), \end{equation}

where $P^{l}_n(\eta )$ is the associated Legendre polynomial of order $n$ and $l$. Expressions for the remaining solid harmonics have been provided in Appendix A, along with expanded expressions for the velocity components. The pressure is given by $p_1 = \sum _{n=1}^{\infty }\mathcal{P}_{-n-1}$. To systematically work out the unknown coefficients $A_{ln}$, $\tilde {A}_{ln}, \ldots$, the velocity on the particle's surface may be written as (Pak & Lauga Reference Pak and Lauga2014)

(3.11)\begin{equation} \boldsymbol{v}_1(r=1,\theta,\phi) = \sum_{l=0}^{\infty}(\boldsymbol{M}_{l}(\theta)\cos l\phi + \tilde{\boldsymbol{M}}_l(\theta)\sin l\phi), \end{equation}

where $\boldsymbol{M}(\theta ) = D_l(\theta )\,\hat {\boldsymbol{e}}_r + E_l(\theta )\, \hat {\boldsymbol{e}}_{\theta } + F_l(\theta )\,\hat {\boldsymbol{e}}_{\phi }$ and $\tilde {\boldsymbol{M}}(\theta ) = \tilde {D}_l(\theta )\,\hat {\boldsymbol{e}}_r + \tilde {E}_l(\theta )\,\hat {\boldsymbol{e}}_{\theta } +$ $\tilde {F}_l(\theta )\,\hat {\boldsymbol{e}}_{\phi }$. Since these velocity components are known from the given boundary condition on the particle surface, one can deduce the coefficients $\boldsymbol{M}$ and $\tilde {\boldsymbol{M}}$ ($l = 0,1,2,\ldots$), from which the unknown coefficients in (3.9) may be jotted down using the orthogonality of the associated Legendre polynomials and following the procedure outlined by Pak & Lauga (Reference Pak and Lauga2014). For instance, the coefficient $A_{ln}$ can be deduced as

(3.12)\begin{equation} A_{ln} = \frac{2n-1}{2}\,\frac{(2n+1)(n-l)!}{(n+1)(n+l)!}\int_{{-}1}^{1}\left[nD_{l} + \frac{\partial}{\partial \eta}(E_{l}\sin\theta) - \frac{l\tilde{F}_{l}}{\sin\theta}\right]P_{n}^{l}(\eta)\,{\rm d}\eta. \end{equation}

The other unknown coefficients may be determined from similar relations, and these have been included in Appendix A. A brief outline of the methodology to compute the functions $\boldsymbol{M}$ and $\tilde {\boldsymbol{M}}$ has been provided in Appendix B.

3.4.2. The leading-order particle velocity

The leading-order particle velocities (translational and angular) may be determined by equating the net hydrodynamic drag and torque on the particle to zero. In terms of the solid harmonics, these force and torque balances may be expressed as (Happel & Brenner Reference Happel and Brenner2012)

(3.13a)$$\begin{gather} \boldsymbol{F} ={-}4{\rm \pi}\,\boldsymbol{\nabla}(r^3\mathcal{P}_{{-}2}) = 0 , \end{gather}$$

(3.13b)$$\begin{gather}\boldsymbol{T} ={-}8{\rm \pi}\,\boldsymbol{\nabla}(r^3\varLambda_{{-}2}) = 0. \end{gather}$$

The above two conditions yield $A_{01} = A_{11} = \tilde {A}_{11} = 0$ for the force balance, and $C_{01} = C_{11} = \tilde {C}_{11} = 0$ for the torque balance – see also (A1b). Hence for a given surface charge distribution, these coefficients may be determined using the procedure described in § 3.4.1, and since they contain $\boldsymbol{U}_1$ and $\boldsymbol \varOmega _1$, those may be evaluated by equating the above coefficients to zero. Analytical examples for specific cases of surface charge distributions have been included in § 3.6.

3.5.1. The generalized reciprocal theorem

The first viscoelastic corrections due to the polymeric stresses contribute at $O(\bar {\zeta }_0^2)$. It is again straightforward to derive the governing equations and the boundary conditions at this order starting from (3.1), and they take the following forms: $\boldsymbol{\nabla } \boldsymbol{\cdot } {\boldsymbol{\mathsf{\Theta}}}_{2} + \boldsymbol{b}_{2} = 0$ and $\boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{v}_2 = 0$, where ${\boldsymbol{\mathsf{\Theta}}}_{2} = -p_2\boldsymbol{\mathsf{I}}+2\boldsymbol{\mathsf{D}}_2$ is the linear part of the stresses at $O(\bar {\zeta }_0^2)$, while $\boldsymbol{b}_2 = \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{\mathsf{\tau}}}_{2}^{exc}$ and ${\boldsymbol{\mathsf{\tau}}}_2^{exc}$ is the excess polymeric stresses at this order. It follows directly from (3.1b) that ${\boldsymbol{\mathsf{\tau}}}_{2}^{exc} = De\,(2\lambda _{2}\boldsymbol{\mathsf{S}}_{2} - \lambda _{1}\boldsymbol{\mathsf{T}}_{2})$ and thus ${\boldsymbol{\mathsf{\tau}}}_2 = 2\boldsymbol{\mathsf{D}}_2+{\boldsymbol{\mathsf{\tau}}}_{2}^{exc}$. Moreover, as a consequence of the regular asymptotic framework, $\boldsymbol{\mathsf{T}}_{2} = 2\boldsymbol{\mathsf{S}}_{2}$ and thus $\boldsymbol{b}_{2} = 2\,De\,(\lambda _2 - \lambda _1)(\boldsymbol{\nabla } \boldsymbol{\cdot } \boldsymbol{\mathsf{S}}_{2})$. On the particle surface, the velocity field satisfies $\boldsymbol{v}_2(r = 1,\theta,\phi ) = \boldsymbol \varOmega _2 \times \hat {\boldsymbol{e}}_r + \boldsymbol{U}_2 +\boldsymbol{V}_{HS}^{(2)}$, where $\boldsymbol{V}_{HS}^{(2)}$ has been specified in (3.8b) and (3.8c). The far-field condition remains unchanged from (3.2b).

The $O(\bar {\zeta }_0^2)$ corrections to the particle velocities may be computed using the generalized reciprocal theorem (Masoud & Stone Reference Masoud and Stone2019), which does not require one to evaluate the complete velocity field at this order. To this end, we consider an auxiliary Stokes flow around a similar spherical particle in a Newtonian medium, wherein the total stress field and the body forces are denoted respectively as $\hat{{\boldsymbol{\mathsf{\Theta}}}}$ ($= -\hat {p} \boldsymbol{\mathsf{I}}+2\hat {\boldsymbol{\mathsf{D}}}$, where $\hat {\boldsymbol{\mathsf{D}}}$ is the auxiliary strain rate tensor) and $\boldsymbol{\hat {b}}$. This auxiliary flow field would satisfy $\boldsymbol{\nabla } \boldsymbol{\cdot } \hat {{\boldsymbol{\mathsf{\Theta}}}} + \hat {\boldsymbol{b}} = 0$ and $\boldsymbol{\nabla } \boldsymbol{\cdot } \hat {\boldsymbol{v}} = 0$, where $\hat {\boldsymbol{v}}$ denotes the velocity field associated with the auxiliary flow. Assuming that $\hat {\boldsymbol{v}}$ vanishes far away from the particle and the stresses of the viscoelastic flow decay at a rate faster than $1/r^{2}$, the generalized reciprocal theorem boils down to (Masoud & Stone Reference Masoud and Stone2019; Li & Koch Reference Li and Koch2020; Ghosh et al. Reference Ghosh, Mukherjee and Chakraborty2021)

(3.14)\begin{equation} \int _{S_p} \hat{\boldsymbol{e}}_r\boldsymbol{\cdot}{\boldsymbol{\mathsf{\Theta}}}_{2} \boldsymbol{\cdot} \hat{\boldsymbol{v}} \,{\rm d}S - \int _{S_p} \hat{\boldsymbol{e}}_r \boldsymbol{\cdot} \hat{{\boldsymbol{\mathsf{\Theta}}}} \boldsymbol{\cdot} \boldsymbol{v}_{2} \,{\rm d}S = \int_{V} \boldsymbol{b}_{2} \boldsymbol{\cdot} \hat{\boldsymbol{v}} \,{\rm d}V - \int_{V} \hat{\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{v}_{2} \,{\rm d}V .\end{equation}

In the above equation, the integrations on the left-hand side are carried out over the particle surface at $r = 1$, and those on the right-hand side are carried out over the entire fluid domain. Hence on $S_{p}$, $\boldsymbol{v}_{2}$ in the second integral on the left-hand side becomes equal to $\boldsymbol{U}_{2} + \boldsymbol \varOmega _{2}\times \hat {\boldsymbol{e}}_{r} + \boldsymbol{V}_{HS}^{(2)}$. The following auxiliary flow fields have been chosen here: (i) to determine the particle's translational velocity ($\boldsymbol{U}_2$), the flow generated by a sphere translating with velocity $\hat {\boldsymbol{e}}$ (unit magnitude, along an arbitrary direction) in a quiescent medium is chosen as the auxiliary flow; (ii) to deduce the particle's rotational velocity ($\boldsymbol \varOmega _2$), we will consider as the auxiliary flow the velocity field of a sphere rotating with angular velocity $\hat {\boldsymbol{e}}$ in a quiescent Newtonian fluid. Therefore, in both the scenarios, $\hat {\boldsymbol{b}} = 0$.

We will first apply the reciprocal theorem to the leading-order problem, which we have already solved in § 3.4, to validate our choices. At this order, $\boldsymbol{b}_{1} = 0$ and hence (3.14) simplifies (replace ${\boldsymbol{\mathsf{\Theta}}}_2$, $\boldsymbol{b}_2$ and $\boldsymbol{v}_2$, respectively, with ${\boldsymbol{\mathsf{\Theta}}}_1$, $\boldsymbol{b}_1 = 0$ and $\boldsymbol{v}_1$) to

(3.15)\begin{equation} \int _{S_p} \hat{\boldsymbol{e}}_r\boldsymbol{\cdot} {\boldsymbol{\mathsf{\Theta}}}_{1} \boldsymbol{\cdot} \hat{\boldsymbol{v}} \,{\rm d}S = \int _{S_p} \hat{\boldsymbol{e}}_r \boldsymbol{\cdot} \hat{{\boldsymbol{\mathsf{\Theta}}}} \boldsymbol{\cdot} \boldsymbol{v}_{1} \,{\rm d}S .\end{equation}

To determine the particle's translation velocity, we chose the auxiliary flow (i) as mentioned above, for which (Pozrikidis & Jankowski Reference Pozrikidis and Jankowski1997) the traction on the particle surface ($r = 1$) is given by $\hat {\boldsymbol{e}}_r\boldsymbol{\cdot }\hat {{\boldsymbol{\mathsf{\Theta}}}} = -\frac {3}{2}\hat {\boldsymbol{e}}$. Now, on $S_{p}$, $\boldsymbol{v}_{1} = \boldsymbol{U}_{1} + \boldsymbol \varOmega _{1}\times \hat {\boldsymbol{e}}_{r} + \boldsymbol{V}_{HS}^{(1)}$ and $\hat {\boldsymbol{v}} = \hat {\boldsymbol{e}}$, whereas the force-free condition requires $\int _{S_p} {\boldsymbol{\mathsf{\Theta}}}_{1}\boldsymbol{\cdot } \hat {\boldsymbol{e}}_r \,{\rm d}S = 0$. As a consequence, (3.15) may be simplified directly to deduce $\boldsymbol{U}_1$ as

(3.16)\begin{equation} \boldsymbol{U}_{1} ={-}\frac{1}{4{\rm \pi}}\int_{Sp}\boldsymbol{V}_{HS}^{(1)}\,{\rm d}S .\end{equation}

On the other hand, the rotational velocity of the particle may be determined by considering the auxiliary flow (ii) stated earlier, so that the traction vector on the particle surface becomes (Pozrikidis & Jankowski Reference Pozrikidis and Jankowski1997) $\hat {\boldsymbol{e}}_{r} \boldsymbol{\cdot } \hat {{\boldsymbol{\mathsf{\Theta}}}} = -3\hat {\boldsymbol{e}}\times \hat {\boldsymbol{e}}_{r}$, while on $S_p$, $\hat {\boldsymbol{v}} = \hat {\boldsymbol{e}}\times \hat {\boldsymbol{e}}_{r}$. The torque-free condition becomes $\int _{S_p} \hat {\boldsymbol{e}}_r\times ({\boldsymbol{\mathsf{\Theta}}}_{1}\boldsymbol{\cdot } \hat {\boldsymbol{e}}_r) \,{\rm d}S = 0$, hence (3.15) may be simplified to determine the expression for the particle's rotational velocity as

(3.17)\begin{equation} \boldsymbol{\varOmega}_{1} ={-}\frac{3}{8{\rm \pi}}\int_{0}^{2{\rm \pi}}{\rm d}\phi\int_{{-}1}^{1}{\rm d}\eta\,(\hat{\boldsymbol{e}}_{r}\times\boldsymbol{V}_{HS}^{(1)}). \end{equation}

This is identical to the expression derived by Ghosh et al. (Reference Ghosh, Mukherjee and Chakraborty2021). Remarkably, application of the reciprocal theorem demonstrates that the leading-order (Newtonian) contribution to the particle's translational velocity is simply the surface average of the Smoluchowski slip velocity, while the angular velocity is simply the average moment of the slip velocity around the particle centre.

3.5.2. Viscoelastic corrections using the reciprocal theorem

We now apply (3.14) to determine $\boldsymbol{U}_2$ and $\boldsymbol \varOmega _2$. To this end, we recall that $\boldsymbol{b}_{2} = 2\,De\,(\lambda _2 - \lambda _1)(\boldsymbol{\nabla } \boldsymbol{\cdot } \boldsymbol{\mathsf{S}}_{2})$ (because $\boldsymbol{\mathsf{T}}_2 = 2\boldsymbol{\mathsf{S}}_2$, a consequence of linearization around the Newtonian solution using regular perturbation). The force-free condition at $O(\bar {\zeta }_0^2)$ reads (see (3.4a,b)) $\int _{S_p} ({\boldsymbol{\mathsf{\Theta}}} _2+{\boldsymbol{\mathsf{\tau}}}_2^{exc})\boldsymbol{\cdot }\hat {\boldsymbol{e}}_r\,{\rm d}S = 0$, and upon enforcing the expression for ${\boldsymbol{\mathsf{\tau}}}_2^{exc}$, it simplifies to

(3.18)\begin{equation} \int_{S_{p}}{\boldsymbol{\mathsf{\Theta}}}_{2} \boldsymbol{\cdot} \hat{\boldsymbol{e}}_{r}\,{\rm d}S ={-}\int_{S_{p}}2\,De\,(\lambda_2 - \lambda_1)\boldsymbol{\mathsf{S}}_{2}\boldsymbol{\cdot}\hat{\boldsymbol{e}}_{r}\,{\rm d}S. \end{equation}

To deduce $\boldsymbol{U}_2$, we use the same auxiliary flow as outlined in § 3.5.1 (see the discussion prior to (3.16)), along with (3.18) and the boundary condition for $\boldsymbol{v}_2$ on $S_p$ to simplify (3.14), which yields

(3.19)\begin{align} \boldsymbol{U}_{2}\boldsymbol{\cdot}\hat{\boldsymbol{e}} ={-}\frac{1}{4{\rm \pi} }\,\hat{\boldsymbol{e}} \boldsymbol{\cdot}\int _{S_p} \boldsymbol{V}_{HS}^{(2)}\,{\rm d}S \!+ \frac{De}{3{\rm \pi}}\,(\lambda_2 - \lambda_1)\left\{ \hat{\boldsymbol{e}}\boldsymbol{\cdot}\int _{S_p} \boldsymbol{\mathsf{S}}_{2} \boldsymbol{\cdot} \hat{\boldsymbol{e}}_{r}\,{\rm d}S +\!\! \int _{V}\hat{\boldsymbol{v}}\boldsymbol{\cdot}(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\mathsf{S}}_{2})\,{\rm d}V\right\}\!. \end{align}

The terms inside the braces on the right-hand side may be simplified further by using Gauss’ divergence theorem and noting that $\hat {\boldsymbol{v}} = \hat {\boldsymbol{e}}$ on $S_p$, to arrive at

(3.20)\begin{equation} \boldsymbol{U}_{2}\boldsymbol{\cdot}\hat{\boldsymbol{e}} ={-}\frac{1}{4{\rm \pi} }\hat{\boldsymbol{e}} \boldsymbol{\cdot}\int _{S_p} \boldsymbol{V}_{HS}^{(2)}\,{\rm d}S - \frac{De}{3{\rm \pi}}\,(\lambda_2 - \lambda_1)\int_{V} \boldsymbol{\mathsf{S}}_2 : \boldsymbol{\nabla}\hat{\boldsymbol{v}}\,{\rm d}V. \end{equation}

The three components of $\boldsymbol{U}_2$ may be evaluated by choosing suitably $\hat {\boldsymbol{e}}$ – the direction of the sphere's velocity in the auxiliary flow field. The above expression for $\boldsymbol{U}_2$ is very similar to that reported by Li & Koch (Reference Li and Koch2020), albeit for axisymmetric flows. For completeness, it is important to note that the velocity field in the auxiliary flow is given by (Pozrikidis & Jankowski Reference Pozrikidis and Jankowski1997)

(3.21)\begin{equation} \hat{\boldsymbol{v}} = \frac{3\hat{\boldsymbol{e}}}{4r}\left(1+\frac{1}{3r^2}\right)+\frac{3(\hat{\boldsymbol{e}}\boldsymbol{\cdot}\boldsymbol{x})\boldsymbol{x}}{4r^3}\left(1-\frac{1}{r^2}\right) . \end{equation}

We note from (3.20) that the $O(\bar {\zeta }_0^2)$ migration velocity of the particle has two major elements: (a) a contribution from the surface average of the $O(\bar {\zeta }_0^2)$ slip velocity, which is very similar to the leading-order particle velocity; and (b) a second contribution from the excess polymeric stresses in the bulk of the flow, as designated by the final integral term. Evaluation of the contributions from (b) requires complete knowledge of $\hat {\boldsymbol{v}}$ as well as $\boldsymbol{v}_1$ – the leading-order flow field.

In order to deduce the particle's rotational velocity, we will use the auxiliary flow (ii) as outlined in § 3.5.1, for which the velocity field is given by

(3.22)\begin{equation} \hat{\boldsymbol{v}} = \boldsymbol{\nabla}\times\left(\frac{\hat{\boldsymbol{e}}\boldsymbol{\cdot}\boldsymbol{x}\boldsymbol{x}}{r^3}\right). \end{equation}

The torque-free condition at this order requires $\int _{S_p} \hat {\boldsymbol{e}}_r\times \{({\boldsymbol{\mathsf{\Theta}}} _2+{\boldsymbol{\mathsf{\tau}}}_2^{exc})\boldsymbol{\cdot }\hat {\boldsymbol{e}}_r\}\,{\rm d}S = 0$, which finally yields (Ghosh et al. Reference Ghosh, Mukherjee and Chakraborty2021)

(3.23)\begin{equation} \int_{S_{p}}\hat{\boldsymbol{e}}_{r}\times\left({\boldsymbol{\mathsf{\Theta}}}_{2}\boldsymbol{\cdot}\hat{\boldsymbol{e}}_{r}\right){\rm d}S = -2\,De\,(\lambda_2 - \lambda_1)\int_{S_{p}}\hat{\boldsymbol{e}}_{r}\times\left(\boldsymbol{\mathsf{S}}_{2}\boldsymbol{\cdot}\hat{\boldsymbol{e}}_{r}\right){\rm d}S. \end{equation}

Using (3.23) and the boundary conditions for $\hat {\boldsymbol{v}}$ and $\boldsymbol{v}_2$ on $S_p$, and following the same procedure as that for $\boldsymbol \varOmega _1$, (3.14) yields (Ghosh et al. Reference Ghosh, Mukherjee and Chakraborty2021) at $O(\bar {\zeta }_0^2)$ that

(3.24)\begin{align} \boldsymbol{\varOmega}_{2}\boldsymbol{\cdot}\hat{\boldsymbol{e}} & ={-}\frac{3}{8{\rm \pi}}\,\hat{\boldsymbol{e}}\boldsymbol{\cdot}\int _{S_p}\left(\hat{\boldsymbol{e}}_{r}\times \boldsymbol{V}_{HS}^{(2)} \right) {\rm d}S + \frac{3\,De\,(\lambda_2 - \lambda_1)}{4{\rm \pi}}\left\{\hat{\boldsymbol{e}}\boldsymbol{\cdot}\int _{S_p} \hat{\boldsymbol{e}}_r \times(\boldsymbol{\mathsf{S}}_{2} \boldsymbol{\cdot} \hat{\boldsymbol{e}}_{r})\,{\rm d}S \right. \nonumber\\ &\quad \left. {}+\int _{V} \hat{\boldsymbol{v}}\boldsymbol{\cdot}(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\mathsf{S}}_{2}) \,{\rm d}V \right\}. \end{align}

The above equation may also be simplified further using the divergence theorem and the fact that on $S_p$, $\hat {\boldsymbol{v}} = \hat {\boldsymbol{e}}\times \hat {\boldsymbol{e}}_r$, to write

(3.25)\begin{equation} \boldsymbol{\varOmega}_{2}\boldsymbol{\cdot}\hat{\boldsymbol{e}} ={-}\frac{3}{8{\rm \pi}}\,\hat{\boldsymbol{e}}\boldsymbol{\cdot}\int _{S_p}\left(\hat{\boldsymbol{e}}_{r}\times \boldsymbol{V}_{HS}^{(2)} \right) {\rm d}S - \frac{3\,De\,(\lambda_2 - \lambda_1)}{4{\rm \pi}} \int_{V} \boldsymbol{\mathsf{S}}_2 : \boldsymbol{\nabla}\hat{\boldsymbol{v}}\,{\rm d}V .\end{equation}

The different components of $\boldsymbol \varOmega _2$ may be evaluated by choosing $\hat {\boldsymbol{e}}$ appropriately. Akin to the translational velocity, the angular velocity at $O(\bar {\zeta }_0^2)$ also has contributions from two analogous elements: (a) the moment of the $O(\bar {\zeta }_0^2)$ slip velocity about the origin; and (b) the excess polymeric stresses in the bulk, as expressed through the second integral in (3.25).

It is important to note from (3.20) and (3.25) that all the terms are linearly proportional to $De\,(\lambda _2-\lambda _1)$ (see (3.8b) and (3.8c) for $\boldsymbol{V}_{HS}^{(2)}$) and therefore the $O(\bar {\zeta }_0^2)$ corrections to the particle velocities originate exclusively from the viscoelasticity of the suspending medium. The leading-order flow stretches the polymers inside the EDL, resulting in excess polymeric stresses therein, and causes changes in the slip velocity at the edge of the EDL and ultimately in the particle's migration velocity. At the same time, the excess polymeric stresses in the bulk also alter the net hydrodynamic drag on the particle, thus leading to a change in its velocity. Although these two effects are not mutually exclusive, they nevertheless may aid or oppose each other to govern the net viscoelastic correction to the particle's motion. Conversely, in many other cases, one or more of these effects may vanish identically, which may cause either the $O(\bar {\zeta }_0^2)$ slip velocity or the bulk polymeric stresses to dominate the viscoelastic corrections.

For completeness, one may note that $\boldsymbol{\mathsf{S}}_2 = {{\rm D}(\boldsymbol{\mathsf{D}}_1)}/{{\rm D}t} - {(\boldsymbol{\nabla } \boldsymbol{v}_1)}^{\rm T} \boldsymbol{\cdot } \boldsymbol{\mathsf{D}}_1 - \boldsymbol{\mathsf{D}}_1 \boldsymbol{\cdot } (\boldsymbol{\nabla } \boldsymbol{v}_1)$, thus this may be computed once $\boldsymbol{v}_1$ is known. Hence once the instantaneous surface charge distribution $\breve {\zeta }(\theta,\phi )$ is known, the leading-order velocity field $\boldsymbol{v}_1$ may be determined using the analysis of § 3.4, and $\boldsymbol{V}_{HS}^{(2)}$ will be known from (3.8b) and (3.8c), which would enable us to evaluate $\boldsymbol{U}_2$ and $\boldsymbol \varOmega _2$ using (3.20) and (3.25). Then, using either (3.16) and (3.17) or (3.13), the total particle velocity $\boldsymbol{U}$ and $\boldsymbol \varOmega$ up to $O(\bar {\zeta }_0^2)$ can be found as

(3.26a)$$\begin{gather} \boldsymbol{U} = \bar{\zeta}_0\boldsymbol{U}_1 + \bar{\zeta}_0^{2}\boldsymbol{U}_2 + O(\bar{\zeta}_0^3) , \end{gather}$$

(3.26b)$$\begin{gather}\boldsymbol\varOmega = \bar{\zeta}_0\boldsymbol\varOmega_{1} + \bar{\zeta}_0^{2}\boldsymbol\varOmega_{2} + O(\bar{\zeta}_0^3) . \end{gather}$$

We reiterate that in the above, the leading-order velocity is purely Newtonian and the $O(\bar {\zeta }_0^2)$ corrections are purely viscoelastic in nature, hence the viscoelasticity always remains subdominant here.

3.6.1. Special case I: axisymmetric surface charge ($a_{0},a_{10} \neq 0$, $a_{11},\tilde {a}_{11} = 0$)

For an axisymmetric surface charge distribution of the form $\breve {\zeta }(\theta ) = a_0+a_1\eta$, we deduce from (3.27) and (3.29) that $\boldsymbol \varOmega = 0$ and $\boldsymbol{U} = (\bar {\zeta }_0 a_0\beta -(3/5)\bar {\zeta }_0^2\beta ^2\,De\,(\lambda _1-\lambda _2)a_0 a_{10})\hat {\boldsymbol{e}}_z$. This is identical to the velocity reported previously by Ghosh et al. (Reference Ghosh, Mukherjee and Chakraborty2021). Evidently, an axisymmetrically charged particle would translate only along $z$ without undergoing any rotation.

3.6.2. Special case II: particle without any net charge ($a_{0} = 0$ and $a_{10},a_{11}, \tilde {a}_{11} \neq 0$)

When $a_0 = 0$, the particle carries zero net charge and $\breve {\zeta }(\theta,\phi ) = a_{10}\eta + \sqrt {1-\eta ^{2}}(a_{11}\cos \phi + \tilde {a}_{11}\sin \phi )$. The particle velocities in (3.27) and (3.29) thus appear as: $\boldsymbol{U} = \frac {3}{40} \bar {\zeta }_0^2\{De\,\beta ^2(\lambda _1-\lambda _2)(a_{10}\tilde {a}_{11}\hat {\boldsymbol{e}}_x-a_{10} a_{11} \hat {\boldsymbol{e}}_y)\}$ and $\boldsymbol \varOmega = ({3\beta }/{4})(\tilde {a}_{11}\hat {\boldsymbol{e}}_x-a_{11} \hat {\boldsymbol{e}}_y)\bar {\zeta }_0 + \{\frac {9}{16}(\lambda _1-\lambda _2)\beta ^2\,De\, (a_{10}\tilde {a}_{11}\hat {\boldsymbol{e}}_x-a_{10} a_{11} \hat {\boldsymbol{e}}_y)\}\bar {\zeta }_0^2$. Remarkably, despite carrying no net charge, $\boldsymbol{U}\neq 0$, although $\boldsymbol{U}_1 = 0$ as expected, and hence the particle will migrate when suspended in a viscoelastic fluid. Further notice that $\boldsymbol{U}_2$ and $\boldsymbol \varOmega _2$ are identical, except for their numeric pre-factor. This is in stark contrast to our intuition, which perhaps dictates that a particle carrying zero net charge should never possess a translational velocity under a uniform external electric field. The reason may be attributed to the breaking of fore–aft symmetry in viscoelastic fluids for non-uniformly charged particles, as shown previously by Ghosh et al. (Reference Ghosh, Mukherjee and Chakraborty2021). The leading-order rotation of the particle drags the polymers around it, thus stretching them in the process, which leads to excess normal and shear stresses on the particle surface. This, coupled with the tendency of the suspending fluid to break fore–aft symmetry, forces the particle to move in a certain direction. Indeed, it has been identified recently (Kroo et al. Reference Kroo, Binagia, Eckman, Prakash and Shaqfeh2022) that a particle lacking fore–aft symmetry may propel itself forwards by undergoing rotation, because of the viscoelasticity of the surrounding medium.

3.7.1. Impact of the shear-thinning nature of the suspending medium

It is well known (Bird et al. Reference Bird, Armstrong and Hassager1987) that most polymeric liquids and aqueous polymer solutions exhibit shear-thinning behaviour to some extent, although the Oldroyd-B model considered in this work is unable to capture such behaviour. It is thus important to highlight how such effects may be potentially incorporated within the current framework. To this end, we note that viscoelastic models such as the Phan-Thien–Tanner model, the finitely extensible nonlinear elastic (FENE) model, and Giesekus constitutive relations have been shown to be capable of capturing the shear-thinning nature of various polymeric liquids, and all of these have a number of similarities with the Oldroyd-B model. As a result, it may be inferred that much of the analytical detail outlined in the present work will remain applicable, with a few notable exceptions.

The general methodology of dividing the flow domain into two distinct regions (see §§ 3.1 and 3.3), namely the EDL (inner layer) and the electroneutral bulk, will still be applicable for shear-thinning constitutive models. Moreover, the orders of magnitude of the various stress components and the overall methodology followed to describe the flow in the outer layer (bulk layer) will also remain unchanged. This entails a leading-order (in $\bar {\zeta }_0$) Newtonian velocity field (§ 3.4), followed by viscoelastic corrections computed using the reciprocal theorem (§ 3.5). However, the scaling patterns for the stresses inside the EDL may change, depending on the specific constitutive relation. In particular, the viscoelastic corrections to the Smoluchowski slip velocity outlined in § 3.3 stem from the normal stress components ($\tau _{\theta \theta }$, $\tau _{\phi \phi }$, etc.) scaling as $O(\delta ^{-2})$ inside the EDL, which is distinct from Newtonian fluids. This scaling pattern as well as those for the shear stress components will be sensitive to the choice of the constitutive relation, and as a consequence, the viscoelastic corrections to the slip velocity ($\boldsymbol{V}_{HS}^{(2)}$) are expected to get altered. This, in turn, will change the electrophoretic mobility of the particle by altering the $O(\bar {\zeta }_0^2)$ flow in the outer region. Several other aspects of the framework, on the other hand, will remain intact. For instance, the relations such as $\tau _{\theta \theta }\sim O(\tau _{r\theta }^2)$ inside the EDL, as is the case here (Ghosh et al. Reference Ghosh, Mukherjee and Chakraborty2021), will remain valid even for the FENE model or other constitutive models of similar type. In addition, for sufficiently low surface potential ($\bar {\zeta }_0 \ll O(1)$), the leading-order velocity within the EDL will be Newtonian, while viscoelastic corrections will arise only at $O(\bar {\zeta }_0^2)$ (Li & Koch Reference Li and Koch2020; Ghosh et al. Reference Ghosh, Mukherjee and Chakraborty2021).

Recently, Ghosh et al. (Reference Ghosh, Mukherjee and Chakraborty2021) have established using numerical simulations that below a critical surface charge density, the predictions from the Oldroyd-B model remain quantitatively accurate when compared to those from shear-thinning viscoelastic relations such as the FENE-P model. Nevertheless, the Oldroyd-B model used here does capture accurately the essential physics governing the flow around the particle, and its general usefulness in describing viscoelastic flows remains undoubted, as is evident from its popularity in the existing literature (Aggarwal & Sarkar Reference Aggarwal and Sarkar2008; Mukherjee & Sarkar Reference Mukherjee and Sarkar2011; Turkoz et al. Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018); see also § 2.2.2.

We reiterate that shear thinning does play a crucial role in governing electrophoretic motion of particles, as shown by recent experimental studies (Li & Xuan Reference Li and Xuan2018; Malekanfard et al. Reference Malekanfard, Ko, Li, Bulloch, Baldwin, Wang, Fu and Xuan2019). For instance, Malekanfard et al. (Reference Malekanfard, Ko, Li, Bulloch, Baldwin, Wang, Fu and Xuan2019) have shown experimentally that particles tend to migrate towards the channel walls in electro-osmotic flows of shear-thinning fluids, which is in stark contrast to Newtonian media, where the particles generally migrate towards the centre line of the channel due to the wall-induced electrical lift and inertia. In general, it is expected that shear-thinning effects will result in an enhancement of the electrophoretic velocity of the particle (Lee, Chen & Hsu Reference Lee, Chen and Hsu2005; Malekanfard et al. Reference Malekanfard, Ko, Li, Bulloch, Baldwin, Wang, Fu and Xuan2019), and these increments can become very significant for particles carrying large surface charge densities or in the presence of very thin EDLs, which naturally facilitate a reduction in the effective viscosity because of increased strain rates.

3.7.2. Impact of bounding walls on particle trajectories

Electrophoresis is always carried out within confined channels in laboratory and commercial applications alike, thus the presence of the bounding walls may influence a particle's motion significantly. Although the effects of bounding surfaces have not been accounted for in our work, it is nevertheless important to shed light on how such effects may be incorporated within the present framework.

Usually, for creeping motion of Newtonian fluids, the presence of walls is accounted for by using the method of reflections (Ho & Leal Reference Ho and Leal1976; Happel & Brenner Reference Happel and Brenner2012). In the present scenario, using reflection, the leading-order flow field ($\boldsymbol{v}_1$ at $O(\bar {\zeta }_0)$) may be written as $\boldsymbol{v}_{1} = \boldsymbol{v}_{1}^{(1)} + \boldsymbol{v}_{1}^{(2)} + \boldsymbol{v}_{1}^{(3)} + \cdots$, where $\boldsymbol{v}_{1}^{(1)}$ is the solution for flow around an unbounded particle (§ 3.4), $\boldsymbol{v}_{1}^{(2)}$ is the first reflection of the flow field $\boldsymbol{v}_{1}^{(1)}$ from the wall, $\boldsymbol{v}_{1}^{(3)}$ is the reflection of $\boldsymbol{v}_{1}^{(2)}$ from the particle, and so on. For instance, $\boldsymbol{v}_{1}^{(2)}$ will satisfy $-\boldsymbol{\nabla } p_{1}^{(2)} + \nabla ^2\boldsymbol{v}_{1}^{(2)} = 0$ and $\boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{v}_1^{(2)}= 0$, subject to $\boldsymbol{v}_{1}^{(2)} = - \boldsymbol{v}_{1}^{(1)}$ on the walls. Accordingly, the leading-order ($O(\bar {\zeta }_0)$, Newtonian) translational and angular velocities of the particle may be estimated as $\boldsymbol{U}_{1} = \boldsymbol{U}_{1}^{(1)} + \boldsymbol{U}_{1}^{(3)} + \cdots$ and $\boldsymbol{\varOmega }_{1} = \boldsymbol{\varOmega }_{1}^{(1)} + \boldsymbol{\varOmega }_{1}^{(3)} + \cdots$, where $\boldsymbol{U}_{1}^{(1)}$, $\boldsymbol{\varOmega }_{1}^{(1)}$ and $\boldsymbol{v}_{1}^{(1)}$ remain the same as outlined in § 3.4, while $\boldsymbol{U}_{1}^{(3)}$ and $\boldsymbol{\varOmega }_{1}^{(3)}$ represent the particle's translational and rotational velocities in response to the first reflected flow field $\boldsymbol{v}_1^{(2)}$. However, for particles carrying an arbitrary surface charge that intrinsically requires a large number of spherical harmonics, accounting for reflections may become algebraically extremely tedious. The viscoelastic corrections to the electrophoretic mobility ($\boldsymbol{U}_{2},\boldsymbol{\varOmega }_{2}$) may still be computed using the reciprocal theorem as described in § 3.5.2, although now the leading-order velocity field $\boldsymbol{v}_1$ will also contain repeated reflections from the walls and the particle; the application of the reciprocal theorem, however, would in principle remain identical regardless of the presence of the walls (Ho & Leal Reference Ho and Leal1976).

When the particle is sufficiently close to the wall, it may not remain force- and torque-free, as enforced in (3.13a) and (3.13b). In fact, Yariv (Reference Yariv2006) has shown that close to a flat bounding surface, the presence of the particle will make the electric field non-uniform, which in turn may result in a non-zero electrical force (and torque), akin to what is observed in dielectrophoresis. In addition, the bounding surfaces themselves may also be charged, which, in the presence of an imposed external electric field, will drive an electro-osmotic flow of their own. Although this background electro-osmotic flow will not alter the slip velocity on the particle and may be accounted for by an appropriate far-field condition in (3.2b), it may, however, aid or oppose the electrophoretic motion of the particle itself, thus changing its trajectory. It is important to note that externally imposed flows are often used in conjunction with electrophoretic transport in laboratory-scale applications such as ‘free-flow electrophoresis’, which can be very effective tools for separating particles with varying electrophoretic mobilities (Westermeier Reference Westermeier2016).

In general, the presence of bounding surfaces may slow down particles in Newtonian media, wherein migration results only from non-zero inertial effects (Ho & Leal Reference Ho and Leal1974). On the other hand, particles suspended in an imposed flow of viscoelastic fluids do migrate even without the presence of inertia. A number of numerical (Ho & Leal Reference Ho and Leal1976; D'Avino et al. Reference D'Avino, Tuccillo, Maffettone, Greco and Hulsen2010; Li, McKinley & Ardekani Reference Li, McKinley and Ardekani2015; Yuan et al. Reference Yuan, Zhao, Yan, Tang, Alici, Zhang and Li2018) and experimental (Li & Xuan Reference Li and Xuan2018; Malekanfard et al. Reference Malekanfard, Ko, Li, Bulloch, Baldwin, Wang, Fu and Xuan2019) studies have been carried out in this regard, which explore primarily the migration of electroneutral (uncharged) particles in confined flows of viscoelastic fluids driven by pressure gradients. These studies generally reveal that inertial and elastic effects together govern the migration in viscoelastic fluids. When inertia is negligible, particle migration depends on the fluid constitution. For instance, in Oldroyd-B or second-order fluids (non-shear-thinning), numerical and theoretical estimates (Ho & Leal Reference Ho and Leal1976; Li et al. Reference Li, McKinley and Ardekani2015) suggest that particles tend to migrate towards the channel centreline. However, when the suspending medium becomes shear-thinning, secondary flows may develop that, depending on the initial position of the particle, may cause it to migrate either towards the wall or towards the centreline of the channel. When inertial effects become important, competition between the elastic and inertial forces governs particle trajectories, and the particles may settle at a location where the two forces balance each other (Li et al. Reference Li, McKinley and Ardekani2015). In the present paper also (see § 5), we establish that the elastic nature of the fluid tends to cause particles carrying non-uniform surface charge to migrate normal to the direction of the applied electric field.