## 1. Introduction

Electrophoresis is characterized as the motion of charged particles under the action of externally imposed electric fields in a liquid medium (Levich Reference Levich1962; Saville Reference Saville1977; Anderson Reference Anderson1985; Ohshima Reference Ohshima1996; Yariv & Brenner Reference Yariv and Brenner2003; Ohshima Reference Ohshima2006; Khair & Squires Reference Khair and Squires2009; Schnitzer *et al.* Reference Schnitzer, Zeyde, Yavneh and Yariv2013; Schnitzer & Yariv Reference Schnitzer and Yariv2014; Goswami *et al.* Reference Goswami, Dhar, Ghosh and Chakraborty2017), which may or may not contain ions in the form of electrolytes (Saville Reference Saville1977; Ohshima Reference Ohshima2006). This phenomenon has been a prominent area of research over decades owing to its wide gamut of applications such as colloid science (Ohshima Reference Ohshima2006; Khair, Posluszny & Walker Reference Khair, Posluszny and Walker2012) and separation of particles (Ghosal Reference Ghosal2006), DNA analysis (Khair & Squires Reference Khair and Squires2009), gel electrophoresis (Babnigg & Giometti Reference Babnigg and Giometti2004), particle deposition (Besra & Liu Reference Besra and Liu2007) among many others. In many instances, electrophoretic motion occurs in complex media (Ramautar, Demirci & de Jong Reference Ramautar, Demirci and de Jong2006), such as bio-fluids (Babnigg & Giometti Reference Babnigg and Giometti2004; Kremser, Blaas & Kenndler Reference Kremser, Blaas and Kenndler2004), polymeric solutions (Li & Koch Reference Li and Koch2020; Barron, Sunada & Blanch Reference Barron, Sunada and Blanch1995) etc., whose constitutive behaviours show strong deviations from the Newtonian paradigm. These facets have been progressively becoming important in recent times (Berli Reference Berli2010; Bandopadhyay & Chakraborty Reference Bandopadhyay and Chakraborty2012*b*; Bandopadhyay, Ghosh & Chakraborty Reference Bandopadhyay, Ghosh and Chakraborty2013; Zhao & Yang Reference Zhao and Yang2013; Ghosh & Chakraborty Reference Ghosh and Chakraborty2015; Ghosh, Chaudhury & Chakraborty Reference Ghosh, Chaudhury and Chakraborty2016), because of their key roles in medical diagnostics (Madou *et al.* Reference Madou, Lee, Daunert, Lai and Shih2001; Groisman, Enzelberger & Quake Reference Groisman, Enzelberger and Quake2003), particle focusing (Lu *et al.* Reference Lu, DuBose, Joo, Qian and Xuan2015), efficient micromixing (Lam *et al.* Reference Lam, Gan, Nguyen and Lie2009) to name a few. Some such specific applications include: gel electrophoresis for proteome sequencing (Babnigg & Giometti Reference Babnigg and Giometti2004), capillary electrophoresis for separation and purification of biological substances (Karger, Cohen & Guttman Reference Karger, Cohen and Guttman1989; Kremser *et al.* Reference Kremser, Blaas and Kenndler2004; Ramautar *et al.* Reference Ramautar, Demirci and de Jong2006), measurements of large aggregates of viruses, bacteria and eukaryotic cells (Kremser *et al.* Reference Kremser, Blaas and Kenndler2004) among others.

It is well established that fluid media of emerging interest in electrically mediated transport of biological entities commonly exhibit viscoelastic behaviour (Skalak, Ozkaya & Skalak Reference Skalak, Ozkaya and Skalak1989; Brust *et al.* Reference Brust, Schaefer, Doerr, Pan, Garcia, Arratia and Wagner2013). However, a majority of studies on electrophoresis (Saville Reference Saville1977; Baygents & Saville Reference Baygents and Saville1991; Ohshima Reference Ohshima2006; Khair & Squires Reference Khair and Squires2009) tend to focus on particle motion in Newtonian fluids, while only a handful of investigations (Hsu, Yeh & Ku Reference Hsu, Yeh and Ku2006; Hsu & Yeh Reference Hsu and Yeh2007; Khair *et al.* Reference Khair, Posluszny and Walker2012; Posluszny Reference Posluszny2014; Li & Koch Reference Li and Koch2020) have addressed the phenomenon in a non-Newtonian medium. Among those, many of the investigations concern Carreau-type constitutive models (Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1987; Khair *et al.* Reference Khair, Posluszny and Walker2012), which are designed to capture the shear dependent viscosity exhibited by many polymeric liquids, but are unable to describe many other important characteristics (such as the ‘normal stress effects’) exhibited by such liquids (Bird *et al.* Reference Bird, Armstrong and Hassager1987; Ghosh *et al.* Reference Ghosh, Chaudhury and Chakraborty2016). On the other hand, viscoelastic constitutive models, such as the retarded motion relations (Chan & Leal Reference Chan and Leal1979) as well as the differential constitutive relations (Bird *et al.* Reference Bird, Armstrong and Hassager1987; Mukherjee & Sarkar Reference Mukherjee and Sarkar2010; Ghosh & Chakraborty Reference Ghosh and Chakraborty2015), hold certain favourable characteristics in capturing many of the essential rheological features of various polymeric (Bird *et al.* Reference Bird, Armstrong and Hassager1987; Barron *et al.* Reference Barron, Sunada and Blanch1995; Li & Koch Reference Li and Koch2020) as well as biological fluids (Yeleswarapu *et al.* Reference Yeleswarapu, Kameneva, Rajagopal and Antaki1998; Brust *et al.* Reference Brust, Schaefer, Doerr, Pan, Garcia, Arratia and Wagner2013). As a result, numerous types of viscoelastic flows have been widely studied (for instance, Vamerzani, Norouzi & Firoozabadi Reference Vamerzani, Norouzi and Firoozabadi2014; Turkoz *et al.* Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018) over the years, although rigorous studies on electrophoretic motion through viscoelastic medium are extremely scarce. Lu *et al.* (Reference Lu, Patel, Zhang, Woo Joo, Qian, Ogale and Xuan2014, Reference Lu, DuBose, Joo, Qian and Xuan2015) have carried out experiments on capillary electrophoresis in viscoelastic fluids, wherein intriguing streamwise particle oscillations were reported that are otherwise not witnessed in Newtonian media. Very recently, Li & Koch (Reference Li and Koch2020) have theoretically investigated electrophoresis in dilute polymer solutions for a uniformly charged particle and predicted a reduction in the electrophoretic mobility because of the polymeric stresses inside the electrical double layer (EDL).

Another important feature of most of the reported theoretical studies on electrophoresis is that they mainly focus on spherical particles with uniform surface charge density (Saville Reference Saville1977; Ohshima Reference Ohshima2006; Khair & Squires Reference Khair and Squires2009; Li & Koch Reference Li and Koch2020). However, in many cases, the moving particles might have non-uniformly distributed charges on their surfaces, either naturally induced or externally imposed. As Anderson (Reference Anderson1985) points out, phenomena like heterocoagulation and crystallinity may lead to non-uniform charge density on a particle's surface. In addition, particles like bacterial cells and DNA (Amro *et al.* Reference Amro, Kotra, Wadu-Mesthrige, Bulychev, Mobashery and Liu2000; Velegol Reference Velegol2002), various inorganic particles (Van Riemsdijk *et al.* Reference Van Riemsdijk, Bolt, Koopal and Blaakmeer1986; Velegol Reference Velegol2002) and particles with adsorbed surfactants (Fleming, Wanless & Biggs Reference Fleming, Wanless and Biggs1999) might naturally exhibit non-uniform surface charge density. At the same time, ‘Janus particles’ (Molotilin, Lobaskin & Vinogradova Reference Molotilin, Lobaskin and Vinogradova2016) might be artificially engineered to have varying surface properties, which may lead to non-uniform surface charges among other features. In view of such prevalence of non-uniformly charged particles, several studies (Anderson Reference Anderson1985; Yoon Reference Yoon1991; Velegol Reference Velegol2002) have been carried out to derive expressions for their electrophoretic mobility, albeit exclusively in Newtonian fluids. Several researchers have also probed the mobility of non-spherical particles (Fair & Anderson Reference Fair and Anderson1989; Solomentsev & Anderson Reference Solomentsev and Anderson1994; Yariv Reference Yariv2005). While these studies clearly predict that the non-uniformity in the surface charge alters the electrophoretic mobility, the underlying implications in non-Newtonian fluid medium have not yet been probed.

The interactions between the viscoelastic polymeric stresses and the Maxwell stresses within the EDL would fundamentally alter the flow dynamics within the same, and hence hold the potential of resulting in substantial alterations in fluid motion around the particle. Considering that perspective and realizing an effective compromise between rheological complexity and analytical tractability, here, we analyse electrokinetics around a non-uniformly charged spherical particle in a viscoelastic medium, whose constitutive behaviour is given by the Oldroyd-B model. This constitutive relation is chosen over the ordered-fluid models, because of its applicability to flows with comparatively higher strain rates (Bird *et al.* Reference Bird, Armstrong and Hassager1987).

We confine our attention to the thin EDL limit (thin EDL) (Ajdari Reference Ajdari1995; Mandal *et al.* Reference Mandal, Ghosh, Bandopadhyay and Chakraborty2015) and first formulate a generalized framework to evaluate the modified Smoluchowski slip velocity around the particle with arbitrary non-uniform surface charge. In order to work with a closed-formed theoretical framework, we assume the particle to be weakly charged; however, the surface charge is otherwise arbitrary and non-uniform. We illustrate that, for weak surface charge, the viscoelastic effects become subdominant in the leading-order asymptotes, which essentially results in a small effective Deborah number ($De_{eff}$, defined later).

As simple applications of the modified Smoluchowski slip thus derived, we consider two specific case studies. First, we explore the electrophoretic translation of a non-uniformly but axisymmetrically charged particle. Second, we probe the pure electrophoretic rotation of a non-uniformly and non-axisymmetrically charged particle in a viscoelastic medium. The analysis is carried out by using a combination of singular and regular asymptotic expansions. Noting that the applicability of Oldroyd-B model becomes questionable at relatively moderate to large dimensionless relaxation times, as expressed in terms of Deborah numbers (defined later), we further report numerical solutions of the physical problem by employing an illustrative nonlinear viscoelastic model, namely, the finitely extensible nonlinear elastic (FENE-P) model, in appropriate limiting scenarios. Subsequent comparisons reveal that the predictions using the Oldroyd-B model agree reasonably well with the FENE-P model as long as the surface charge is weak and the maximum allowed polymer lengths are large. One of the main goals of the case studies stated above is to bring out the unique effects of non-uniform charge density on the particle's mobility in a viscoelastic fluid. To this end, it is demonstrated that the electrophoretic mobility (both translational and angular velocities) of the particle strongly depends on the non-uniformity of the surface charge, the relaxation and retardation times scales as well as the particle's curvature and may lead to breaking of fore–aft symmetry even for viscous flows in the low surface charge limits, which is in stark contrast to what is witnessed in Newtonian fluids. We successfully demonstrate that, under appropriate limiting conditions, previously reported results for Newtonian as well as viscoelastic fluids can be recovered as special cases of our theoretical framework.

The rest of the paper is organized as follows. In § 2, we provide a basic outline of the problem statement, the fundamental governing equations and the essential force and torque balance around the particle. In § 3, the modified Smoluchowski slip velocity is computed in the thin EDL limit for effectively weak viscoelastic flows. A more complete theoretical foundation, which includes the effects of rotation at higher orders of expansion, is presented in the supplementary material available at https://doi.org/10.1017/jfm.2021.643. In § 4, a model example of a non-uniformly charged particle is considered and its electrophoretic mobility is obtained based on the developed theory. Further, subsequent comparisons with reported results are carried out, to benchmark the theoretical framework as well as to bring out some of the novel insights specific to our study. In § 5, we explore how viscoelasticity changes the particle's angular velocity at a given instant, wherein the particle carries a non-axisymmetric surface charge. In § 6, we shed light on a potential experimental set-up, which may be used to validate some of the key predictions of our analysis. Finally, in § 7, concluding remarks are presented.

## 2. The physical paradigm and the governing equations

### 2.1. Description of the system

The prototypical system, as shown schematically in figure 1, consists of a spherical particle of radius $a$, carrying an arbitrary non-uniform surface charge of density $\sigma '(\theta ,\varphi )$, suspended in a viscoelastic medium of viscosity $\eta$, permittivity $\varepsilon$, relaxation time (Bird *et al.* Reference Bird, Armstrong and Hassager1987) $\lambda '_1$ and a retardation time $\lambda '_2$. The fluid also contains dissolved electrolytes, which dissociate into ions and form an EDL around the particle surface (Ohshima Reference Ohshima2006; Poddar *et al.* Reference Poddar, Maity, Bandopadhyay and Chakraborty2016), as shown in the schematic. The electrolyte concentration away from the particle is $c'_0$. The fluid is assumed to obey the Oldroyd-B constitutive relation. A uniform electric field (i.e. uniform only far away from the particle) of magnitude $E_0$ is applied to actuate motion. Without loss of generality, we may choose the direction of the applied field to be the $z$-axis. The variables around the particle are expressed using a spherical polar coordinate ($r',\theta ,\varphi$), with the origin fixed at the particle's centre, translating but not rotating with it.

As a result of the electrical forces acting on the particle, it will move with a velocity $\mathcal {U}'{\hat {\boldsymbol {e}}}_u$, where ${\hat {\boldsymbol {e}}}_u$ is the unit vector pointing towards the direction of the particle's motion. In general, ${\hat {\boldsymbol {e}}}_u$ is not constrained to be oriented along the $z$-axis (Anderson Reference Anderson1985). In addition, because of non-uniform surface charge, the particle may also undergo rotational motion (Anderson Reference Anderson1985; Yoon Reference Yoon1991; Solomentsev & Anderson Reference Solomentsev and Anderson1994) with angular velocity $\boldsymbol {\varOmega }'$, much like other anisotropic particles.

Noting that analysis of electrophoretic motion is quite a formidable problem (Saville Reference Saville1977), several physically realistic assumptions may be made (Ohshima Reference Ohshima1996, Reference Ohshima2015; Schnitzer *et al.* Reference Schnitzer, Zeyde, Yavneh and Yariv2013), which can simplify the underlying governing equations significantly, improving analytical tractability in the process. In the same spirit, we also apply a few simplifying assumptions, keeping in mind that our main goal here is to assess the role of viscoelasticity in the presence of non-uniform surface charge.

First, we shall assume the surface charge density to be weak for analytical tractability without sacrificing the essential physics of interest –we quantify the specific regime qualifying this ‘weak’ surface charge limit later. Second, we assume the characteristic EDL thickness ($\lambda _D$) to be much smaller than the particle's radius ($a$), which is a reasonably valid assumption for most practical considerations, since the EDL thickness typically lies in the tune of ${\sim }1 - 100\ \textrm {nm}$ (Ajdari Reference Ajdari1995; Bandopadhyay & Chakraborty Reference Bandopadhyay and Chakraborty2012*a*). Third, we assume that the ionic Péclet number (denoted by $Pe = u_c a/D'$; $u_c$ being the characteristic velocity, as defined later, $D'$ being the ionic diffusivity), as well as the dimensionless applied external field (as defined later), to be of the order of unity or less. This is in line with the consideration that for most particles, $Pe\sim O(1)$ (Saville Reference Saville1977). Fourth, the dimensionless relaxation time of the fluid medium as expressed in terms of the nominal Deborah number ($De$, see § 2.2 for definition) remains $O(1)$ or less. Despite invoking this consideration, later in § 3.4 we establish that, because of the weak surface charge limit (see the first assumption), the dimensionless velocity scale itself is less than $O(1)$, and hence the effective Deborah number (denoted as $De_{eff}$ later) is rendered much smaller than unity. As a consequence, the overall flow is only weakly viscoelastic, i.e. linear contributions dominate the leading-order flow stresses. Hence, the mathematical analysis carried out herein is very similar to an ordered-fluid expansion around the Newtonian limit (Chan & Leal Reference Chan and Leal1979). However, the Oldroyd-B constitutive relation is used because of its larger range of applicability as compared with an ordered fluid, which remains valid only for low strain rates (Bird *et al.* Reference Bird, Armstrong and Hassager1987). Fifth, we shall disregard the presence of a depletion layer (Pranay, Henríquez-Rivera & Graham Reference Pranay, Henríquez-Rivera and Graham2012; Mukherjee *et al.* Reference Mukherjee, Das, Dhar, Chakraborty and DasGupta2017) and assume the entire medium to exhibit viscoelastic characteristics. This requires that the characteristic EDL thickness be much larger than the polymer characteristic dimensions so that the continuum approximation remains valid inside the EDL. Towards assessing the validity of this assumption, one may refer to the radius of gyration ($R_g$) of the polymer as its representative length scale, which is a measure of the length scale of its random coil shapes (Israelachvili Reference Israelachvili2011). To compare physically realistic values of $R_g$ with those of $\lambda _D$ (the characteristic EDL thickness), we consider the example of poly-ethylene glycol in water (Cruje & Chithrani Reference Cruje and Chithrani2014). With a monomer length of ${\approx }0.35\ \textrm {nm}$ for $n = 100$ segments, one arrives at $R_g \approx 1.43\ \textrm {nm}$; for $n = 1000$ segments, $R_g \approx 4.5\ \textrm {nm}$. Therefore, for $\lambda _D \sim O(10\ \text {nm})$, the EDL thickness is an order of magnitude larger than the radius of gyration, which should safely allow us to assume a continuum distribution of polymers inside the Debye layer, as has been considered in a number of previous studies (Afonso, Alves & Pinho Reference Afonso, Alves and Pinho2009; Li & Koch Reference Li and Koch2020). Further, in practice, the polymer molecule includes only a tiny fraction of the random coil structure, which further justifies this continuum approximation. Finally, for a rotating particle (i.e. $\boldsymbol {\varOmega } \neq 0$), its surface charge distribution may be time variant. Although in § 5 we shall consider rotating particles, we will neglect such a dynamically evolving surface charge for the theoretical derivations. Nevertheless, an outline of time variation of the surface charge due to rotation has been provided in § S1 of the supplementary material accompanying this manuscript.

### 2.2. The characteristic scales

Electrophoretic motion entails several important characteristic scales for a number of different variables, dictating the flow dynamics. For any variable $\xi '$, we denote it's characteristic scale by $\xi _c$. As such, the following characteristic scales are chosen: the characteristic potential: $\psi _c = kT/e$ (where *k* is the Boltzmann constant, *T* is the absolute temperature and *e* the protonic charge) – this is the thermal potential; the characteristic length: $r_c = a$; the characteristic velocity: (Saville Reference Saville1977), $u_c = {\varepsilon k^{2} T^{2}}/{e^{2}\eta a}$; the characteristic concentration: $c_c = c'_0$. Based on these choices, the characteristic stress reads: $\tau _c = \eta u_c/a$, whereas the characteristic surface charge is $\sigma _c$. Finally, the characteristic relaxation time may be equated to $\lambda _0 =\max (\lambda '_1,\lambda '_2)$.

A number of important non-dimensional numbers emerge from the above characteristic scales, which strongly influence the flow around the particle. These include: (i) the characteristic EDL thickness, $\lambda _D = \sqrt {{\varepsilon k T}/{2c'_0 e^{2}}}$, with $\delta = \lambda _D/a$ – as the non-dimensional EDL thickness; (ii) the characteristic surface potential of the particle (Ajdari Reference Ajdari1995), $\zeta _c = \sigma _c \lambda _D/\varepsilon$, with $\bar {\zeta }_0 = e\zeta _c/kT$ as the characteristic potential relative to the thermal potential; (iii) the relative strength of the external field may be expressed through (Saville Reference Saville1977) $\beta = e E_0 a/kT$; (iv) the ionic Péclet number may be defined as (Saville Reference Saville1977) $Pe = u_c a/D$; (v) the nominal Deborah number, which indicates the extent of departure from Newtonian behaviour, may be defined as $De = u_c \lambda _0/a$. While various alternative approaches to quantify the extent of viscoelasticity have been reported (Li & Koch Reference Li and Koch2020); also see § 4.3 for further discussion, the advantage of defining $De$ using $u_c$ mentioned as above lies in the fact that various important asymptotic limits (such as weak surface charge, weak field limit, thin EDL limit etc.) may be independently explored without necessitating alteration to the magnitude of the Deborah number under consideration. In other words, the effect of electrokinetics may be elegantly isolated without requiring us to alter the quantitative representation of the fluid constitution.

Various assumptions stated in § 2.1 are now quantified in terms of the relevant non-dimensional numbers. ‘Weak surface charge’ indicates, $\bar {\zeta }_0 \ll 1$; ‘thin EDL’ implicates $\delta \ll 1$; the strengths of advection and applied field are constrained by: $\beta \sim O(1)$ and $Pe\sim O(1)$. It needs to be clarified here that, despite $De \sim O(1)$, the limit of weak surface charge ($\bar {\zeta }_0\ll 1$) implies that the dimensionless electrokinetic velocity would scale as $O(\bar {\zeta }_0)$. Therefore, the confluence of electromechanics and hydrodynamics would lead to an effective Deborah number given by $De_{eff}\sim O(\bar {\zeta }_0 De) \ll 1$ – see §§ 3.4 and 3.5. Accordingly, the effect of viscoelasticity remains subdominant in the theoretical derivations of the flow field under weak surface charge limits, despite the nominal $De$ being ${\sim }O(1)$. As a consequence, it is not necessary to further assume low polymer concentration (Li & Koch Reference Li and Koch2020), which would require $\lambda '_1/\lambda '_2-1\ll 1$ – see § 3.5 for further discussion of this.

### 2.3. The governing equations

The transport processes are governed by the Poisson–Nernst–Planck–Cauchy momentum equations (Saville Reference Saville1977; Ghosh *et al.* Reference Ghosh, Chaudhury and Chakraborty2016), along with the continuity equation for conservation of mass and the Oldroyd-B constitutive equations. Since these equations are well established, we directly start with their dimensionless forms, wherein the non-dimensional version of the any variable $\xi '$ is expressed as $\xi = \xi '/\xi _c$. The Nernst–Planck equations are reformulated in terms of the net charge density and concentration, defined as (Schnitzer & Yariv Reference Schnitzer and Yariv2014) $c = c_{+}+c_{-}$ and $\rho = c_{+}-c_{-}$ respectively, where $c_{+(-)}$ is the concentration of the positive (negative) ions. In view of the characteristic scales outlined in § 2.2, the governing equations are written as (Saville Reference Saville1977; Bird *et al.* Reference Bird, Armstrong and Hassager1987)

*a*)\begin{gather} Pe({\boldsymbol{v}}\boldsymbol{\cdot}\boldsymbol{\nabla} c) = \nabla^{2} c+\boldsymbol{\nabla} \boldsymbol{\cdot}\left\{\rho\boldsymbol{\nabla}\psi\right\}, \end{gather}

*b*)\begin{gather}Pe({\boldsymbol{v}}\boldsymbol{\cdot}\boldsymbol{\nabla} \rho) = \nabla^{2} \rho+\boldsymbol{\nabla} \boldsymbol{\cdot}\left\{c\boldsymbol{\nabla}\psi\right\} , \end{gather}

*c*)\begin{gather}\delta^{2}\nabla^{2}\psi ={-}\tfrac{1}{2}\rho, \end{gather}

*d*)\begin{gather}-\boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tau} +\nabla^{2}\psi\boldsymbol{\nabla} \psi = 0 \quad \text{and} \quad \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}} = 0, \end{gather}

*e*)\begin{gather}\boldsymbol{\tau} + De\lambda_1 \overset{\nabla}{\boldsymbol{\tau}} = 2{\boldsymbol{D}} + 2\lambda_2 De \overset{\nabla}{{\boldsymbol{D}}} . \end{gather}

In the above equations, $\boldsymbol {\tau }$ is the stress field, $p$ is the pressure, $\psi$ is the electrical potential, ${\boldsymbol {D}}$ is the rate of strain tensor (${=}\frac {1}{2}[\boldsymbol {\nabla }{\boldsymbol {v}} + (\boldsymbol {\nabla }{\boldsymbol {v}})^\textrm {T}]$) and ${\boldsymbol {v}}$ is the velocity field. Further, $\overset {\nabla }{\boldsymbol {\tau }}$ and $\overset {\nabla }{{\boldsymbol {D}}}$ indicate the first convected derivatives of the stress and strain rate. For any second-rank tensor ${\boldsymbol {A}}$, it's convected derivative $\overset {\nabla }{{\boldsymbol {A}}}$ is expressed as (Bird *et al.* Reference Bird, Armstrong and Hassager1987) $\overset {\nabla }{{\boldsymbol {A}}} = {\textrm {D}{\boldsymbol {A}}}/{\textrm {D}t}-(\boldsymbol {\nabla }{\boldsymbol {v}})^\textrm {T}\boldsymbol {\cdot }{\boldsymbol {A}}-{\boldsymbol {A}}\boldsymbol {\cdot }\boldsymbol {\nabla }{\boldsymbol {v}}$. The above equations are subject to the following boundary conditions:

*a*)\begin{gather} {\hat{\boldsymbol{e}}}_r\boldsymbol{\cdot}\left[Pe{\boldsymbol{v}} c - \boldsymbol{\nabla} c -\rho\boldsymbol{\nabla}\psi\right]_{r=1} = {\hat{\boldsymbol{e}}}_r\boldsymbol{\cdot}\left[Pe{\boldsymbol{v}} \rho - \boldsymbol{\nabla} \rho -c\boldsymbol{\nabla}\psi\right]_{r=1} = 0 \end{gather}

*b*)\begin{gather}{\hat{\boldsymbol{e}}}_r\boldsymbol{\cdot}[\boldsymbol{\nabla}\psi]_{r=1} ={-}\frac{\zeta(\theta,\varphi)}{\delta};\quad \zeta = \frac{\sigma'\lambda_D e}{\varepsilon kT}\quad \text{and}\quad [\boldsymbol{\nabla}\psi]_{r\rightarrow\infty} ={-}\beta{\hat{\boldsymbol{e}}}_z \end{gather}

*c*)\begin{gather}{\boldsymbol{v}}(r=1,\theta,\varphi) = \boldsymbol{\varOmega}\times {\hat{\boldsymbol{e}}}_r\quad \text{and}\quad {\boldsymbol{v}}(r\rightarrow\infty,\theta,\varphi) ={-}\mathcal{U}{\hat{\boldsymbol{e}}}_u \end{gather}

*d*)\begin{gather}\text{at}\ r\rightarrow \infty, \ c = 2 \quad \text{and}\quad \rho = 0. \end{gather}

Note that, in (2.2*c*), the boundary condition for velocity at the particle surface accounts for its rotation. In (2.2*b*), $\zeta (\theta ,\varphi )$ is analogous to the surface potential and $\zeta \sim O(\bar {\zeta }_0)$ is it's characteristic magnitude. For convenience, we define, $\zeta (\theta ,\varphi ) = \bar {\zeta }_0\bar {\zeta }(\theta ,\varphi )$, such that $\bar {\zeta } \sim O(1)$. For electrophoretic motion, $\mathcal {U}{\hat {\boldsymbol {e}}}_u$ and $\boldsymbol {\varOmega }$, i.e. the migration velocity of the particle and it's angular velocity are the primary unknowns. To compute them, we simply need to balance the force and the moment around the origin, acting on the particle at steady state. The resulting balances are as follows (Anderson Reference Anderson1985; Ohshima Reference Ohshima2006; Goswami *et al.* Reference Goswami, Dhar, Ghosh and Chakraborty2017):

*a*)\begin{gather} -\frac{\bar{\zeta}_0}{\delta}\int_{S_p} \bar{\zeta}\boldsymbol{\nabla}\psi \,\textrm{d}S + \int_{S_p} \boldsymbol{\tau}\boldsymbol{\cdot}{\hat{\boldsymbol{e}}}_r \,\textrm{d}S = 0, \end{gather}

*b*)\begin{gather}-\frac{\bar{\zeta}_0}{\delta}\int_{S_p} \bar{\zeta}({\hat{\boldsymbol{e}}}_r\times\boldsymbol{\nabla}\psi) \,\textrm{d}S + \int_{S_p} {\hat{\boldsymbol{e}}}_r\times\left(\boldsymbol{\tau}\boldsymbol{\cdot}{\hat{\boldsymbol{e}}}_r\right) \textrm{d}S = 0 . \end{gather}

It is important to note here that the above ‘force-free’ conditions only remain valid in an unbounded medium, as previously noted by Yariv (Reference Yariv2006). Presence of a wall sufficiently close to the particle may change the force balance on the same.

The first terms in the above equations represent the force and the moment, respectively, due to the electric field. The second terms are essentially contributed by the stresses in the fluid. The integration is carried out over the particle surface $S_p$. In the limit of a thin EDL, i.e. for $\delta \rightarrow 0$, the EDL may also be included within $S_p$ so that the integration is carried out on a surface just outside the EDL (as $\delta \rightarrow 0$, the EDL effectively has zero thickness). Since the net charge in the EDL and on the particle surface are opposite and equal, the net charge on an imaginary surface just outside of the EDL would be zero. Therefore, the first terms in both (2.3*a*) and (2.3*b*) vanish (Ye *et al.* Reference Ye, Sinton, Erickson and Li2002; Chen & Keh Reference Chen and Keh2014) and, as a result, the net force and moment caused by the stresses in the fluid on the imaginary sphere lying just outside the EDL (i.e. on the object $\textrm {particle}+\textrm {the}$ EDL) becomes zero (Ye *et al.* Reference Ye, Sinton, Erickson and Li2002; Chen & Keh Reference Chen and Keh2014). Equations (2.3*a*) and (2.3*b*) may then be re-written as

*a*,

*b*)\begin{equation} \int_{S_p} \boldsymbol{\tau}\boldsymbol{\cdot}{\hat{\boldsymbol{e}}}_r \,\textrm{d}S = 0\quad \text{and}\quad \int_{S_p} {\hat{\boldsymbol{e}}}_r\times (\boldsymbol{\tau}\boldsymbol{\cdot}{\hat{\boldsymbol{e}}}_r) \,\textrm{d}S = 0. \end{equation}

Notice that the integration is to be performed over the surface $S_p$, which has a radius $r \sim 1+\delta$ and includes the EDL as well.

Finally, we note that the potential may be conveniently split into two parts as follows (Bahga, Vinogradova & Bazant Reference Bahga, Vinogradova and Bazant2010; Ghosh, Mandal & Chakraborty Reference Ghosh, Mandal and Chakraborty2017): $\psi = \phi + \phi _{ext}$, where $\phi _{ext}$ is the contribution from only the externally imposed electric field and $\phi$ is the contribution of the particle's surface charge to the total electrostatic potential. It may then be deduced that $\phi _{ext}$ satisfies the equation (Ghosh *et al.* Reference Ghosh, Mandal and Chakraborty2017) $\nabla ^{2}\phi _{ext} = 0$, subject to ${\hat {\boldsymbol {e}}}_r\boldsymbol {\cdot }[\boldsymbol {\nabla }\phi _{ext}]_{r=1} = 0$ and $[\boldsymbol {\nabla }\phi _{ext}]_{r\rightarrow \infty } = -\beta {\hat {\boldsymbol {e}}}_z$. Hence, $\phi _{ext}$ has the solution

where, $\mu = \cos \theta$ and $P_{n}(x)$ is the Legendre polynomial of the first kind of order $n$. As a consequence, the component $\phi$ would satisfy the following equation (derived from (2.1*c*)):

subject to the following conditions:

*a*,

*b*)\begin{equation} {\hat{\boldsymbol{e}}}_r\boldsymbol{\cdot}[\boldsymbol{\nabla}\phi]_{r=1} ={-}\frac{\zeta(\theta,\varphi)}{\delta}\quad \text{and}\quad \phi(r\rightarrow\infty,\theta,\varphi) = 0 .\end{equation}

The Nernst–Planck equations (2.1*a*) and (2.1*b*) and the corresponding boundary conditions (2.2*a*) may also be similarly modified using the above split, which is not explicitly carried out here for the sake of brevity.

When particle rotation is present, it will cause the surface potential ($\zeta$) to change dynamically. The relation between the particle's angular velocity and evolution of $\zeta$ has been included in § S1.1 of the supplementary material.

## 3. The modified Smoluchowski slip in the thin EDL limit

In this section, the electro-hydrodynamics around the particle is analysed asymptotically. When the EDL is thin, the fluid only experiences an electrical force in a very small region near the particle surface, wherein all the diffuse charges are located. At the same time, the effect of the particle's surface charge and motion in the EDL is transmitted into the bulk through the Smoluchowski slip velocity (Ajdari Reference Ajdari1995; Squires & Bazant Reference Squires and Bazant2004; Ghosh *et al.* Reference Ghosh, Chaudhury and Chakraborty2016). Note that this Smoluchowski slip velocity is nothing but the velocity tangent to the solid surface at the outer edge of the EDL. This slip velocity is what drives the motion in the bulk and therefore dictates the viscous resistance exerted by the bulk at the edge of the EDL. As a result, it becomes necessary to first analyse the transport within the EDL.

From (2.1*c*), we note that the thin EDL limit ($\delta \ll 1$) is a singular problem, as also pointed out by a number of earlier studies (Yariv Reference Yariv2009; Schnitzer & Yariv Reference Schnitzer and Yariv2012; Ghosh *et al.* Reference Ghosh, Chaudhury and Chakraborty2016). This calls for the application of a matched asymptotic expansion (Leal Reference Leal2007; Bender & Orszag Reference Bender and Orszag2013), wherein the fluid domain is split into two regions with two distinct length scales, (*a*) the EDL, or ‘inner layer’, which has a characteristic length scale $\delta$ and lies next to the particle surface and (*b*) the bulk, or ‘outer layer’, with characteristic length scale $r\sim O(1)$. In both the layers, all variables may be expanded in an asymptotic series of $\delta$ as follows (Leal Reference Leal2007; Yariv Reference Yariv2009):

Note that this expansion also applies to $\mathcal {U}$ and $\boldsymbol {\varOmega }$. The flow variables have to be matched asymptotically at the edge of the EDL, where the two regions meet. Keeping in mind that our aim is to first deduce the modified Smoluchowski slip, it would suffice to only consider the leading term in the expansion (3.1).

To keep our analysis structured and generic, we shall first outline the governing equations in the outer region in § 3.1 and the rescaled equations in the inner region in § 3.2, followed by appropriate matching conditions in § 3.3. This analysis is presented without any restriction on the surface charge. However, we consider $\delta \ll 1$ and $De\sim O(1)$ for writing the governing equations in the two layers. As we show later, the inner layer equations thus derived are distinct as compared with Newtonian fluids and they vividly bear the consequences of viscoelasticity. Subsequently in § 3.4, we shall invoke the assumption of weak surface charge ($\bar {\zeta }_0\ll 1$), which will enable us to pin down regular asymptotic solutions for the flow field inside the EDL and will lead to closed-form expressions for the modified Smoluchowski slip for an arbitrary distribution of surface charge. Solutions to the flow field in the outer region for specific instances of weak but non-uniform surface charge are reported in §§ 4 and 5, respectively. Because we will only consider the leading-order terms in the expansion (3.1) for all variables, we will drop the ‘0’ subscript from here onwards.

### 3.1. Leading-order equations in the outer layer

Leading-order outer layer equations may be derived by inserting the expansion (3.1) into the set of equations as outlined in (2.1*a*)–(2.1*e*), which take the following form:

*a*)\begin{gather} \rho = 0;\quad \boldsymbol{\nabla} \boldsymbol{\cdot}\left(c\boldsymbol{\nabla}(\phi+\phi_{ext})\right)= 0 \end{gather}

*b*)\begin{gather}{\boldsymbol{v}}\boldsymbol{\cdot}\boldsymbol{\nabla} c = {\nabla}^{2} c \end{gather}

*c*)\begin{gather}-\boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tau} +\nabla^{2}\phi\boldsymbol{\nabla} (\phi + \phi_{ext}) = 0\quad \text{and}\quad \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}} = 0 \end{gather}

*d*)\begin{gather}\boldsymbol{\tau} + De\lambda_1 \boldsymbol{\mathcal{T}} = 2{\boldsymbol{D}} + 2\lambda_2 De \boldsymbol{{\mathcal{S}}}. \end{gather}

In (3.2*d*), $\boldsymbol {\mathcal {T}} = \overset {\nabla }{\boldsymbol {\tau }}$ and $\boldsymbol {{\mathcal {S}}} = \overset {\nabla }{{\boldsymbol {D}}}$. Detailed expressions for $\boldsymbol {\mathcal {T}}$ and $\boldsymbol {{\mathcal {S}}}$ in spherical coordinates may be found in Bird *et al.* (Reference Bird, Armstrong and Hassager1987). Note that, $\phi _{ext}$ has already been determined and thus there is no need to expand it in $\delta$. The above equations are subject to the following far field conditions:

The boundary conditions at the edge of the EDL (where two layers meet) will be given in terms of the matching conditions later. Note that the outer layer is electroneutral at the leading order of $\delta$, and hence it would not experience any force due to the externally applied electric field (as shown later).

### 3.2. Leading-order equations in the EDL

A thorough rescaling of almost all quantities is required in the inner layer, so that the rescaled variables are $O(1)$ in the EDL. To this end, we shall first introduce the rescaled radial coordinate (Schnitzer & Yariv Reference Schnitzer and Yariv2012), $R = (r-1)/\delta$, which is $O(1)$ inside the inner layer, since $r - 1 \sim \delta$. Accordingly, the velocity components would rescale as follows: $u_{\theta } \rightarrow U \sim O(1)$, $u_{\varphi } \rightarrow W \sim O(1)$ and from the continuity equation, $u_{r}\rightarrow \delta V \sim O(\delta )$, with $V\sim O(1)$. In the above, $u_{r}$, $u_{\theta }$ and $u_{\varphi }$ are the radial, polar and azimuthal components of velocity, respectively. There are no changes in the scaling of the potential, charge and salt concentration: $\varphi \rightarrow \varPhi \sim O(1)$, $c\rightarrow C\sim O(1)$ and $\rho \rightarrow \varPi \sim O(1)$. The rescaling of the stresses and the strain rates are slightly more involved (Saprykin, Koopmans & Kalliadasis Reference Saprykin, Koopmans and Kalliadasis2007; Ghosh & Chakraborty Reference Ghosh and Chakraborty2015; Ghosh *et al.* Reference Ghosh, Chaudhury and Chakraborty2016) and need to be worked out based on the constitutive model (2.1*e*). The order of magnitude of the variables and their rescaled versions in the inner layer (the EDL) are summarized in table 1. Here, the rescaled stress, strain and their convected derivatives are denoted by a ‘tilde’ overhead in the inner layer, whereas the rescaled versions of the primitive variables (velocity, potential, concentration etc.) are denoted by an uppercase symbol.

There is significant difference between the rescaling mentioned in table 1 and that of Newtonian fluids. In a Newtonian medium, the shear stress components (such as $\tau _{r\theta }$) would scale as $O(\delta ^{-1})$. This remains unchanged for viscoelastic fluids. However, the normal stress components always scale as $O(1)$ inside the EDL for Newtonian fluids. In stark contrast, for viscoelastic fluids, the normal stresses (such as $\tau _{\theta \theta }$) scale as $O(\delta ^{-2})$ in the inner layer. In addition, similar scaling is also observed for the component $\tau _{\theta \varphi }$. These unusual scalings would have strong implications for the flow dynamics inside the EDL, as is discussed later in more detail. The scaling mentioned above is also supported by a few of the previous studies (Saprykin *et al.* Reference Saprykin, Koopmans and Kalliadasis2007; Ghosh & Chakraborty Reference Ghosh and Chakraborty2015; Ghosh *et al.* Reference Ghosh, Chaudhury and Chakraborty2016), where asymptotically large normal stresses as compared with the shear stresses were reported, albeit for significantly more restrictive flat geometries. Furthermore, note that the components of $\boldsymbol {{\mathcal {S}}}$ and ${\boldsymbol {D}}$ do not exhibit the same scaling, contrary to the stress ($\boldsymbol {\tau }$) and its convected derivative ($\boldsymbol {\mathcal {T}}$). We insert the rescaled variables as described in table 1 in the governing equations, from which the leading-order inner layer equations may be obtained as follows:

*a*)\begin{gather} \frac{\partial^{2}\varPhi}{\partial R^{2}} ={-}\frac{1}{2}\varPi \end{gather}

*b*)\begin{gather} \frac{\partial^{2} C}{\partial R^{2}} + \frac{\partial}{\partial R}\left(\varPi\frac{\partial\varPhi}{\partial R}\right) = \frac{\partial^{2} \varPi}{\partial R^{2}} + \frac{\partial}{\partial R}\left(C\frac{\partial\varPhi}{\partial R}\right)=0\end{gather}

*c*)\begin{gather} \frac{\partial V}{\partial R} - \frac{\partial}{\partial \mu}\left(U\sqrt{1-\mu^{2}}\right) +\frac{1}{\sqrt{1-\mu^{2}}}\frac{\partial W}{\partial\varphi}= 0\\ -\frac{\partial P}{\partial R} + \frac{\partial^{2}\varPhi}{\partial R^{2}}\frac{\partial\varPhi}{\partial R} = 0 \end{gather}

*d*)\begin{gather} \sqrt{1-\mu^{2}}\frac{\partial P}{\partial\mu}+\frac{\partial\tilde{\tau}_{r\theta}}{\partial R}-\frac{\partial}{\partial \mu}\left(\tilde{\tau}_{\theta\theta}\sqrt{1-\mu^{2}}\right) +\frac{1}{\sqrt{1-\mu^{2}}}\frac{\partial\tilde{\tau}_{\theta\varphi}}{\partial \varphi} - \frac{\mu\tilde{\tau}_{\varphi\varphi}}{\sqrt{1-\mu^{2}}}\nonumber\\ +\frac{3}{2}\beta\sqrt{1-\mu^{2}}\frac{\partial^{2}\varPhi}{\partial R^{2}} - \sqrt{1-\mu^{2}}\frac{\partial^{2}\varPhi}{\partial R^{2}}\frac{\partial\varPhi}{\partial \mu} = 0 \end{gather}

*f*)\begin{gather} -\frac{1}{\sqrt{1-\mu^{2}}}\frac{\partial P}{\partial\varphi}+\frac{\partial\tilde{\tau}_{r\varphi}}{\partial R} - \frac{\partial}{\partial\mu}\left(\tilde{\tau}_{\theta\varphi}\sqrt{1-\mu^{2}}\right) + \frac{1}{\sqrt{1-\mu^{2}}}\frac{\partial\tilde{\tau}_{\varphi\varphi}}{\partial\varphi}\nonumber\\ + \frac{\mu\tilde{\tau}_{\theta\varphi}}{\sqrt{1-\mu^{2}}} + \frac{\partial^{2}\varPhi}{\partial R^{2}} \frac{\partial\varPhi}{\partial \varphi} = 0. \end{gather}

Recall that $\mu = \cos \theta$. Also, from (2.5), in the inner layer, $E_{r}^{ext} = 0$ and $E_{\theta }^{ext} = -\frac {3}{2}\beta \sqrt {1-\mu ^{2}}$ to the leading order in $\delta$, where $\boldsymbol {E}^{ext} = -\boldsymbol {\nabla }\phi _{ext}$. This is why the external field does not influence the charge and salt transport at the leading order. The inner layer momentum equations (3.4*d*)–(3.4*f*) indicate the consequences of the scaling outlined earlier.

In order to better understand how the different stress components influence the flow within the EDL, we may try to appeal to some of the fundamental properties exhibited by viscoelastic fluids in simple flows. First note that the flow inside the EDL is locally unidirectional and hence is characterized by strong shear strain rates, because of its small thickness ($\delta \ll 1$). It is well documented (Bird *et al.* Reference Bird, Armstrong and Hassager1987) that, in perfectly unidirectional shear flows (such as Couette or Poiseuille flows), the shear stress ($\tau _{xy}$, $x$ being the direction of flow) varies as $\tau _{xy} = \eta _*(\dot {\gamma })\dot {\gamma }$, where $\dot {\gamma }$ is the rate of strain and $\eta _*(\dot {\gamma })$ is the appropriately non-dimensionalized shear dependent viscosity. On the other hand, the first normal stress difference varies as $\tau _{xx}-\tau _{yy} = \eta _{N,1}(\dot {\gamma })\dot {\gamma }^{2}$, where $\eta _{N,1}$ is the (unitless) first normal stress coefficient. Because we have considered the Oldroyd-B constitutive model, it follows that both $\eta _*$ and $\eta _{N,1}$ are constants (more discussion on this is provided in § 3.6). Since the shear rate in the EDL is $\dot {\gamma }\sim \delta ^{-1}$, it immediately follows that in the EDL on a perfectly plane surface $\tau _{xy} \sim \delta ^{-1}$ and $\tau _{xx} \sim \delta ^{-2}$, i.e. the stresses along the streamwise directions become very large, because the unidirectional flow with strong shear rate stretches the polymers along those directions. Now, in an EDL adhering to a spherical particle, locally, the streamwise directions are $\theta$ and $\varphi$, using spherical coordinates. Thus the stresses in the streamwise directions are $\tau _{\theta \theta }$, $\tau _{\varphi \varphi }$ and $\tau _{\theta \varphi }$, while stresses equivalent to $\tau _{xy}$ are $\tau _{r\theta }$ and $\tau _{r\varphi }$. As a result, we expect that $\tau _{r\theta },\tau _{r\varphi } \sim O(\delta ^{-1})$ and $\tau _{\theta \theta },\tau _{\varphi \varphi },\tau _{\theta \varphi }\sim O(\delta ^{-2})$, as is indeed verified from table 1. We would like to point out here that similar scalings of stress components in viscoelastic flows have been previously shown in earlier studies, albeit only for motion over flat surfaces (Saprykin *et al.* Reference Saprykin, Koopmans and Kalliadasis2007; Ghosh *et al.* Reference Ghosh, Chaudhury and Chakraborty2016).

On a perfectly flat surface with uniform surface charge, the streamwise stresses ($\tau _{xx}$ etc.) would be uniform and hence they would not affect the flow. However, on the surface of a spherical particle, the streamwise stresses (such as $\tau _{\theta \theta },\tau _{\varphi \varphi }$ etc.) would vary on a length scale of $O(1)$, provided that the surface charge on the particle also varies at a scale of $O(a)$ (dimensionless $O(1)$). Therefore, the gradients of the streamwise stresses would vary as $O(\delta ^{-2})$. At the same time, because the cross-stream gradients of $\tau _{r\theta }$ and $\tau _{r\varphi }$ govern the velocity field, these gradients also scale as $\delta ^{-2}$, despite $\tau _{r\theta }$ and $\tau _{r\varphi }$ scaling as $\delta ^{-1}$. As a result, the gradients of the extensional stresses and the shear stress both have the same order of magnitude inside the EDL and thus both together dictate the velocity field therein. The reasoning presented above physically justifies the equations (3.4) that govern the leading-order motion inside the EDL. This is distinct from Newtonian fluids, where the streamwise stresses do not contribute to the velocity field at the leading order of $\delta$. Further, note that, even when the surface charge is uniform, the streamwise gradients in the extensional stress components remain non-zero because of the particle's curvature. This argument indicates that a Newtonian fluid cannot ‘feel’ the curvature of the particle's surface at the leading order of $\delta$. The only way curvature affects the flow is through the $\theta$ component of the external electric field. In contrast, a viscoelastic fluid is able to ‘see’ the curvature of the particle surface even at the leading order of $\delta$, because of its asymptotically large normal stresses. One of the major consequences of such behaviour is that the particle's shape plays a crucial role in modifying the Smoluchowski slip velocity at the edge of the EDL, as shown later. The changes thus brought about in the slip velocity also alter the overall electrophoretic velocity of the particle, as discussed in detail in the next section.

The rescaled Oldroyd-B constitutive relation may be expressed in the inner layer as follows:

*a*)\begin{gather} \tilde{\tau}_{ij} + \lambda_1 De \tilde{\mathcal{T}}_{ij} = 2[\tilde{D}_{ij} + \lambda_2 De \tilde{{\mathcal{S}}}_{ij}],\quad \text{when}\ ij\equiv rr,\ r\theta\ \text{and}\ r\varphi, \end{gather}

*b*)\begin{gather}\tilde{\tau}_{ij} + \lambda_1 De \tilde{\mathcal{T}}_{ij} = 2\lambda_2 De \tilde{{\mathcal{S}}}_{ij},\quad \text{when}\ ij\equiv \theta\theta, \theta\varphi\ \text{and}\ \varphi\varphi . \end{gather}

The detailed expressions for the rescaled convected derivatives ($\tilde {\mathcal {T}}$ and $\tilde {{\mathcal {S}}}$) are given in Appendix A. Derivation of the above equations for a selected few stress components is outlined in § S1.3 in the supplementary material. The equations in the inner layer are subject to the following boundary conditions at the particle surface:

*a*)\begin{gather} \frac{\partial C}{\partial R}+\varPi\frac{\partial\varPhi}{\partial R} = \frac{\partial\varPi}{\partial R}+C\frac{\partial\varPhi}{\partial R} = 0 \end{gather}

*b*)\begin{gather}\frac{\partial\varPhi}{\partial R} ={-}\bar{\zeta}_0\bar{\zeta}(\theta,\phi) \end{gather}

*c*)\begin{gather}U = \varGamma(\theta,\varphi);\quad W = \chi(\theta,\varphi);\quad V = 0 \ \text{at},\quad R = 0. \end{gather}

In (3.6*c*), $\varGamma = (\boldsymbol {\varOmega }\times {\hat {\boldsymbol {e}}}_r)\boldsymbol {\cdot }{\hat {\boldsymbol {e}}}_{\theta } = \varOmega _y\cos (\varphi )-\varOmega _x\sin (\varphi )$ and $\chi = (\boldsymbol {\varOmega }\times {\hat {\boldsymbol {e}}}_r)\boldsymbol {\cdot }{\hat {\boldsymbol {e}}}_{\varphi } = -(\varOmega _x \cos (\varphi )+\varOmega _y\sin (\varphi ))\cos (\theta )+\varOmega _z\sin (\theta )$ are respectively the $\theta$ and $\varphi$ components of velocity at the particle surface due to it's rotation.

### 3.3. Asymptotic matching and the Smoluchowski slip

The matching conditions for the primitive variables (such as velocity, potential, charge density etc.) at the edge of the EDL (where the inner and outer regions overlap) are given by (Ghosh *et al.* Reference Ghosh, Chaudhury and Chakraborty2016, Reference Ghosh, Mandal and Chakraborty2017)

*a*)\begin{gather} \lim_{R\rightarrow\infty} [U,\delta V,W] = \lim_{r\rightarrow 1} [u_{\theta},u_{r},u_{\varphi}], \end{gather}

*b*)\begin{gather}\lim_{R\rightarrow\infty} [\delta^{{-}2}P,\varPhi,\varPi,C] = \lim_{r\rightarrow 1} [p,\varphi,\rho,c]. \end{gather}

In addition, the net salt and charge fluxes across the EDL–bulk interface also need to be matched (Yariv Reference Yariv2009; Ghosh *et al.* Reference Ghosh, Chaudhury and Chakraborty2016, Reference Ghosh, Mandal and Chakraborty2017) to ensure that the EDL does not lose or gain net charge or salt at steady state. It may be shown (see § S1.2 in the supplementary material for further details) that, at the leading order, the matching conditions mentioned above predict the following boundary conditions for bulk salt concentration and potential at the edge of the EDL (Yariv Reference Yariv2009; Ghosh *et al.* Reference Ghosh, Chaudhury and Chakraborty2016):

Finally, it is important to note that the quantities, $\lim _{R\rightarrow \infty } U$ and $\lim _{R\rightarrow \infty } W$ may be combined to write $\boldsymbol {v}_S = \lim _{R\rightarrow \infty }[U{\hat {\boldsymbol {e}}}_{\theta } + W{\hat {\boldsymbol {e}}}_{\phi }]-\boldsymbol {\varOmega }\times {\hat {\boldsymbol {e}}}_r$. We identify the quantity ${\boldsymbol {v}}_S$ as the ‘modified Smoluchowski slip’ at the edge of the EDL and it denotes the tangential slip velocity experienced by the outer layer fluid, owing to the presence of electrokinetic flow inside the EDL. Note that the slip velocity is defined relative to the particle surface. Once the slip velocity is known, the outer layer momentum and continuity equations, i.e. (3.2*c*) may be solved subject to (3.3) in the far field and ${\boldsymbol {v}} = \boldsymbol {\varOmega }\times {\hat {\boldsymbol {e}}}_r + {\boldsymbol {v}}_S$ and $u_r = 0$ at $r = 1$.

### 3.4. Analysis for weak surface charge

The analysis within the EDL would ultimately lead to the ‘modified Smoluchowski slip’, for which one needs to first solve the inner layer equations, subject to the relevant boundary and matching conditions. In order to derive the closed-form analytical solutions, it is necessary to assume the surface charges to be weak ($\bar {\zeta }_0 \ll 1$). Accordingly, all the variables (both in the inner and outer layers) may be further expanded (Ghosh *et al.* Reference Ghosh, Mandal and Chakraborty2017) in a regular asymptotic series of $\bar {\zeta }_0$ as

Here, $\xi$ may represent any variable such as $U,V,{\boldsymbol {v}},\phi ,\ldots$ and so on. Recall that the above expansion is applied to the variables which already denote leading-order terms in $\delta$. We re-emphasize that, in the regular expansion, $\bar {\zeta }_0$ is defined in terms of the characteristic surface charge (see § 2.2), which is assumed to be small here. Using (3.9), in the subsequent subsections we shall determine the modified Smoluchowski slip for arbitrary distribution of surface charge. In the next section, the slip velocity thus derived will be used in a representative example of a non-uniformly charged particle to determine its electrophoretic mobility, by solving the outer layer equations.

#### 3.4.1. Simplified outer layer equations

From the conditions in (3.8), it is easy to deduce using (3.2*a*) and (3.2*b*) that the solutions for the potential and concentration in the outer layer are (at the leading order of $\delta$): $c = 2$ and $\phi = 0$. The equations governing the fluid motion in the outer layer then take the following form:

*a*)\begin{gather} -\boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tau} = 0\quad \text{and}\quad \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}} = 0, \end{gather}

*b*)\begin{gather}\boldsymbol{\tau} + De\lambda_1 \boldsymbol{\mathcal{T}} = 2{\boldsymbol{D}} + 2\lambda_2 De \boldsymbol{{\mathcal{S}}}. \end{gather}

These are subject to the following boundary conditions:

*a*)\begin{gather} \text{at}\ r\rightarrow \infty,\quad {\boldsymbol{v}} ={-}\mathcal{U}{\hat{\boldsymbol{e}}}_u, \end{gather}

*b*)\begin{gather}\text{at}\ r = 1,\quad u_r = 0;\quad u_{\theta} = \lim_{R\rightarrow\infty} U;\quad\text{and}\quad u_{\varphi} = \lim_{R\rightarrow\infty} W . \end{gather}

#### 3.4.2. Solutions to the inner layer equations

The leading order (in $\delta$) inner layer equations may also be solved using the asymptotic expansion mentioned in (3.9), subject to conditions (3.6*a*)–(3.6*c*). To this end, we note that the inner layer variables exhibit the following expansions in $\bar {\zeta }_0$:

*a*)\begin{gather} C = 2 + \bar{\zeta}_0 C_1 + \bar{\zeta}_0^{2} C_2 + \ldots\,;\quad (\varPi,\varPhi) = \bar{\zeta}_0 (\varPi_1,\varPhi_1) + \bar{\zeta}_0^{2}(\varPi_2,\varPhi_2) + \ldots \end{gather}

*b*)\begin{gather}(U,V,W) = \bar{\zeta}_0 (U_1,V_1,W_1) + \bar{\zeta}_0^{2} (U_2,V_2,W_2) + \bar{\zeta}_0^{3} (U_3,V_3,W_3) + \ldots \end{gather}

*c*)\begin{gather}(P,\boldsymbol{\tilde{\mathcal{T}}},\boldsymbol{\tilde{{\mathcal{S}}}}) = \bar{\zeta}_0^{2}(P_2, \boldsymbol{\tilde{\mathcal{T}}}_2,\boldsymbol{\tilde{{\mathcal{S}}}}_2) + \bar{\zeta}_0^{3}(P_3, \boldsymbol{\tilde{\mathcal{T}}}_3,\boldsymbol{\tilde{{\mathcal{S}}}}_3) + \ldots. \end{gather}

Implications and the physical basis of the asymptotic expansions mentioned above deem further elaboration. Notice that, at the leading order of $\bar {\zeta }_0$, the nonlinear polymeric contributions to the stresses are absent, which indicates that at $O(\bar {\zeta }_0)$, the flow inside the EDL remains asymptotically Newtonian, despite the viscoelastic stresses appearing in (3.4) and (3.5). This apparent contradiction may be resolved by noting that in the weak surface charge limits, the dimensionless velocity ($U$ and $W$) inside the EDL scales as $O(\bar {\zeta }_0)$. As a result, the rescaled Newtonian part of the stresses (also refer to table 1) would scale as $O(\bar {\zeta }_0)$ and the convected derivatives ($\tilde {\mathcal {T}}$ and $\tilde {{\mathcal {S}}}$), which arise from the polymeric contributions to the stresses, would scale as $O(\bar {\zeta }_0^{2})$, as evident from their expressions in Appendix A. Here we would like to clarify that the ‘Newtonian part’ of the stresses as mentioned above refers to the linear part of the constitutive relation given in (2.1*e*), obtained by substituting $\lambda _1 = \lambda _2 = 0$. Therefore, at the leading order of $\bar {\zeta }_0$, i.e. at $O(\bar {\zeta }_0)$, the flow is asymptotically Newtonian and the nonlinearities arising from the polymeric stresses only contribute from $O(\bar {\zeta }_0^{2})$ onwards, stemming from the corresponding convected derivatives.

Following the discussion after (3.4), we note that inside the EDL, the linear (or Newtonian) parts of the shear stress gradients (such as ${\partial \tau _{r\theta }}/{\partial r}$) scale as $O(\bar {\zeta }_0\delta ^{-2})$ for weak surface charge, since inside the EDL, $r\sim 1+\delta$ and $\tau _{r\theta },$ etc. $\sim O(\delta ^{-1})$. On the other hand, the normal stress derivatives (terms like ${\partial }/{\partial \mu }(\tau _{\theta \theta }\sqrt {1-\mu ^{2}})$) scale as $O(\bar {\zeta }_0^{2}\delta ^{-2})$ under the same conditions, since the normal stresses themselves are $O(\delta ^{-2})$ (see table 1), whereas $\tilde {\mathcal {T}}$ and $\tilde {{\mathcal {S}}} \sim O(\bar {\zeta }_0^{2})$ – see (3.5). Since $\bar {\zeta }_0\ll 1$ (weak surface charge), the linear (i.e. Newtonian) parts of the stresses dominate both inside and outside the EDL. From the above discussion, it immediately follows that, despite nominally $De$ being $O(1)$, the effective flow is only weakly viscoelastic, since the Newtonian, i.e. the linear, component of the stress–strain relation dominates. Below, we report the order-wise velocity components in the inner layer and the resulting slip velocity.

(i) The $O(\bar {\zeta }_0)$ velocity field: at $O(\bar {\zeta }_0)$, $\varPhi _1 = \bar {\zeta }(\theta ,\varphi ) e^{-R}$, while $C_1 = 0$ and $\varPi _1 = -2\varPhi _1$. This is the first approximation and essentially amounts to the Debye–Huckel linearization. The charge density, concentration and the potential can be found by solving the inner layer Nernst–Planck and Poisson equations, as outlined in (3.4*a*)–(3.4*b*), subject to the boundary conditions, (3.6*a*) and (3.6*b*). As already mentioned, at this order, the velocity field is identical to that of a Newtonian fluid, when nominally $De\sim O(1)$. The velocity field thus has the solution

*a*)\begin{gather} U_1 = \varGamma_1(\varphi)-\frac{3\beta\omega_1(\mu,\varphi)}{\sqrt{1-\mu^{2}}}(1-\textrm{e}^{{-}R}), \end{gather}

*b*)\begin{gather}W_1 = \chi_1(\mu,\varphi), \end{gather}

*c*)\begin{gather}V_1 = 3\beta\omega_{1,\mu}(\mu,\varphi)(1-R-\textrm{e}^{{-}R})-\omega_2(\mu,\varphi)R. \end{gather}

In the above, $\omega _1(\mu ,\varphi ) = \bar {\zeta }(\theta ,\varphi )Q_1(\mu )$, $\omega _2(\mu ,\varphi ) = {(\chi _{1,\varphi }+\mu \varGamma _1)}/{\sqrt {1-\mu ^{2}}}$ and $\omega _{1,\mu } = {\partial \omega _1}/{\partial \mu }$, $\omega _{1,\mu \mu } = {\partial ^{2}\omega _1}/{\partial \mu ^{2}}$, $\chi _{1,\varphi } = {\partial \chi _1}/{\partial \varphi }$ etc., while $Q_n(x)$ is the Gegenbauer polynomial of the first kind and order $n$ (Leal Reference Leal2007), expressed as $Q_n(x) = \int _{-1}^{x} P_n(u) \,\textrm {d}u$. For example (Leal Reference Leal2007), $Q_1(x) = (x^{2}-1)/{2}$, $Q_2(x) = xQ_1(x)$, etc. Further, we have defined $\varGamma _{k}(\varphi ) = \varOmega _y^{(k)}\cos (\varphi )-\varOmega _x^{(k)}\sin (\varphi )$ and $\chi _k(\mu ,\varphi ) = -(\varOmega _x^{(k)}\cos (\theta )+\varOmega _y^{(k)}\sin (\theta )) \sin (\varphi )$. The leading-order slip velocity without rotation is thus given by

(ii) The $O(\bar {\zeta }_0^{2})$ velocity field: at $O(\bar {\zeta }_0^{2})$, the nonlinear components of the viscoelastic stresses play a key role in altering the velocity field, through the convected derivatives. The leading-order terms in the expansion of the convected derivatives may be easily evaluated from their expressions given in Appendix A. Further note that, at $O(\bar {\zeta }_0^{2})$, $\varPhi _2 = \varPi _2 = 0$, $C_2 = \varPhi _1^{2}$ and from the $r$-momentum equation, $P_2 = \frac {1}{2}\bar {\zeta }^{2}e^{-2R}$. Combining (3.4*c*)–(3.4*f*), along with (3.5), it may be shown that, $U_2$ and $W_2$ are governed by the following equations:

*a*)\begin{gather} \frac{\partial^{2} U_2}{\partial R^{2}} = 2(\lambda_2-\lambda_1)De\left[\frac{\partial}{\partial\mu}\left\{\sqrt{1-\mu^{2}}\tilde{{\mathcal{S}}}_{\theta\theta}^{(2)}\right\}-\frac{\partial\tilde{{\mathcal{S}}}_{r\theta}^{(2)}}{\partial R}\right], \end{gather}

*b*)\begin{gather}\frac{\partial^{2} W_2}{\partial R^{2}} = (\lambda_1-\lambda_2)De \left\{\omega_3(\mu,\varphi)\frac{\partial^{2} U_1}{\partial R^{2}}\right\}. \end{gather}

In the above, $\omega _3 = \sqrt {1-\mu ^{2}}\chi _{1,\varphi }+{\mu \chi _1}/{\sqrt {1-\mu ^{2}}}$. Equation (3.15) may be solved for $U_2$ and $W_2$ subject to the no-slip condition at the particle surface and the constraint that both the velocity components remain bounded, to obtain

*a*)\begin{gather} U_2 = De(\lambda_1-\lambda_2)[\mathcal{A}_1(R,\mu)\omega_1^{2}+\mathcal{A}_2(R,\mu)\omega_1+\mathcal{A}_3(R,\mu)\omega_{1,\mu}+\mathcal{A}_4(R,\mu)\omega_{1,\varphi}]+\varGamma_2, \end{gather}

*b*)\begin{gather}W_2 = \chi_2-De(\lambda_1-\lambda_2)\omega_3(\mu,\varphi)\left(U_1-\varGamma_1\right). \end{gather}

Expressions for $\mathcal {A}_1,\mathcal {A}_2,\ldots$ etc. are included in Appendix B. The continuity equation may be used to obtain $V_2$

The $O(\bar {\zeta }_0^{2})$ contribution to the slip velocity is therefore,

*a*)\begin{gather} {\boldsymbol{v}}_{S}^{(2)} = \lim_{R\rightarrow\infty}[U_2{\hat{\boldsymbol{e}}}_{\theta} + W_2{\hat{\boldsymbol{e}}}_{\phi}] - \boldsymbol{\varOmega}^{(2)}\times{\hat{\boldsymbol{e}}}_r = v_{S,\theta}^{(2)}{\hat{\boldsymbol{e}}}_{\theta} + v_{S,\varphi}^{(2)}{\hat{\boldsymbol{e}}}_{\varphi},\end{gather}

*b*)\begin{gather} v_{S,\theta}^{(2)} = {\boldsymbol{v}}_{S}^{(2)}\boldsymbol{\cdot}{\hat{\boldsymbol{e}}}_{\theta} = De(\lambda_1-\lambda_2)\left[-\frac{9\mu\beta^{2}}{2\left(1-\mu^{2}\right)^{3/2}}\omega_1^{2}+\omega_1 \left\{\frac{3\beta\mu}{1-\mu^{2}}\varGamma_1-\frac{27\beta^{2}}{\sqrt{1-\mu^{2}}}\omega_{1,\mu}\right.\right.\nonumber\\ \quad \left.\left. -\frac{3\beta}{\sqrt{1-\mu^{2}}}\omega_2\right\}+3\beta\varGamma_1\omega_{1,\mu}+\frac{3\beta}{1-\mu^{2}}\chi_1\omega_{1,\varphi}\right], \end{gather}

*c*)\begin{gather} v_{S,\varphi}^{(2)} = {\boldsymbol{v}}_{S}^{(2)}\boldsymbol{\cdot}{\hat{\boldsymbol{e}}}_{\varphi} = De(\lambda_1-\lambda_2)\frac{3\beta\omega_1\omega_3}{\sqrt{1-\mu^{2}}}. \end{gather}

More discussion on the nature of the velocity profiles at $O(\bar {\zeta }_0^{2})$ have been included in § 3.5.

(iii) The $O(\bar {\zeta }_0^{3})$ velocity field: at $O(\bar {\zeta }_0^{3})$, $\varPhi _3 = \frac {1}{48}\bar {\zeta }^{3}(\theta ,\phi )(\textrm {e}^{-3R}-3\textrm {e}^{-R})$ and the charge density has the form, $\varPi _3 = -2\varPhi _3-\frac {1}{3}\varPhi _1^{3}$. The governing equations for the velocity components may be derived by inserting the $O(\bar {\zeta }_0^{3})$ stresses into the governing equations for the inner layer. It may be verified that the $O(\bar {\zeta }_0^{3})$ velocities satisfy equations of the following form:

*a*)\begin{gather} \frac{\partial^{2} U_3}{\partial R^{2}} = 2(\lambda_1-\lambda_2)De\tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\left[\left(\boldsymbol{\tilde{{\mathcal{S}}}} -\lambda_1 De\,\Im \right)\boldsymbol{\cdot}{\hat{\boldsymbol{e}}}_{\theta}\right]-\frac{3}{2}\beta\sqrt{1-\mu^{2}} \frac{\partial^{2}\varPhi_3}{\partial R^{3}}, \end{gather}

*b*)\begin{gather}\frac{\partial^{2} W_3}{\partial R^{2}} = 2(\lambda_1-\lambda_2)De\tilde{\boldsymbol{\nabla}}_{*}\boldsymbol{\cdot}\left[\left(\boldsymbol{\tilde{{\mathcal{S}}}} -\lambda_1 De\,\Im \right)\boldsymbol{\cdot}{\hat{\boldsymbol{e}}}_{\varphi}\right], \end{gather}

where, $\tilde {\boldsymbol {\nabla }} = {\hat {\boldsymbol {e}}}_r({\partial }/{\partial R}) -{\hat {\boldsymbol {e}}}_{\theta }({\partial }/{\partial \mu })\sqrt {1-\mu ^{2}}+{\hat {\boldsymbol {e}}}_{\varphi }({\partial }/{\partial \varphi })$ and $\tilde {\boldsymbol {\nabla }}_{*} = {\hat {\boldsymbol {e}}}_r({\partial }/{\partial R})-{\hat {\boldsymbol {e}}}_{\theta }({\partial }/{\partial \mu}) \sqrt {1-\mu ^{2}}+{\hat {\boldsymbol {e}}}_{\varphi }({\mu }/{\sqrt {1-\mu ^{2}}})$. In the above, components of $\Im $ can be deduced from (3.5) in combination with (A1) and (A2) in Appendix A. Detailed expressions for various components of $\Im $ have been included in the supplementary material – see § S1.4 therein.

The solutions for $U_3$ and $W_3$ may be computed by integrating the above equations twice, subject to the boundary conditions, $U_3 = \varGamma _3$ and $W_3 = \chi _3$ at $R = 0$ and both remain bounded as $R\rightarrow \infty$. Although the process is straightforward, the results are rather cumbersome to represent and hence we do not give them here; they will be made available upon request to the authors. That said, it is worth noting that the $O(\bar {\zeta }_0^{3})$ contribution to the slip velocity without particle rotation, takes the following form:

*a*)\begin{gather} {\boldsymbol{v}}_{S}^{(3)} = \lim_{R\rightarrow\infty}[U_3{\hat{\boldsymbol{e}}}_{\theta} + W_3{\hat{\boldsymbol{e}}}_{\phi}] = v_{S,\theta}^{(3)}{\hat{\boldsymbol{e}}}_{\theta}, \end{gather}

*b*)\begin{gather} v_{S,\theta}^{(3)} = \frac{De^{2}\beta^{3}(\lambda_1-\lambda_2)}{\sqrt{1-\mu^{2}}}\left[27(11\lambda_1-18\lambda_2)\omega_1\omega_{1,\mu}^{2}+ \frac{81(4\lambda_1-5\lambda_2)}{2(1-\mu^{2})}\omega_1^{2}\omega_{1,\mu}\right.\nonumber\\ \left.\quad + \frac{27}{4}(21\lambda_1-34\lambda_2)\omega_1^{2}\omega_{1,\mu\mu}+ \frac{3(\lambda_1-\lambda_2)(37\mu^{2}+29)}{2(1-\mu^{2})^{2}}\omega_1^{3}\right] +\frac{\beta\bar{\zeta}^{3} Q_1}{8\sqrt{1-\mu^{2}}}. \end{gather}

The complete modified Smoluchowski slip velocity at $O(\bar {\zeta }_0^{3})$ accounting for particle rotation, has been included in the supplementary material – see § S1.4 therein. The modified Smoluchowski slip for a viscoelastic fluid in the presence of arbitrary surface charge is thus given by (up to $O(\bar {\zeta }_0^{3})$)

The contributions at the individual orders of $\bar {\zeta }_0$ are given in (3.14), (3.18) and (3.20).

### 3.5. Effect of viscoelasticity on the Smoluchowski slip: the key features

There are several interesting points to note from the modified slip velocity derived above. First, recall that the regular expansion to incorporate the effects of viscoelastic stresses is in $\bar {\zeta }_0$, instead of $De$, which in many cases is the usual choice (Bird *et al.* Reference Bird, Armstrong and Hassager1987). In this regard, it should be noted here that, although nominally $De\sim O(1)$, because the velocity is $O(\bar {\zeta }_0)$ (on account of weak surface charge), the effective Deborah number becomes $De_{eff} = \bar {\zeta }_0 De \ll 1$. Therefore, the regular expansion in $\bar {\zeta }_0$ may also be treated as a regular expansion in $De_{eff}$. Note that the fluid actually ‘sees’ the effective Deborah number ($De_{eff}$) to manifest the interplay of the electro-mechanics and viscoelastic hydrodynamics and not the nominal one ($De$) and hence the overall flow here is only weakly viscoelastic in nature (as also stated in § 2.2). As a consequence, the expansion in (3.21) is exactly equivalent to an ordered expansion around the Newtonian limit, carried out for an Oldroyd-B fluid. For example, the $O(\bar {\zeta }_0^{2})$ solution is equivalent to the $O(De_{eff})$ correction in an ordered-fluid expansion. This indicates that one does need to further impose the limit of small polymer concentration ($\mathcal {C}\ll 1$), which translates into $\lambda _1/\lambda _2-1\ll 1$, to derive closed-form analytical solutions. Later, in § 4.3.2, we explore the limit of low polymer concentration to compare our solutions with previously reported studies in the literature (Li & Koch Reference Li and Koch2020).

The basic physics behind the $O(\bar {\zeta }_0^{2})$ and the $O(\bar {\zeta }_0^{3})$ equations may be appreciated as follows. The leading-order flow (at $O(\bar {\zeta }_0)$) is effectively Newtonian, which stretches the polymers in the EDL, thus giving rise to excess polymeric stresses at $O(\bar {\zeta }_0^{2})$. These excess polymeric stresses are non-uniform along the particle surface, because of its curvature and variations in $\bar {\zeta }(\theta ,\varphi )$. As a result, the gradients of these excess stresses drive a flow at $O(\bar {\zeta }_0^{2})$, which is observed in (3.15). Because of the nonlinear constitutive relation of the fluid itself, this $O(\bar {\zeta }_0^{2})$ velocity field and the leading-order Newtonian contribution interact between each other and give rise to additional polymeric stretching, which results in higher-order stresses, as reflected in the $O(\bar {\zeta }_0^{3})$ equations in (3.19).

From the discussion in § 3.2 following table 1, it is evident that, inside the EDL, the polymeric stresses play a major role in governing the local flow patterns. This can be better understood by defining a characteristic Deborah number inside the EDL as $De_{EDL} = u_c\lambda _0/\lambda _D = \delta ^{-1}De$, since the effective length scale inside the EDL is $O(\delta )$. If $De\sim O(1)$, it follows that $De_{EDL}\gg 1$, which indicates that, for $U\sim O(1)$, the polymeric stresses play a key role in governing the motion inside the EDL. This is exactly equivalent to the fact that the stress components $\tau _{\theta \theta }$, $\tau _{\varphi \varphi }$ and $\tau _{\theta \varphi }$ all scale as $\delta ^{-2}$ inside the EDL, which are also signatures of strong polymeric stretching therein. This assertion becomes clear from a detailed rescaling of the constitutive equations inside the EDL, which has been provided in § S1.3 of the supplementary material. The same may also understood from (3.5), which for the $r\theta$ component (as an example) may be rewritten as

where $\tilde {\mathcal {T}}_{r\theta }$ and $\tilde {\mathcal {S}}_{r\theta }$ are given in Appendix A and they denote the effects of polymeric stresses inside the EDL. Since $De_{EDL}\sim \delta ^{-1} \gg 1$, $\delta De_{EDL} \sim O(1)$ and hence these terms cannot be neglected, although they are multiplied with $\delta$. In the outer region (addressed in the next section), the Oldroyd-B constitutive model should be applicable as long as the effective Deborah number is small ($De_{eff} \ll 1$). As a result, it may be inferred that the analysis in § 3.4 is indeed valid when $De\sim O(1)$, which would imply $De_{eff} = \bar {\zeta }_0 De \ll 1$ in the outer region, for weakly charged particles. In other words, in the outer region, effectively the Deborah number is small and hence, despite $De$ (the nominal Deborah number) being $O(1)$, the Oldroyd-B constitutive relation should apply. On the other hand, when the motion inside the EDL is considered, our analysis is valid when $\delta De_{EDL} \sim O(1)$, as indicated in (3.22) above.

Further note from (3.14), (3.18) and (3.20) that the Smoluchowski slip depends on the quantity $\omega _1 = \bar {\zeta }(\theta ,\varphi )Q_1(\mu )$ and its derivatives, which encodes information about variations in surface charge density as well as the particle's curvature. While $\omega _1$ itself determines the Smoluchowski slip at the leading order, its derivatives and their products appear in the higher-order corrections, as may be observed in (3.18) and (3.20). The derivatives of $\omega _1$ contain $\bar {\zeta }(\theta ,\varphi )$ and its derivatives, which underscore the effects of variations in surface charge density, while the factor $Q_1(\mu )/\sqrt {1-\mu ^{2}}$ bears the signature of the particle's curvature. Notice that, even for $\bar {\zeta } =$ a constant, i.e. for a uniformly charged particle, the derivatives of $\omega _1$ are in general non-zero because of $Q_1(\mu )$ appearing in them, which indicates that the curvature will affect the slip velocity, appearing in the $O(\bar {\zeta }_0^{2})$ and $O(\bar {\zeta }_0^{3})$ terms. Recall that the effect of the particle's curvature influences the velocity inside the EDL because of the anomalous scaling shown by the normal stresses, as discussed in § 3.2. As a consequence, the modified Smoluchowski slip would depend on the particle's radius ($a$) in a viscoelastic fluid. All of these features are in stark contrast to what is observed in Newtonian fluids (Ajdari Reference Ajdari1995; Yariv Reference Yariv2009). Although it is difficult to understand this facet from the general expression given in (3.21) for an arbitrary surface charge, the explicit effect of particle curvature becomes clear when one considers a uniformly charged particle, as done in the next section.

Further, note that the linear relation observed between the applied external fields and the Smoluchowski slip for Newtonian media gets lost in viscoelastic fluids, because of the appearance of nonlinear terms in $\beta$, as evident from (3.18) and (3.20). The third important point to note is that, for $\lambda _1 = \lambda _2$, one recovers the slip velocity for a Newtonian fluid – this is of course an intrinsic property of the Oldroyd-B constitutive relation (Bird *et al.* Reference Bird, Armstrong and Hassager1987). The same limit may also be recovered by enforcing $De = 0$.

Finally, we shall quickly summarize the most interesting outcome of particle rotation on the modified Smoluchowski slip, as noted in (3.18) and § S1.4 of the supplementary material, where the $O(\bar {\zeta }_0^{3})$ contributions have been reported. Notice from (3.18) that, in presence of particle rotation (when $\omega _3 \neq 0$), the slip velocity at $O(\bar {\zeta }_0^{2})$ also has a component along the ${\hat {\boldsymbol {e}}}_{\varphi }$ direction, resulting in anisotropic motion. When the medium is Newtonian, i.e. $De = 0$, this component vanishes. We therefore conclude that the slip velocity (${\boldsymbol {v}}_S$) in a viscoelastic medium does depend on the particle's angular velocity; this is again in stark contrast to Newtonian fluids. Because of the nonlinear rheology of the fluid combined with the fact that rotational motion leads to differential velocity at various points on the surface, the angular velocity of the particle gets embedded in the motion within the EDL, from $O(\bar {\zeta }_0^{2})$ onwards.

Although the expression of the Smoluchowski slip velocity derived above applies to an arbitrary distribution of $\bar {\zeta }$, we shall now look into two specific examples of non-uniformly charged particles and apply the analysis carried out herein to derive their electrophoretic motion. In the first instance, we shall consider pure translation of a particle, carrying non-uniform but axisymmetric charge. The second example will explore pure electrophoretic rotation of a particle with non-axisymmetric surface charge density.

### 3.6. On the choice of the constitutive model

In our analysis, the Oldroyd-B constitutive model has been used, considering a specific vision. This model does not suffer from many of the limitations of its predecessors (e.g. the second-order model only applies to small strain rates). At the same time, it is able to capture many critical behaviours exhibited by polymeric fluids, which include, for instance, growth of shear stress at low to moderate shear rates; it also correctly predicts that the first normal stress difference varies as $\dot {\gamma }^{2}$ ($\dot {\gamma }$ is the shear strain rate) in polymeric liquids. In addition, some of the key physics of confluence between non-uniformity in the surface charge distribution and the viscoelastic rheology may indeed be probed by appealing to this model, as evidenced by the results presented later.

That said, the Oldroyd-B model is not without its limitations (Bird *et al.* Reference Bird, Armstrong and Hassager1987). Perhaps two of its biggest shortcomings are that it fails to account for shear dependent viscosity and normal stress coefficients, exhibited by polymeric liquids. Moreover, the Oldroyd-B model also predicts infinite extensional viscosity, beyond a critical rate of strain. Therefore, this model should be used with caution when the flow is strongly extensional in nature. It is generally accepted that, for flows with large elongational stresses, use of the Oldroyd-B model requires special care to ensure that the elongational viscosity remains finite. For further insight, one may refer to the book of Bird *et al.* (Reference Bird, Armstrong and Hassager1987).

In view of the above, one may therefore conclude that the strong stretching rates inside the EDL call into question the accuracy of the Oldroyd-B model used here. As a result, it is perhaps judicious to compare its predictions with those from more robust nonlinear viscoelastic models (Bird *et al.* Reference Bird, Armstrong and Hassager1987; Afonso *et al.* Reference Afonso, Alves and Pinho2009) that remain valid for strong polymer stretching. To this end, in § 4.4 we have carried out numerical simulations for electrophoretic motion of a non-uniformly (but axisymmetric) charged particle in a FENE-P fluid (Afonso *et al.* Reference Afonso, Alves and Pinho2009). This model is one of the simplest nonlinear viscoelastic models that is able to physically account for shear thinning in polymeric liquids. In the appropriate limiting case (see § 4.4 for further details), the numerical solutions have been compared with the analytical ones from the Oldroyd-B model. Overall, we observe that the analytical predictions for the electrophoretic velocity using the Oldroyd-B model remain reasonably accurate, when the surface charge is sufficiently small and shear thinning effects are subdominant. However, in spite of its well-known limitations, the importance of Oldroyd-B model in developing insights into the dynamics of complex fluids remains undisputed, as evidenced from the vast body of existing literature (see for instance 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) premised on this model. As a consequence, it is perhaps justified to assert that the predictions based on the Oldroyd-B model can indeed capture much of the essential physics of some of the complex fluids.

## 4. Case study-I: electrophoretic translation of a non-uniformly charged particle

### 4.1. Overview of the analysis

As the first representative example, we take up the case of a spherical particle carrying non-uniform but axisymmetric surface charge and derive an expression of its electrophoretic mobility in the thin EDL limit. This particular example has been chosen for its relatively simpler demand of algebraic manipulation as well as the light it sheds on the essential physics governing electrophoresis in a viscoelastic medium even without bringing the particle's rotation into purview. A generic axisymmetric but otherwise non-uniform charge density on a particle may be expressed as $\bar {\zeta }(\mu ,\varphi ) = \sum _{n=0}^{\infty } a_n P_n(\mu )$. Here, to keep the algebra tractable, we choose $a_0 \neq 0$ and $a_1 \neq 0$, while $a_n = 0$, $\forall n\geqslant 2$; hence, the charge density is given by $\bar {\zeta }(\mu ,\varphi ) = a_0 P_0(\mu ) + a_1 P_1(\mu ) = a_0 + a_1\mu$. It may be checked that, in this case, the particle carries a net charge (appropriately non-dimensionalized) amounting to $4{\rm \pi} a_0$. Since we are only considering axisymmetric flow, we may straight away make a couple of assertions as follows: (i) the particle will not rotate, i.e. $\boldsymbol {\varOmega } = 0$ and (ii) the particle will translate along the $z$-axis, i.e. ${\hat {\boldsymbol {e}}}_u = {\hat {\boldsymbol {e}}}_z$.

It may be easily verified from (3.14), (3.18) and (3.20) that, for a particle with axisymmetric charge density, the ‘modified Smoluchowski slip’ takes the following form:

*a*,

*b*)\begin{equation} {\boldsymbol{v}}_S = v_{S,\theta}{\hat{\boldsymbol{e}}}_{\theta},\quad\text{and}\quad v_{S,\theta} = \bar{\zeta}_0 v_{S,\theta}^{(1)} + \bar{\zeta}_0^{2} v_{S,\theta}^{(2)} + \bar{\zeta}_0^{3} v_{S,\theta}^{(3)} + \ldots. \end{equation}

Inserting $\bar {\zeta } = a_0 + a_1\mu$ into (4.1*a*,*b*), we deduce that the slip velocity for this particular case looks as follows:

*a*)\begin{gather} v_{S,\theta}^{(1)} ={-}\frac{3\beta}{\sqrt{1-\mu^{2}}}\left(a_0 Q_1 + a_1 Q_2\right) \end{gather}

*b*)\begin{gather}v_{S,\theta}^{(2)} = \frac{De\beta^{2}(\lambda_1-\lambda_2)}{\sqrt{1-\mu^{2}}} \left[\frac{9}{10}a_0 a_1 Q_1 - \frac{9(77a_0^{2} + 9 a_1^{2})}{28} Q_2-\frac{252}{5}a_0 a_1 Q_3 - \frac{153}{7} a_1^{2} Q_4\right] \end{gather}

*c*)\begin{gather} v_{S,\theta}^{(3)} = \left[\frac{1}{8}a_0^{3} \beta + \frac{3}{40}\beta a_0 a_1^{2} + De^{2}\beta^{3} (\lambda_1-\lambda_2)\left\{-\frac{3}{20}(\lambda_1+8\lambda_2)a_0^{3}\right.\right.\nonumber\\ \quad\left.\left. -\frac{3}{20}\left(\frac{57}{7}\lambda_1-\frac{48}{7}\lambda_2\right)a_0 a_1^{2}\right\}\right] \frac{Q_1}{\sqrt{1-\mu^{2}}} + \text{h.o.t in}\ Q_n\text{'s}. \end{gather}

In (4.2*c*), components up to $Q_6$ contribute to the modified Smoluchowski slip, although we only require the term associated with $Q_1$ to compute the electrophoretic mobility at this order. Therefore, we do not mention the other higher-order terms here and their mathematical expressions will be made available upon request. Note that the special case of a uniformly charged particle may be recovered by inserting $a_0 = 1$ and $a_1 = 0$ (i.e. $\bar {\zeta } = 1$). It may be easily verified that, in such a case,

The components $v_{S,\theta }^{(1)}$ and $v_{S,\theta }^{(2)}$ may also be evaluated by inserting $a_0 = 1$ and $a_1 = 0$ into (4.2*a*) and (4.2*b*), respectively. Here, we would like to emphasize that the special case of a uniformly charged particle explicitly filters out the influence of its curvature on the liquid motion as well as on its electrophoretic mobility. Since the derivatives of the surface charge with respect to the polar angle ($\mu = \cos \theta$) vanish, the only alterations to the slip velocity comes from the curvature of the particle, as denoted by the factor $Q_1(\mu )/\sqrt {1-\mu ^{2}}$.

The electrophoretic velocity $\mathcal {U}$ may be ascertained by solving the governing equations in the outer layer, i.e. (3.10*a*) and (3.10*b*), subject to the far field condition (3.11*a*) and the matching condition (3.11*b*) at the particle surface, wherein ${\boldsymbol {v}}_S$ (${=}v_{S,\theta }{\hat {\boldsymbol {e}}}_{\theta }$) is given by (4.2) (also see the discussion after (3.8)). The final step is to apply the force balance at the edge of the EDL, as outlined in (2.4*a*,*b*). Here, the force balance effectively reduces to the following form (Leal Reference Leal2007): $\mathbb {F}_z = 2{\rm \pi} \int _{0}^{{\rm \pi} } [(-p+\tau _{rr})\mu -\sqrt {1-\mu ^{2}}\tau _{r\theta }]_{r=1} \,\textrm {d}\mu = 0$. Following the analysis for the EDL, the outer layer variables are also expanded in a regular asymptotic series of $\bar {\zeta }_0$. Therefore, the electrophoretic velocity and the net force on the particle are expanded as

*a*)\begin{gather} \mathcal{U} = \bar{\zeta}_0\mathcal{U}_1 + \bar{\zeta}_0^{2}\mathcal{U}_2 + \bar{\zeta}_0^{3}\mathcal{U}_3 + \ldots \end{gather}

*b*)\begin{gather}\mathbb{F}_z = \bar{\zeta}_0\mathbb{F}_z^{(1)} + \bar{\zeta}_0^{2}\mathbb{F}_z^{(2)} + \bar{\zeta}_0^{3}\mathbb{F}_z^{(3)} + \ldots\,, \end{gather}

where,

*c*)\begin{equation} \mathbb{F}_z^{(k)} = 2{\rm \pi}\int_{0}^{\rm \pi} \left[({-}p_k+\tau_{rr}^{(k)})\mu-\sqrt{1-\mu^{2}} \tau_{r\theta}^{(k)}\right]_{r=1} \,\textrm{d}\mu = 0. \end{equation}

All other variables may also be expanded in a similar way. The $k$th-order (in $\bar {\zeta }_0$) outer layer equations may be expressed as (for $k = 1, 2,3,\ldots$ etc.)

*a*)\begin{gather} -\boldsymbol{\nabla} p_k + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tau}^{(k)} = 0\quad \text{and}\quad \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}_k = 0, \end{gather}

*b*)\begin{gather}\boldsymbol{\tau}^{(k)} + De\lambda_1 \boldsymbol{\mathcal{T}}^{(k)} = 2{\boldsymbol{D}}_{k} + 2\lambda_2 De \boldsymbol{{\mathcal{S}}}^{(k)}. \end{gather}

These are subject to the following boundary conditions:

*a*)$$\begin{gather} \text{at}\ r\rightarrow \infty,\quad {\boldsymbol{v}}_k ={-}\mathcal{U}_k{\hat{\boldsymbol{e}}}_z, \end{gather}$$

*b*)$$\begin{gather}\text{at}\ r = 1,\quad u_{r}^{(k)} = 0;\quad u_{\theta}^{(k)} = v_{S,\theta}^{(k)}. \end{gather}$$

The force balance at every order of $\bar {\zeta }_0$ reads $\mathbb {F}_z^{(k)} = 0$, which may be used to derive $\mathcal {U}_k$. Note that, using (4.5*b*), the momentum equation (4.5*a*) may be rewritten as follows: $-\boldsymbol {\nabla } p_k + \boldsymbol {\nabla }\boldsymbol {\cdot }[\boldsymbol {\nabla }{\boldsymbol {v}}_k+(\boldsymbol {\nabla }{\boldsymbol {v}}_k)^\textrm {T}] + De\boldsymbol {\nabla }\boldsymbol {\cdot }(2\lambda _2\boldsymbol {{\mathcal {S}}}^{(k)}-\lambda _1\boldsymbol {\mathcal {T}}^{(k)}) = 0$ and $\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {v}}_k = 0$. For an axisymmetric flow, (4.5*a*) may be solved using a streamfunction (Leal Reference Leal2007), defined as (for $k$th order of $\bar {\zeta }_0$) $u_r^{(k)} = -({1}/{r^{2}})({\partial \varPsi _k}/{\partial \mu })$ and $u_{\theta }^{(k)} = -({1}/{r\sqrt {1-\mu ^{2}}})({\partial \varPsi _k}/{\partial r})$, where $\varPsi$ is the streamfunction. As such, the momentum equation at any order of $\bar {\zeta }_0$ may be expressed in terms of $\varPsi$ as follows (Leal Reference Leal2007; Goswami *et al.* Reference Goswami, Dhar, Ghosh and Chakraborty2017):

where (Leal Reference Leal2007), $\mathcal {E}^{2} = {\partial ^{2}}/{\partial r^{2}} + (({1-\mu ^{2}})/{r^{2}})({\partial ^{2}}/{\partial \mu ^{2}})$; $f_{\theta }^{(k)} = De\boldsymbol {\nabla }\boldsymbol {\cdot }(2\lambda _2\boldsymbol {{\mathcal {S}}}^{(k)}-\lambda _1\boldsymbol {\mathcal {T}}^{(k)})\boldsymbol {\cdot }{\hat {\boldsymbol {e}}}_{\theta }$ and $f_{r}^{(k)} = De\boldsymbol {\nabla }\boldsymbol {\cdot }(2\lambda _2\boldsymbol {{\mathcal {S}}}^{(k)}-\lambda _1\boldsymbol {\mathcal {T}}^{(k)})\boldsymbol {\cdot }{\hat {\boldsymbol {e}}}_{r}$.

### 4.2. Solution for the streamfunction and the particle velocity

Following the outline of § 4.1, it is straightforward to write the solutions to the streamfunction, apply the force balance and subsequently find the electrophoretic velocity at successive orders of $\bar {\zeta }_0$. At $O(\bar {\zeta }_0)$ the flow behaves exactly the same as that of a Newtonian fluid. The streamfunction satisfies the equation $\mathcal {E}^{4}\varPsi _1 = 0$, subject to (4.6*a*) at the far field and (4.6*b*) at the particle surface (with $k = 1$), where $v_{S,\theta }^{(1)}$ is given by (4.2*a*). The resulting solutions for the streamfunction and the particle velocity read

*a*)\begin{gather} \varPsi_1 = \beta a_0 \left(r^{2}-\frac{1}{r}\right)Q_1(\mu) + \frac{3}{2}\beta a_1\left(1-\frac{1}{r^{2}}\right)Q_2(\mu), \end{gather}

*b*)\begin{gather}\mathcal{U}_1 = \beta a_0. \end{gather}

At $O(\bar {\zeta }_0^{2})$, the streamfunction satisfies the equation $\mathcal {E}^{4}\varPsi _2 = \beta ^{2} De(\lambda _1-\lambda _2)\sum _{n=1}^{4}\mathcal {N}_n^{(2)} (r) Q_n(\mu )$, subject to (4.6*a*) and (4.6*b*), where $v_{S,\theta }^{(2)}$ is given by (4.2*b*). Here, $\mathcal {N}_1^{(2)} = -{216 a_0 a_1}/{r^{8}}$, $\mathcal {N}_2^{(2)} = {1296 a_1^{2}}/{7r^{9}}(r^{2}-{25}/{6})$, $\mathcal {N}_3^{(2)} = -{648 a_0 a_1}/{r^{8}}$ and $\mathcal {N}_4^{(2)} = 3240 a_1^{2}/{7 r^{9}}(r^{2}-3)$. Carrying out the force balance yields the following solutions for the electrophoretic velocity and the streamfunction at $O(\bar {\zeta }_0^{2})$:

*a*)\begin{gather} \varPsi_2 = \beta^{2} De (\lambda_1-\lambda_2)\sum_{n=1}^{4} \mathcal{M}_n^{(2)}(r) Q_n(\mu), \end{gather}

*b*)\begin{gather}\mathcal{U}_2 ={-}\frac{3}{5} \beta^{2} De (\lambda_1-\lambda_2) a_0 a_1 . \end{gather}

Expressions for $\mathcal {M}^{(2)}_k$ values ($k = 1,2,\ldots$) are included in Appendix B. The analysis at $O(\bar {\zeta }_0^{3})$ also follows the same pattern as that of the previous orders. It may be shown that the streamfunction, $\varPsi _3$ satisfies the following equation: $\mathcal {E}^{4}\varPsi _3 = De^{2}\beta ^{3} (\lambda _1-\lambda _2)\sum _{n=1}^{6}\mathcal {N}_n^{(3)}(r) Q_n(\mu )$. As already mentioned, to compute the electrophoretic mobility, only the contribution from the mode $n = 1$ will be adequate. As such, the expressions for $\mathcal {N}^{(3)}_1$ have been included in the supplementary material (see § S2.1 therein). The streamfunction has the form $\varPsi _3 = \sum _{n=1}^{6} \mathcal {M}_n^{(3)}(r) Q_n(\mu )$; expressions are given in Appendix-B. Carrying out the force balance yields the following expressions for the electrophoretic mobility:

*a*)\begin{align} \mathcal{U}_3 \!&=\!{-}\frac{20\,819}{5720}\beta^{3} De^{2} (\lambda_1\!-\!\lambda_2)\left[\left(\lambda_1\!-\!\frac{16\,445}{20\,819} \lambda_2\right) a_0^{3} \!+\! \frac{55\,749}{104\,095}\left(\lambda_1\!-\!\frac{41\,101}{130\,081}\lambda_2\right)a_0 a_1^{2}\right]\nonumber\\ &\quad -\beta a_0\left(\frac{a_0^{2}}{24}+\frac{a_1^{2}}{40}\right). \end{align}

Recall that the total electrophoretic mobility is given by $\mathcal {U} = \bar {\zeta }_0 \mathcal {U}_1 + \bar {\zeta }_0^{2} \mathcal {U}_2 + \bar {\zeta }_0^{3} \mathcal {U}_3 + \ldots\,$, where $\mathcal {U}_1$, $\mathcal {U}_2$ and $\mathcal {U}_3$ are given in (4.8*b*), (4.9*b*) and (4.10), respectively. More insight may be gained from the dimensional form of the electrophoretic velocity, expressed as

*a*)\begin{equation} \mathcal{U}' = \bar{\zeta}_0 \alpha + \bar{\zeta}_0^{2} \alpha^{2} \left(\frac{\lambda_0}{a}\right)\mathcal{G}_1 + \bar{\zeta}_0^{3}\left\{ \alpha^{3} \left(\frac{\lambda_0}{a}\right)^{2} \mathcal{G}_2 - \alpha a_0 \left(\frac{a_0^{2}}{24}+\frac{a_1^{2}}{40} \right)\right\} + O(\bar{\zeta}_0^{4}),\end{equation}

where,

*b*)\begin{gather} \alpha = \frac{\epsilon}{\eta}\left(\frac{kT}{e}\right) E_0,\quad\text{and}\quad \mathcal{G}_1 ={-}\frac{3}{5}(\lambda_1-\lambda_2)a_0 a_1, \end{gather}

*c*)\begin{gather}\mathcal{G}_2 ={-}\frac{20\,819}{5720}(\lambda_1-\lambda_2)\left[\left(\lambda_1-\frac{16\,445}{20\,819}\lambda_2\right)a_0^{3} + \frac{55\,749}{104\,095}\left(\lambda_1-\frac{41\,101}{130\,081}\lambda_2\right)a_0 a_1^{2}\right]. \end{gather}

For the special case of a uniformly charged particle (ucp), $a_0 = 1$ and $a_1 = 0$ and hence $\mathcal {G}_1 = 0$ and $\mathcal {G}_2 = -(\lambda _1-\lambda _2)(\frac {20\,819}{5720}\lambda _1-\frac {23}{8})\lambda _2)$ and thus

The complete solution for the $O(\bar {\zeta }_0^{3})$ streamfunction for a uniformly charged particle has been included in the supplementary material (see § S2 therein). Further discussion is included in the next subsection.

### 4.3. Comparison with previous studies

We shall first compare our results for the electrophoretic velocity ($\mathcal {U}'$) with two of the previous studies: (i) non-uniformly charged particle in a Newtonian medium, carried out by Anderson and coworkers (Anderson Reference Anderson1985; Fair & Anderson Reference Fair and Anderson1989) and (ii) uniformly charged particle in a viscoelastic medium with weak polymeric viscosity, investigated by Li & Koch (Reference Li and Koch2020). This is followed by a comparison with numerical simulations using the FENE-P constitutive model in § 4.4, where the accuracy of the Oldroyd-B model itself is assessed.

#### 4.3.1. Comparison with results for Newtonian fluids

Anderson and coworkers (Anderson Reference Anderson1985; Fair & Anderson Reference Fair and Anderson1989) have investigated electrophoresis of non-uniformly charged spherical particles in Newtonian fluids and reported general analytical solutions for the mobility for arbitrary ‘zeta’ potentials in the thin EDL limit. We shall use a different notation to represent their results herein. The zeta potential of the particle is denoted as $\phi '_P$ (non-dim. $\phi _P = \phi '_P/\psi _c$), while the electrophoretic velocity (with units) in Newtonian fluids is denoted by $\mathcal{U}'_N$. Anderson (Reference Anderson1985) showed that the electrophoretic velocity for an arbitrary distribution of $\phi '_P$ on the particle surface is given by

where $\boldsymbol {E}_{\infty }$ is the far-field imposed electric field (uniform), $\boldsymbol {H}'_2 = \langle \phi '_P\boldsymbol {Y}_2\rangle$ is the quadrupole moment of $\phi '_P$ and $\boldsymbol {Y}_2 = 3\hat {{\boldsymbol {n}}}\hat {{\boldsymbol {n}}}-\boldsymbol {I}$ is the second spherical harmonic. Also, $\langle \boldsymbol {\cdot }\rangle$ denotes the average over the particle surface. Referring to the description in § 4.1, we note that, here, $\boldsymbol {E}_{\infty } = E_0{\hat {\boldsymbol {e}}}_z$, ${\mathcal {U}}'_N = \mathcal {U}'_N{\hat {\boldsymbol {e}}}_z$, $\hat {{\boldsymbol {n}}} = {\hat {\boldsymbol {e}}}_r$ and $\phi '_P$ is related to the dimensionless surface charge density ($\sigma = \bar {\zeta }_0\bar {\zeta }$) as $\phi '_P = 2\psi _c\sinh ^{-1}(\frac {1}{2}\bar {\zeta }_0\bar {\zeta }(\mu ))$. As a result, for $\bar {\zeta }_0\ll 1$, $\phi '_P$ has the expansion, $\phi '_P/\psi _c = \bar {\zeta }_0\bar {\zeta }-\frac {1}{24}\bar {\zeta }_0^{3}\bar {\zeta }^{3} + \ldots\,$, where $\bar {\zeta } = a_0+a_1\mu$. Equation (4.13) may then be simplified to $\mathcal {U}'_N = {3\varepsilon E_0}/{2\eta }[\langle \phi '_P\rangle \boldsymbol {I}-\langle {\hat {\boldsymbol {e}}}_r{\hat {\boldsymbol {e}}}_r\phi '_P\rangle ]:{\hat {\boldsymbol {e}}}_z{\hat {\boldsymbol {e}}}_z$ and hence $\mathcal {U}'_N$ may be expanded as $\mathcal {U}'_N = \bar {\zeta }_0\mathcal {U}'^{(1)}_N+\bar {\zeta }_0^{2}\mathcal {U}'^{(2)}_N + \bar {\zeta }_0^{3}\mathcal {U}'^{(3)}_N + \ldots\,$, when $\bar {\zeta }_0\ll 1$. It can be verified that $\mathcal {U}'^{(1)}_N = {3\varepsilon \psi _c E_0}/{2\eta }[\langle \bar {\zeta }\rangle \boldsymbol {I}-\langle {\hat {\boldsymbol {e}}}_r{\hat {\boldsymbol {e}}}_r\bar {\zeta }\rangle ]:{\hat {\boldsymbol {e}}}_z{\hat {\boldsymbol {e}}}_z$, $\mathcal {U}'^{(2)}_N = 0$ and $\mathcal {U}'^{(3)}_N = {\varepsilon \psi _c E_0}/{16\eta }[-\langle \bar {\zeta }^{3}\rangle \boldsymbol {I}+\langle {\hat {\boldsymbol {e}}}_r{\hat {\boldsymbol {e}}}_r\bar {\zeta }^{3}\rangle ]:{\hat {\boldsymbol {e}}}_z{\hat {\boldsymbol {e}}}_z$. Taking $\bar {\zeta } = a_0+a_1\mu$, we obtain $\langle \zeta \rangle = a_0$ and $\langle {\hat {\boldsymbol {e}}}_r{\hat {\boldsymbol {e}}}_r\zeta \rangle = \frac {1}{3}a_0\boldsymbol {I}$ and hence, $\mathcal {U}'^{(1)}_N = \alpha a_0$. At the same time, $\langle \bar {\zeta }^{3}\rangle = a_0(a_0^{2}+a_1^{2})$ and $\langle {\hat {\boldsymbol {e}}}_r{\hat {\boldsymbol {e}}}_r\bar {\zeta }^{3}\rangle = a_0({a_0^{2}}/{3}+{a_1^{2}}/{5})\boldsymbol {I}+\frac {2}{5}a_0a_1^{2}{\hat {\boldsymbol {e}}}_z{\hat {\boldsymbol {e}}}_z$, which implies, $\mathcal {U}'^{(3)}_N = -\alpha a_0({a_0^{2}}/{24}+{a_1^{2}}/{40})$. Combining the $O(\bar {\zeta }_0)$ and $O(\bar {\zeta }_0^{3})$ results, Anderson's analysis yields the following electrophoretic velocity for the given surface charge distribution:

It may be noted that, by inserting $\lambda _0 = 0$ (Newtonian limit) in (4.11*a*), the velocity reported in (4.14) above is exactly recovered. We therefore conclude that our analysis can correctly reproduce the mobility of non-uniformly charged particles in the Newtonian limit.

#### 4.3.2. Comparison with results for viscoelastic fluids

Li & Koch (Reference Li and Koch2020) have analysed the electrophoretic mobility of a uniformly charged sphere in a Giesekus fluid. They considered the limit of small polymeric viscosity ($\mathcal {C} \ll 1$, defined later), weak surface charge ($\bar {\zeta }_0 \ll 1$) and thin EDL. Note that we have used different symbols (to those used by Li & Koch) to express their results. Li & Koch (Reference Li and Koch2020) derive their results for $\mathcal {C} = \eta ^{P}/\eta ^{S} \ll 1$, where $\eta ^{S}$ is the solvent viscosity and $\eta ^{P}$ is the polymeric viscosity. The Oldroyd-B limit in Li and Koch's work is obtained by equating the ‘mobility factor’ (denoted by $\alpha$ in their work) to zero (Bird *et al.* Reference Bird, Armstrong and Hassager1987). Then, it may shown that the polymeric stresses ($\boldsymbol {\tau }'_P$) are governed by $\boldsymbol {\tau }'_P + \lambda _1\overset {\nabla }{\boldsymbol {\tau '}}_P = 2\eta ^{P}\boldsymbol {D}'$, whereas the solvent stress ($\boldsymbol {\tau }'_S$) satisfies $\boldsymbol {\tau }'_S = 2\eta ^{S}\boldsymbol {D}'$, while the total stress is $\boldsymbol {\tau }' =\boldsymbol {\tau }'_P+\boldsymbol {\tau }'_S$. Noting that $\eta = \eta ^{S}+\eta ^{P}$ and non-dimensionalizing stresses with the characteristic scales of § 2.2, it may be shown that the total stress satisfies $\boldsymbol {\tau }+De\overset {\nabla }{\boldsymbol {\tau }} = 2[\boldsymbol {D}+De(1+\mathcal {C})^{-1}\overset {\nabla }{{\boldsymbol {D}}}]$. Comparing it with the constitutive relation of our study, as mentioned in (2.1*e*), we note that $\mathcal {C} = {\lambda _1}/{\lambda _2}-1$. Assuming $\mathcal {C} = {\lambda _1}/{\lambda _2}-1\ll 1$, Li and Koch expressed the electrophoretic velocity ($\mathcal {U}_{*}$) as an asymptotic expansion in $\mathcal {C}$ as: $\mathcal {U}_{*} = \mathcal {U}^{(0)}_{*}+\mathcal {C} \mathcal {U}^{(1)}_{*} + \ldots\,$, where $\mathcal {U}^{(0)}_{*} = 1$ and $\mathcal {U}^{(1)}_{*} = -1-\frac {2187}{2860}De^{2}_{*}$, enforcing the ‘mobility factor’ to be zero. Here, $De_*$ is the Deborah number defined in Li and Koch's work and it is related to our work as: $De_{*} = \beta \phi _P (1+\mathcal {C}) De$ and $\phi _P$ is the dimensionless potential on the particle surface, defined in § 4.3.1. For a uniformly charged particle with $a_0 = 1$ and $a_1 = 0$, we write, $\phi _P = \bar {\zeta }_0-\frac {1}{24}\bar {\zeta }_0^{3} + \ldots\,$. Also, in Li and Koch's work, the characteristic velocity was taken as $u_{c,LK} = \alpha (1+\mathcal {C})\phi _P$. Hence the dimensional form of $\mathcal {U}_{*}$ becomes

Further, inserting the expansion for $\phi _P$ as mentioned above with $De_* = \beta \phi _P De(1+\mathcal {C})$ and only retaining terms up to $O(\mathcal {C})$, (4.15) simplifies to

At the same time, in (4.12), enforcing $\lambda _1 = 1$, $\lambda _2 = (1+\mathcal {C})^{-1}$ and retaining terms only up to $O(\mathcal {C})$, it may be shown that the velocity reported in (4.16) is exactly recovered. We have therefore shown that our analysis is also able to successfully reproduce the reported results for uniformly charged particles, in viscoelastic fluids.

### 4.4. Comparison with numerical simulations of nonlinear viscoelastic models

As discussed in § 3.6, the Oldroyd-B constitutive model is known to produce quantitatively inaccurate results in the presence of strong elongational stresses. Here, such stresses are likely to be present within the EDL, wherein $\tau _{\theta \theta }$, $\tau _{\varphi \varphi }$ and $\tau _{\theta \varphi }$ all scale as $\delta ^{-2}$. Therefore, it is instructive to assess the accuracy of Oldroyd-B model itself, by comparing our analytical results with numerical solutions, carried out for more robust nonlinear viscoelastic models, which do not necessarily suffer from the same shortcomings. To this end, in this subsection we shall compare our results for the Oldroyd-B model with the numerical solutions obtained using the FENE-P constitutive model (Li & Koch Reference Li and Koch2020), which is one of the simplest nonlinear viscoelastic models that has reasonable mathematical tractability without sacrificing the essential physics of interest.

#### 4.4.1. The simplified governing equations for FENE-P model

We shall directly start with the dimensionless version of the governing equations, where all the variables are non-dimensionalized using the characteristic scales described in § 2.2. We reiterate that the electrostatic potential ($\phi$) is governed by the Poisson–Boltzmann equation (Kilic, Bazant & Ajdari Reference Kilic, Bazant and Ajdari2007) and the velocity field is governed by the Cauchy momentum equation along with the continuity equation for mass conservation. These equations are expressed as

*a*)\begin{gather} \nabla^{2}\phi = \delta^{{-}2}\sinh\phi, \end{gather}

*b*)\begin{gather}-\boldsymbol{\nabla} p +\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tau}+ \delta^{{-}2}\sinh(\phi)\boldsymbol{\nabla}\phi_{ext} = 0\quad \text{and}\quad \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}} = 0. \end{gather}

Here, $\phi _{ext}$ is the externally imposed potential and is given by (2.5). The total stress $\boldsymbol {\tau }$ for a FENE-P fluid may be written as $\boldsymbol {\tau } = \boldsymbol {\tau }^{S} + \boldsymbol {\tau }^{P}$, where $\boldsymbol {\tau }^{S}$ and $\boldsymbol {\tau }^{P}$ are the solvent and polymeric stresses respectively (Li & Koch Reference Li and Koch2020; Bird *et al.* Reference Bird, Armstrong and Hassager1987). Following § 4.3.2, we may define $\mathcal {C} = \eta ^{P}/\eta ^{S}$ as the ratio of polymeric and solvent viscosity and $\eta = \eta ^{S}+\eta ^{P}$ as the total viscosity, appearing in §§ 2.1 and 2.2. Then, $\boldsymbol {\tau }^{S} = 2(1+\mathcal {C})^{-1}\boldsymbol {D}$, while the polymeric stress satisfies the equation (Oliveira Reference Oliveira2002; Li & Koch Reference Li and Koch2020)

where $\ell$ is called the extensibility parameter, which is the ratio of the maximum allowed length of the springs in a bead-spring model of polymers to their equilibrium length (Oliveira Reference Oliveira2002). Therefore, the total stress ($\boldsymbol {\tau }$) is governed by the following equation:

which indicates that $\lambda _1 = 1$ and $\lambda _2 = (1+\mathcal {C})^{-1}$. It is well established that the FENE-P model can account for shear thinning of polymeric liquids, owing to the parameter $\ell$. Equation (4.19) reduces to the Oldroyd-B model of (2.1*e*), when $\ell \rightarrow \infty$. Since our analysis is valid for $O(\lambda _1)\sim O(\lambda _2)$, we shall compare the analytical solutions of § 4.2 with the numerical solutions of (4.17) in the limiting case of $\ell \gg 1$ and $\mathcal {C}\sim O(1)$.

The electrophoretic velocity of the particle ($\mathcal {U}{\hat {\boldsymbol {e}}}_u$) may be evaluated by holding the particle stationary and letting the surrounding fluid flow over it. Then, the far-field velocity will be equal to $-\mathcal {U}{\hat {\boldsymbol {e}}}_u$. Here, we shall consider an axisymmetric flow, where ${\hat {\boldsymbol {e}}}_u = {\hat {\boldsymbol {e}}}_z$. Therefore, in a frame fixed at the particle centre (see § 2.1), the following boundary conditions are satisfied:

*a*)\begin{gather} {\hat{\boldsymbol{e}}}_r\boldsymbol{\cdot}\boldsymbol{\nabla}\phi ={-}\delta^{{-}1}\zeta(\mu);\quad {\boldsymbol{v}} = 0,\quad \text{at},\ r=1, \end{gather}

*b*)\begin{gather}\boldsymbol{\tau}\rightarrow 0,\ p\rightarrow 0,\ \phi\rightarrow 0,\quad \text{as}\ r\rightarrow\infty. \end{gather}

The desired quantity may be obtained by noting that $\mathcal {U} = |{\boldsymbol {v}}(r\rightarrow \infty )|$.

#### 4.4.2. The numerical model

The numerical model is depicted in figures 2(*a*) and 2(*b*). Since we are looking into an axisymmetric flow, it is sufficient to consider any one half of the flow domain. The particle is held fixed with its centre at the origin. A cylindrical coordinate system ($z,\rho$) is used to compute the numerical solutions, where $r = \sqrt {z^{2}+\rho ^{2}}$. The particle carries a surface charge density of the form $\zeta (\theta ,\varphi ) = \bar {\zeta }_0 (a_0+a_1\mu )$. The relevant boundary conditions on all the surfaces are pictorially represented in the schematic 2(*a*). For the rest of this subsection, we have fixed the following parameters: $a_0 = a_1 = 0.5$, $De = \beta = 1$ and $\delta = 0.005$.

The numerical simulations are carried out in COMSOL Multiphysics 5.6, which uses finite element-based discretization. To ensure that the far-field free boundary is accurately represented, we have chosen $L_z = 7.5$ and $L_{\rho } = 10$. The flow domain is discretized using an unstructured triangular mesh. As our analysis suggests, the velocity and the stresses vary rapidly within the EDL and the polymeric stresses can be extremely large therein. Therefore, to accurately capture the variations within the EDL, an extremely fine mesh is taken within a distance $5\delta$ from the particle surface at $r = 1$. Outside the EDL, a relatively coarser mesh has been used. In particular, within the EDL, a ‘scaling factor’ of 30 is chosen while outside the same it is taken as 2. As a result, the minimum element size is $3\times 10^{-4}$ and the maximum element size (close to the boundaries) becomes $0.15$. In the region $1< r<1+5\delta$, there are 8562 domain elements while in aggregate there are 98 706 domain elements after discretization. Panel (*b*) represents a close-up view of the mesh around the particle surface. It may be clearly observed that the mesh size is extremely fine close to the particle and gradually becomes coarser away from it. The linear direct solver ‘MUMPS’ with relative tolerance $10^{-3}$ has been used to compute the numerical solutions. A representative example of the velocity contours obtained from the numerical simulations has been included in § S3 of the supplementary material.

#### 4.4.3. Comparison with analytical solutions

Figure 3(*a*,*b*) compares the electrophoretic velocities obtained from the analytical solutions using the Oldroyd-B model as reported in (4.11*a*) and the numerical solutions obtained using the FENE-P model. Panel (*a*) compares $\mathcal {U}$ as a function of $\bar {\zeta }_0$, whereas in panel(*b*) $\mathcal {U}$ vs $\ell ^{-2}$ is plotted; values of other relevant parameters are mentioned in the caption. As mentioned in § 4.4.1, the numerical simulations have been carried out in the limit of $\ell ^{2} \gg 1$ in order to recover the Oldroyd-B model.

From panel (*a*), we observe that the Oldroyd-B model can accurately predict the electrophoretic velocity until approximately $\bar {\zeta }_0\approx 0.15$ when $\ell ^{2} = 200$, given the particular choices of other relevant parameters. Recall that $\mathcal {C} = 0.5$ implies $\lambda _1 = 1$ and $\lambda _2 = 2/3$. This indicates that the Oldroyd-B model can reasonably describe the polymer stretching within the EDL only when two conditions are simultaneously met: (i) the surface charge density is small and (ii) the maximum allowable length of the polymers is large. When $\bar {\zeta }_0 > 0.15$, the Oldroyd-B model underpredicts the electrophoretic velocity. This is expected, because the FENE-P as well as many other nonlinear viscoelastic models predict shear thinning, which becomes especially important within the EDL, where the stretching rates are large. Therefore, as $\bar {\zeta }_0$ becomes larger, the shear thinning behaviour becomes increasingly important, which the Oldroyd-B model fails of capture. Noting that shear thinning essentially reduces the effective viscosity, the resulting electrophoretic velocity will naturally be larger as compared with that predicted by the Oldroyd-B model, as evident from panel (*a*). We have further observed that the comparison between the FENE-P and the Oldroyd-B model becomes more favourable when $\mathcal {C}$ is small (${\sim }0.1$). The reason may be attributed to the fact that, when $\mathcal {C} = 0$, Newtonian behaviour is recovered from (4.19) and hence, for small values of $\mathcal {C}$, both the Oldroyd-B and FENE-P models exhibit flow characteristics close to those of a Newtonian fluid, which results in a relatively better comparison between the two.

Panel (*b*), on the other hand, shows the comparison between the Oldroyd-B and FENE-P models for the electrophoretic velocity as a function of $\ell ^{-2}$, for different choices of $\bar {\zeta }_0$. Of course, $\mathcal {U}$ calculated using the Oldroyd-B model is independent of $\ell$ and hence these are represented by the horizontal straight lines in the figure. We observe that, for very small values of $\bar {\zeta }_0$ ($0.025$ and $0.05$), the Oldroyd-B and the FENE-P models agree closely for all values of $\ell \gg 1$. In other words, the polymer's maximum allowable length has very little influence on the translational velocity when the surface charge is weak. However, for larger values of $\bar {\zeta }_0$, especially when $\bar {\zeta }_0 \geqslant 0.15$, large differences are observed between the Oldroyd-B and the FENE-P models, which underlines the limitations of the former. The main reason may be attributed to the shear thinning behaviour embedded in the FENE-P model that eventually dominates when $\bar {\zeta }_0$ increases beyond a critical value (here $0.15$). In essence, one may therefore conclude that the Oldroyd-B model can predict the electrophoretic translational velocity in the thin EDL limit in a polymeric liquid with reasonable accuracy, when the surface charge is weak and the maximum permissible polymer lengths are very large. In other words, the effective Deborah number in the bulk, $De_{eff} = \bar {\zeta }_0 De$, needs to be small (close to $0.1$ or smaller) for favourable comparison between the Oldroyd-B and other nonlinear viscoelastic models. The condition related to the permissible length of the polymers stems from the fact that both the FENE-P and Oldroyd-B models are derived by treating the polymers as beads connected by springs (Van Heel, Hulsen & Van den Brule Reference Van Heel, Hulsen and Van den Brule1998). However, the Oldroyd-B model assumes the springs to be Hookean and infinitely extensible, while in FENE-P the springs are nonlinear and can only be stretched up to a maximum length (Van Heel *et al.* Reference Van Heel, Hulsen and Van den Brule1998). Therefore, it is easily realized that when this maximum allowable length is very large ($\ell \gg 1$), the FENE-P and the Oldroyd-B models should exhibit similar characteristics, as is indeed observed in figure 3. Finally, we note from figure 3 that the accuracy of the predictions using the Oldroyd-B model is far more sensitive to $\bar {\zeta }_0$ than $\ell$.

The differences between the Oldroyd-B and the FENE-P models, particularly inside the EDL, can be better understood by looking into a simple example of unidirectional electroosmotic flow over a flat plate, carrying uniform surface/zeta potential $\bar {\zeta }_0$. For analytical simplicity, we shall consider the special case of $\mathcal {C}\gg 1$, i.e. very large polymeric viscosity as compared with the solvent viscosity. For this particular example, the $x$-axis runs along the plate and the $y$-axis runs vertical to the plate. Choosing the characteristic length scale as $H$ (e.g. the plate length), taking $\delta = \lambda _D/H$ and keeping all other characteristic scales and the non-dimensionalization scheme intact as in § 2.2, the FENE-P constitutive model simplifies to the following form for the flow under consideration (Oliveira Reference Oliveira2002):