Electrophoretic motion of a non-uniformly charged particle in a viscoelastic medium in thin electrical double layer limit

Abstract The electrophoretic motion of a non-uniformly charged particle in an Oldroyd-B fluid is analysed here in the limit of thin electrical double layers. To this end, we analytically derive expressions for the modified Smoluchowski slip velocity around the particle, carrying weak but otherwise arbitrary surface charge. Our analysis reveals that the modified Smoluchowski slip around a particle differs significantly in a viscoelastic medium as compared with Newtonian fluids. The flow field thus derived is applied to two special cases of non-uniformly charged particles to obtain a closed-form expression for their electrophoretic translational and rotational velocities. We show that the particle's velocity strongly depends on its size in a viscoelastic medium, even for weakly charged surfaces, which is in stark contrast to the well-established theory for Newtonian fluids for weakly charged particles with negligible surface conduction. We further demonstrate that the presence of non-uniform surface charge enhances the influence of the medium's viscoelasticity on the particle's translational as well as angular velocity and this effect strongly depends on the nature of surface charge distribution. Such a physical paradigm, which leads to a breaking of fore–aft symmetry that is unique to complex fluids despite operating in the regime of creeping flows. Our study provides new theoretical framework for understanding electrophoresis of charged entities (such as DNA or active matter) in complex fluids, including biologically relevant fluidic media.

It is well established that fluid media of emerging interest in electrically mediated transport of biological entities commonly exhibit viscoelastic behaviour (Skalak, Ozkaya & Skalak 1989;Brust et al. 2013). However, a majority of studies on electrophoresis (Saville 1977;Baygents & Saville 1991;Ohshima 2006;Khair & Squires 2009) tend to focus on particle motion in Newtonian fluids, while only a handful of investigations (Hsu, Yeh & Ku 2006;Hsu & Yeh 2007;Khair et al. 2012;Posluszny 2014;Li & Koch 2020) have addressed the phenomenon in a non-Newtonian medium. Among those, many of the investigations concern Carreau-type constitutive models (Bird, Armstrong & Hassager 1987;Khair et al. 2012), 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. 1987;Ghosh et al. 2016). On the other hand, viscoelastic constitutive models, such as the retarded motion relations (Chan & Leal 1979) as well as the differential constitutive relations (Bird et al. 1987;Mukherjee & Sarkar 2010;, hold certain favourable characteristics in capturing many of the essential rheological features of various polymeric (Bird et al. 1987; Barron et al. 1995;Li & Koch 2020) as well as biological fluids (Yeleswarapu et al. 1998;Brust et al. 2013). As a result, numerous types of viscoelastic flows have been widely studied (for instance, Vamerzani, Norouzi & Firoozabadi 2014;Turkoz et al. 2018) over the years, although rigorous studies on electrophoretic motion through viscoelastic medium are extremely scarce. Lu et al. (2014Lu et al. ( , 2015 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 (2020) 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 1977;Ohshima 2006;Khair & Squires 2009;Li & Koch 2020). However, in many cases, the moving particles might have non-uniformly distributed charges on their surfaces, either naturally induced or externally imposed. As Anderson (1985) 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. 2000;Velegol 2002), various inorganic particles (Van Riemsdijk et al. 1986;Velegol 2002) and particles with adsorbed surfactants (Fleming, Wanless & Biggs 1999) might naturally exhibit non-uniform surface charge density. At the same time, 'Janus particles' (Molotilin, Lobaskin & Vinogradova 2016) 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 1985;Yoon 1991;Velegol 2002) 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 1989;Solomentsev & Anderson 1994;Yariv 2005). 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. 1987).
We confine our attention to the thin EDL limit (thin EDL) (Ajdari 1995;Mandal et al. 2015) 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 σ (θ, ϕ), suspended in a viscoelastic medium of viscosity η, permittivity ε, relaxation time (Bird et al. 1987) λ 1 and a retardation time λ 2 . The fluid also contains dissolved electrolytes, which dissociate into ions and form an EDL around the particle surface (Ohshima 2006;Poddar et al. 2016), 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 , θ, ϕ), 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 U ê u , whereê u is the unit vector pointing towards the direction of the particle's motion. In general,ê u is not constrained to be oriented along the z-axis (Anderson 1985). In addition, because of non-uniform surface charge, the particle may also undergo rotational motion (Anderson 1985;Yoon 1991;Solomentsev & Anderson 1994) with angular velocity Ω , much like other anisotropic particles.

Electrophoresis in viscoelastic medium
x Radius: a Figure 1. Schematic of a spherical particle of radius a, carrying an arbitrary surface charge density σ (θ, ϕ). The particle is translating with velocity U ê u in a viscoelastic medium, subject to an externally applied electric field (magnitude far from the particle E 0 ) along the z direction. The fluid has viscosity η, relaxation time λ 1 , retardation time λ 2 and permittivity ε. Schnitzer et al. 2013), 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 (λ 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 ∼1-100 nm (Ajdari 1995;Bandopadhyay & Chakraborty 2012a). 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 ∼ O(1) (Saville 1977). 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 1979). 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. 1987). Fifth, we shall disregard the presence of a depletion layer (Pranay, Henríquez-Rivera & Graham 2012;Mukherjee et al. 2017) 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 2011). To compare physically realistic values of R g with those of λ D (the characteristic EDL thickness), we consider the example of poly-ethylene glycol in water (Cruje & Chithrani 2014). With a monomer length of ≈0.35 nm for n = 100 segments, one arrives at R g ≈ 1.43 nm; for n = 1000 segments, R g ≈ 4.5 nm. Therefore, for λ D ∼ O(10 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 2009;Li & Koch 2020). 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. Ω / = 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.

The characteristic scales
Electrophoretic motion entails several important characteristic scales for a number of different variables, dictating the flow dynamics. For any variable ξ , we denote it's characteristic scale by ξ c . As such, the following characteristic scales are chosen: the characteristic potential: ψ 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 1977), u c = εk 2 T 2 /e 2 ηa; the characteristic concentration: c c = c 0 . Based on these choices, the characteristic stress reads: τ c = ηu c /a, whereas the characteristic surface charge is σ c . Finally, the characteristic relaxation time may be equated to λ 0 = max(λ 1 , λ 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, λ D = εkT/2c 0 e 2 , with δ = λ D /a -as the non-dimensional EDL thickness; (ii) the characteristic surface potential of the particle (Ajdari 1995), ζ c = σ c λ D /ε, withζ 0 = eζ c /kT as the characteristic potential relative to the thermal potential; (iii) the relative strength of the external field may be expressed through (Saville 1977) β = eE 0 a/kT; (iv) the ionic Péclet number may be defined as (Saville 1977) 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 λ 0 /a. While various alternative approaches to quantify the extent of viscoelasticity have been reported (Li & Koch 2020); 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,ζ 0 1; 'thin EDL' implicates δ 1; the strengths of advection and applied field are constrained by: β ∼ O(1) and Pe ∼ O(1). It needs to be clarified here that, despite De ∼ O(1), the limit of weak surface charge (ζ 0 1) implies that the dimensionless electrokinetic velocity would scale as O(ζ 0 ). Therefore, the confluence of electromechanics and hydrodynamics would lead to an effective Deborah number given by De eff ∼ O(ζ 0 De) 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 ∼O(1). As a consequence, it is not necessary to further assume low polymer concentration (Li & Koch 2020), which would require λ 1 /λ 2 − 1 1 -see § 3.5 for further discussion of this.

The governing equations
The transport processes are governed by the Poisson-Nernst-Planck-Cauchy momentum equations (Saville 1977;Ghosh et al. 2016), 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 ξ is expressed as ξ = ξ /ξ c . The Nernst-Planck equations are reformulated in terms of the net charge density and concentration, defined as (Schnitzer & Yariv 2014) 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 1977;Bird et al. 1987) (2.1e) In the above equations, τ is the stress field, p is the pressure, ψ is the electrical potential, D is the rate of strain tensor (= 1 2 [∇v + (∇v) T ]) and v is the velocity field. Further, ∇ τ and ∇ D indicate the first convected derivatives of the stress and strain rate. For any second-rank tensor A, it's convected derivative ∇ A is expressed as (Bird et al. 1987) The above equations are subject to the following boundary conditions: Note that, in (2.2c), the boundary condition for velocity at the particle surface accounts for its rotation. In (2.2b), ζ(θ, ϕ) is analogous to the surface potential and ζ ∼ O(ζ 0 ) is it's characteristic magnitude. For convenience, we define, ζ(θ, ϕ) =ζ 0ζ (θ, ϕ), such thatζ ∼ O(1). For electrophoretic motion, Uê u and Ω, 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 1985;Ohshima 2006;Goswami et al. 2017): It is important to note here that the above 'force-free' conditions only remain valid in an unbounded medium, as previously noted by Yariv (2006). 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 δ → 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 δ → 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.3a) and (2.3b) vanish (Ye et al. 2002;Chen & Keh 2014) 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 particle + the EDL) becomes zero (Ye et al. 2002;Chen & Keh 2014). Equations (2.3a) and (2.3b) may then be re-written as S p τ ·ê r dS = 0 and S pê r × (τ ·ê r ) dS = 0. (2.4a,b) Notice that the integration is to be performed over the surface S p , which has a radius r ∼ 1 + δ and includes the EDL as well. Finally, we note that the potential may be conveniently split into two parts as follows (Bahga, Vinogradova & Bazant 2010;Ghosh, Mandal & Chakraborty 2017): ψ = φ + φ ext , where φ ext is the contribution from only the externally imposed electric field and φ is the contribution of the particle's surface charge to the total electrostatic potential. It may then be deduced that φ ext satisfies the equation ) ∇ 2 φ ext = 0, subject tô e r · [∇φ ext ] r=1 = 0 and [∇φ ext ] r→∞ = −βê z . Hence, φ ext has the solution φ ext = −β r + 1 2r 2 P 1 (μ), (2.5) where, μ = cos θ and P n (x) is the Legendre polynomial of the first kind of order n. As a consequence, the component φ would satisfy the following equation (derived from (2.1c)): subject to the following conditions: The Nernst-Planck equations (2.1a) and (2.1b) and the corresponding boundary conditions (2.2a) 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 (ζ ) to change dynamically. The relation between the particle's angular velocity and evolution of ζ has been included in § S1.1 of the supplementary material.

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 1995;Squires & Bazant 2004;Ghosh et al. 2016). 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.1c), we note that the thin EDL limit (δ 1) is a singular problem, as also pointed out by a number of earlier studies (Yariv 2009;Schnitzer & Yariv 2012;Ghosh et al. 2016). This calls for the application of a matched asymptotic expansion (Leal 2007;Bender & Orszag 2013), 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 δ and lies next to the particle surface and (b) the bulk, or 'outer layer', with characteristic length scale r ∼ O(1). In both the layers, all variables may be expanded in an asymptotic series of δ as follows (Leal 2007;Yariv 2009): (3.1) Note that this expansion also applies to U and Ω. 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 δ 1 and De ∼ 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 (ζ 0 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.

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.1a)-(2.1e), which take the following form: In (3.2d), T = ∇ τ and S = ∇ D. Detailed expressions for T and S in spherical coordinates may be found in Bird et al. (1987). Note that, φ ext has already been determined and thus there is no need to expand it in δ. 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 δ, 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 2012), R = (r − 1)/δ, which is O(1) inside the inner layer, since r − 1 ∼ δ. Accordingly, the velocity components would rescale as follows: (1) and from the continuity equation, u r → δV ∼ O(δ), with V ∼ O(1). In the above, u r , u θ and u ϕ are the radial, polar and azimuthal components of velocity, respectively. There are no changes in the scaling of the potential, charge and salt concentration: The rescaling of the stresses and the strain rates are slightly more involved (Saprykin, Koopmans & Kalliadasis 2007;Ghosh et al. 2016) and need to be worked out based on the constitutive model (2.1e). 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 τ rθ ) would scale as O(δ −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 τ θθ ) scale as O(δ −2 ) in the inner layer. In addition, similar scaling is also observed for the component τ θϕ . 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. 2007;Ghosh et al. 2016), 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 S and D do not exhibit the same scaling, contrary to the stress (τ ) and its convected derivative (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: Recall that μ = cos θ. Also, from (2.5), in the inner layer, E ext r = 0 and E ext θ = − 3 2 β 1 − μ 2 to the leading order in δ, where E ext = −∇φ 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.4d)-(3.4f ) 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 (δ 1). It is well documented (Bird et al. 1987) that, in perfectly unidirectional shear flows (such as Couette or Poiseuille flows), the shear stress (τ xy , x being the direction of flow) varies as τ xy = η * (γ )γ , whereγ is the rate of strain and η * (γ ) is the appropriately non-dimensionalized shear dependent viscosity. On the other hand, the first normal stress difference varies as τ xx − τ yy = η N,1 (γ )γ 2 , where η N,1 is the (unitless) first normal stress coefficient. Because we have considered the Oldroyd-B constitutive model, it follows that both η * and η N,1 are constants (more discussion on this is provided in § 3.6). Since the shear rate in the EDL isγ ∼ δ −1 , it immediately follows that in the EDL on a perfectly plane surface τ xy ∼ δ −1 and τ xx ∼ δ −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 θ and ϕ, using spherical coordinates. Thus the stresses in the streamwise directions are τ θθ , τ ϕϕ and τ θϕ , while stresses equivalent to τ xy are τ rθ and τ rϕ . As a result, we expect that τ rθ , τ rϕ ∼ O(δ −1 ) and τ θθ , τ ϕϕ , τ θϕ ∼ O(δ −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. 2007;Ghosh et al. 2016). On a perfectly flat surface with uniform surface charge, the streamwise stresses (τ 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 τ θθ , τ ϕϕ 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(δ −2 ). At the same time, because the cross-stream gradients of τ rθ and τ rϕ govern the velocity field, these gradients also scale as δ −2 , despite τ rθ and τ rϕ scaling as δ −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 δ. 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 δ. The only way curvature affects the flow is through the θ 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 δ, 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:τ ij + λ 1 DeT ij = 2[D ij + λ 2 DeS ij ], when ij ≡ rr, rθ and rϕ, (3.5a) τ ij + λ 1 DeT ij = 2λ 2 DeS ij , when ij ≡ θθ, θϕ and ϕϕ. (3.5b) The detailed expressions for the rescaled convected derivatives (T andS) 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: are respectively the θ and ϕ components of velocity at the particle surface due to it's rotation.

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. 2016 In addition, the net salt and charge fluxes across the EDL-bulk interface also need to be matched (Yariv 2009;Ghosh et al. 2016Ghosh et al. , 2017 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 2009;Ghosh et al. 2016): Finally, it is important to note that the quantities, lim R→∞ U and lim R→∞ W may be combined to write v S = lim R→∞ [Uê θ + Wê φ ] − Ω ×ê r . We identify the quantity 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.2c) may be solved subject to (3.3) in the far field and v = Ω ×ê r + v S and u r = 0 at r = 1.

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 (ζ 0 1). Accordingly, all the variables (both in the inner and outer layers) may be further expanded ) in a regular asymptotic series ofζ 0 as Here, ξ may represent any variable such as U, V, v, φ, . . . and so on. Recall that the above expansion is applied to the variables which already denote leading-order terms in δ. We re-emphasize that, in the regular expansion,ζ 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.

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

Solutions to the inner layer equations
The leading order (in δ) inner layer equations may also be solved using the asymptotic expansion mentioned in (3.9), subject to conditions (3.6a)-(3.6c). To this end, we note that the inner layer variables exhibit the following expansions inζ 0 : Implications and the physical basis of the asymptotic expansions mentioned above deem further elaboration. Notice that, at the leading order ofζ 0 , the nonlinear polymeric contributions to the stresses are absent, which indicates that at O(ζ 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(ζ 0 ). As a result, the rescaled Newtonian part of the stresses (also refer to table 1) would scale as O(ζ 0 ) and the convected derivatives (T andS), which arise from the polymeric contributions to the stresses, would scale as O(ζ 2 0 ), 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.1e), obtained by substituting λ 1 = λ 2 = 0. Therefore, at the leading order ofζ 0 , i.e. at O(ζ 0 ), the flow is asymptotically Newtonian and the nonlinearities arising from the polymeric stresses only contribute from O(ζ 2 0 ) 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 ∂τ rθ /∂r) scale as O(ζ 0 δ −2 ) for weak surface charge, since inside the EDL, r ∼ 1 + δ and τ rθ , etc. ∼ O(δ −1 ). On the other hand, the normal stress derivatives (terms like ∂/∂μ(τ θθ 1 − μ 2 )) scale as O(ζ 2 0 δ −2 ) under the same conditions, since the normal stresses themselves are O(δ −2 ) (see table 1), whereasT andS ∼ O(ζ 2 0 ) -see (3.5). Sinceζ 0 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. ( . 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.4a)-(3.4b), subject to the boundary conditions, (3.6a) and (3.6b). As already mentioned, at this order, the velocity field is identical to that of a Newtonian fluid, when nominally De ∼ O(1). The velocity field thus has the solution is the Gegenbauer polynomial of the first kind and order n (Leal 2007), expressed as Q n (x) = x −1 P n (u) du. For example (Leal 2007 y sin(θ )) sin(ϕ). The leading-order slip velocity without rotation is thus given by v (1) (3.14) (ii) The O(ζ 2 0 ) velocity field: at O(ζ 2 0 ), 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(ζ 2 0 ), Φ 2 = Π 2 = 0, C 2 = Φ 2 1 and from the r-momentum equation, P 2 = 1 2ζ 2 e −2R . Combining (3.4c)-(3.4f ), along with (3.5), it may be shown that, U 2 and W 2 are governed by the following equations: In the above, ω 3 = 1 − μ 2 χ 1,ϕ + μχ 1 / 1 − μ 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 More discussion on the nature of the velocity profiles at O(ζ 2 0 ) have been included in § 3.5.
The governing equations for the velocity components may be derived by inserting the O(ζ 3 0 ) stresses into the governing equations for the inner layer. It may be verified that the O(ζ 3 0 ) velocities satisfy equations of the following form: In the above, components of can be deduced from (3.5) in combination with (A1) and (A2) in Appendix A. Detailed expressions for various components of 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 = Γ 3 and W 3 = χ 3 at R = 0 and both remain bounded as R → ∞. 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(ζ 3 0 ) contribution to the slip velocity without particle rotation, takes the following form: The complete modified Smoluchowski slip velocity at O(ζ 3 0 ) 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 The contributions at the individual orders ofζ 0 are given in (3.14), (3.18) and (3.20).

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ζ 0 , instead of De, which in many cases is the usual choice (Bird et al. 1987). In this regard, it should be noted here that, although nominally De ∼ O(1), because the velocity is O(ζ 0 ) (on account of weak surface charge), the effective Deborah number becomes De eff =ζ 0 De 1. Therefore, the regular expansion inζ 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(ζ 2 0 ) 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 (C 1), which translates into λ 1 /λ 2 − 1 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 2020).
The basic physics behind the O(ζ 2 0 ) and the O(ζ 3 0 ) equations may be appreciated as follows. The leading-order flow (at O(ζ 0 )) is effectively Newtonian, which stretches the polymers in the EDL, thus giving rise to excess polymeric stresses at O(ζ 2 0 ). These excess polymeric stresses are non-uniform along the particle surface, because of its curvature and variations inζ(θ, ϕ). As a result, the gradients of these excess stresses drive a flow at O(ζ 2 0 ), which is observed in (3.15). Because of the nonlinear constitutive relation of the fluid itself, this O(ζ 2 0 ) 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(ζ 3 0 ) 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 1, which indicates that, for U ∼ 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 τ θθ , τ ϕϕ and τ θϕ all scale as δ −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θ component (as an example) may be rewritten as whereT rθ andS rθ are given in Appendix A and they denote the effects of polymeric stresses inside the EDL. Since De EDL ∼ δ −1 1, δDe EDL ∼ O(1) and hence these terms cannot be neglected, although they are multiplied with δ. 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 1). As a result, it may be inferred that the analysis in § 3.4 is indeed valid when De ∼ O(1), which would imply De eff =ζ 0 De 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 δDe EDL ∼ 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 ω 1 =ζ(θ, ϕ)Q 1 (μ) and its derivatives, which encodes information about variations in surface charge density as well as the particle's curvature. While ω 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 ω 1 containζ(θ, ϕ) and its derivatives, which underscore the effects of variations in surface charge density, while the factor Q 1 (μ)/ 1 − μ 2 bears the signature of the particle's curvature. Notice that, even forζ = a constant, i.e. for a uniformly charged particle, the derivatives of ω 1 are in general non-zero because of Q 1 (μ) appearing in them, which indicates that the curvature will affect the slip velocity, appearing in the O(ζ 2 0 ) and O(ζ 3 0 ) 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 1995;Yariv 2009). 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 β, as evident from (3.18) and (3.20). The third important point to note is that, for λ 1 = λ 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. 1987). 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(ζ 3 0 ) contributions have been reported. Notice from (3.18) that, in presence of particle rotation (when ω 3 / = 0), the slip velocity at O(ζ 2 0 ) also has a component along theê ϕ direction, resulting in anisotropic motion. When the medium is Newtonian, i.e. De = 0, this component vanishes. We therefore conclude that the slip velocity (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(ζ 2 0 ) onwards. Although the expression of the Smoluchowski slip velocity derived above applies to an arbitrary distribution ofζ , 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γ 2 (γ 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. 1987). 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. (1987).
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. 1987;Afonso et al. 2009) 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. 2009). 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 1983;Bird et al. 1987;Tan & Masuoka 2005;Aggarwal & Sarkar 2008;Mukherjee & Sarkar 2011) 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ζ(μ, ϕ) = ∞ n=0 a n P n (μ). Here, to keep the algebra tractable, we choose a 0 / = 0 and a 1 / = 0, while a n = 0, ∀n 2; hence, the charge density is given byζ(μ, ϕ) = a 0 P 0 (μ) + a 1 P 1 (μ) = a 0 + a 1 μ. It may be checked that, in this case, the particle carries a net charge (appropriately non-dimensionalized) amounting to 4π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. Ω = 0 and (ii) the particle will translate along the z-axis, i.e.ê u =ê 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: Insertingζ = a 0 + a 1 μ into (4.1a,b), we deduce that the slip velocity for this particular case looks as follows: In (4.2c), 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. ζ = 1). It may be easily verified that, in such a case, v The components v (1) S,θ and v S,θ may also be evaluated by inserting a 0 = 1 and a 1 = 0 into (4.2a) and (4.2b), 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 (μ = cos θ ) vanish, the only alterations to the slip velocity comes from the curvature of the particle, as denoted by the factor Q 1 (μ)/ 1 − μ 2 .
The electrophoretic velocity U may be ascertained by solving the governing equations in the outer layer, i.e. (3.10a) and (3.10b), subject to the far field condition (3.11a) and the matching condition (3.11b) at the particle surface, wherein v S (=v S,θêθ ) 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.4a,b). Here, the force balance effectively reduces to the following form (Leal 2007) Following the analysis for the EDL, the outer layer variables are also expanded in a regular asymptotic series ofζ 0 . Therefore, the electrophoretic velocity and the net force on the particle are expanded as where, All other variables may also be expanded in a similar way. The kth-order (inζ 0 ) outer layer equations may be expressed as (for k = 1, 2, 3, . . . etc.) −∇p k + ∇ · τ (k) = 0 and ∇ · v k = 0, (4.5a) (4.5b) These are subject to the following boundary conditions: The force balance at every order ofζ 0 reads F (k) z = 0, which may be used to derive U k . Note that, using (4.5b), the momentum equation (4.5a) may be rewritten as follows: −∇p k + ∇ · [∇v k + (∇v k ) T ] + De∇ · (2λ 2 S (k) − λ 1 T (k) ) = 0 and ∇ · v k = 0. For an axisymmetric flow, (4.5a) may be solved using a streamfunction (Leal 2007), defined as (for kth order ofζ 0 ) u (k) r = −(1/r 2 )(∂Ψ k /∂μ) and u (k) θ = −(1/r 1 − μ 2 )(∂Ψ k /∂r), where Ψ is the streamfunction. As such, the momentum equation at any order ofζ 0 may be expressed in terms of Ψ as follows (Leal 2007;Goswami et al. 2017): where (Leal 2007),
(4.12) The complete solution for the O(ζ 3 0 ) streamfunction for a uniformly charged particle has been included in the supplementary material (see § S2 therein). Further discussion is included in the next subsection.

Comparison with previous studies
We shall first compare our results for the electrophoretic velocity (U ) with two of the previous studies: (i) non-uniformly charged particle in a Newtonian medium, carried out by Anderson and coworkers (Anderson 1985;Fair & Anderson 1989) and (ii) uniformly charged particle in a viscoelastic medium with weak polymeric viscosity, investigated by Li & Koch (2020). 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.

Comparison with results for Newtonian fluids
Anderson and coworkers (Anderson 1985;Fair & Anderson 1989) 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 φ P (non-dim. φ P = φ P /ψ c ), while the electrophoretic velocity (with units) in Newtonian fluids is denoted by U N . Anderson (1985) showed that the electrophoretic velocity for an arbitrary distribution of φ P on the particle surface is given by where E ∞ is the far-field imposed electric field (uniform), H 2 = φ P Y 2 is the quadrupole moment of φ P and Y 2 = 3nn − I is the second spherical harmonic. Also, · denotes the average over the particle surface. Referring to the description in § 4.1, we note that, here, E ∞ = E 0êz , U N = U Nê z ,n =ê r and φ P is related to the dimensionless surface charge density (σ =ζ 0ζ ) as φ P = 2ψ c sinh −1 ( 1 2ζ 0ζ (μ)). As a result, forζ 0 1, φ P has the expansion, φ P /ψ c =ζ 0ζ − 1 24ζ 3 0ζ 3 + . . . , whereζ = a 0 + a 1 μ. Equation (4.13) may then be simplified to U N = 3εE 0 /2η[ φ P I − ê rêr φ P ] :ê zêz and hence U N may be expanded as 3 ] :ê zêz . Takingζ = a 0 + a 1 μ, we obtain ζ = a 0 and ê rêr ζ = 1 3 a 0 I and hence, U (1) N = αa 0 . At the same time, ζ 3 = a 0 (a 2 0 + a 2 1 ) and ê rêrζ 3 = a 0 (a 2 0 /3 + a 2 1 /5)I + (4.14) It may be noted that, by inserting λ 0 = 0 (Newtonian limit) in (4.11a), 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.

Comparison with results for viscoelastic fluids
Li & Koch (2020) have analysed the electrophoretic mobility of a uniformly charged sphere in a Giesekus fluid. They considered the limit of small polymeric viscosity (C 1, defined later), weak surface charge (ζ 0 1) and thin EDL. Note that we have used different symbols (to those used by Li & Koch) to express their results. Li & Koch (2020) derive their results for C = η P /η S 1, where η S is the solvent viscosity and η P is the polymeric viscosity. The Oldroyd-B limit in Li and Koch's work is obtained by equating the 'mobility factor' (denoted by α in their work) to zero (Bird et al. 1987). Then, it may shown that the polymeric stresses (τ P ) are governed by τ P + λ 1 ∇ τ P = 2η P D , whereas the solvent stress (τ S ) satisfies τ S = 2η S D , while the total stress is τ = τ P + τ S . Noting that η = η S + η P and non-dimensionalizing stresses with the characteristic scales of § 2.2, it may be shown that the total stress satisfies τ + De Comparing it with the constitutive relation of our study, as mentioned in (2.1e), we note that C = λ 1 /λ 2 − 1. Assuming C = λ 1 /λ 2 − 1 1, Li and Koch expressed the electrophoretic velocity (U * ) as an asymptotic expansion in C as: 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 * = βφ P (1 + C)De and φ 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, φ P =ζ 0 − 1 24ζ 3 0 + . . . . Also, in Li and Koch's work, the characteristic velocity was taken as u c,LK = α(1 + C)φ P . Hence the dimensional form of U * becomes U * = αφ P − αφ P 2187 2860 De 2 * C + . . . . At the same time, in (4.12), enforcing λ 1 = 1, λ 2 = (1 + C) −1 and retaining terms only up to O(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.

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 τ θθ , τ ϕϕ and τ θϕ all scale as δ −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 2020), which is one of the simplest nonlinear viscoelastic models that has reasonable mathematical tractability without sacrificing the essential physics of interest.

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 (φ) is governed by the Poisson-Boltzmann equation (Kilic, Bazant & Ajdari 2007) and the velocity field is governed by the Cauchy momentum equation along with the continuity equation for mass conservation. These equations are expressed as (4.17a) −∇p + ∇ · τ + δ −2 sinh(φ)∇φ ext = 0 and ∇ · v = 0. (4.17b) Here, φ ext is the externally imposed potential and is given by (2.5). The total stress τ for a FENE-P fluid may be written as τ = τ S + τ P , where τ S and τ P are the solvent and polymeric stresses respectively (Li & Koch 2020;Bird et al. 1987). Following § 4.3.2, we may define C = η P /η S as the ratio of polymeric and solvent viscosity and η = η S + η P as the total viscosity, appearing in § § 2.1 and 2.2. Then, τ S = 2(1 + C) −1 D, while the polymeric stress satisfies the equation (Oliveira 2002;Li & Koch 2020) (4.18) where 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 2002). Therefore, the total stress (τ ) is governed by the following equation: (4.19) which indicates that λ 1 = 1 and λ 2 = (1 + C) −1 . It is well established that the FENE-P model can account for shear thinning of polymeric liquids, owing to the parameter . Equation (4.19) reduces to the Oldroyd-B model of (2.1e), when → ∞. Since our analysis is valid for O(λ 1 ) ∼ O(λ 2 ), we shall compare the analytical solutions of § 4.2 with the numerical solutions of (4.17) in the limiting case of 1 and C ∼ O(1). The electrophoretic velocity of the particle (Uê 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 −Uê u . Here, we shall consider an axisymmetric flow, whereê u =ê z . Therefore, in a frame fixed at the particle centre (see § 2.1), the following boundary conditions are satisfied: e r · ∇φ = −δ −1 ζ(μ); v = 0, at, r = 1, (4.20a) The desired quantity may be obtained by noting that U = |v(r → ∞)|.

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, ρ) is used to compute the numerical solutions, where r = z 2 + ρ 2 . The particle carries a surface charge density of the form ζ(θ, ϕ) =ζ 0 (a 0 + a 1 μ). 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 = β = 1 and δ = 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 ρ = 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δ 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 × 10 −4 and the maximum element size (close to the boundaries) becomes 0.15.  Figure 2. (a) Schematic of the numerical model is depicted for an axisymmetric flow, with a cylindrical coordinate system (z, ρ) fixed at the particle centre. The particle carries a surface charge density of the form ζ(μ) =ζ 0 (a 0 + a 1 μ). The domain size is 2L z × L ρ , where L z = 7.5 and L ρ = 10. All the stresses and the potential vanish on the boundaries and along ρ = 0, all the derivatives vanish. (b) Close-up of the mesh around the particle is shown. The domain is discretized using a triangular unstructured mesh. A large number of grid points are taken within a distance 1 + 5δ from the particle, to accurately capture the variations within the EDL.
In the region 1 < r < 1 + 5δ, 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. Figure 3(a,b) compares the electrophoretic velocities obtained from the analytical solutions using the Oldroyd-B model as reported in (4.11a) and the numerical solutions obtained using the FENE-P model. Panel (a) compares U as a function ofζ 0 , whereas in panel(b) U vs −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 2 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ζ 0 ≈ 0.15 when 2 = 200, given the particular choices of other relevant parameters. Recall that C = 0.5 implies λ 1 = 1 and λ 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ζ 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ζ 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 C is small (∼0.1). The reason may be attributed to the fact that, when C = 0, Newtonian behaviour is recovered from (4.19) and hence, for small values of 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.

Comparison with analytical solutions
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 −2 , for different choices ofζ 0 . Of course, U calculated using the Oldroyd-B model is independent of and hence these are represented by the horizontal straight lines in the figure. We observe that, for very small values ofζ 0 (0.025 and 0.05), the Oldroyd-B and the FENE-P models agree closely for all values of 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ζ 0 , especially whenζ 0 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ζ 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 =ζ 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 1998). 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. 1998). Therefore, it is easily realized that when this maximum allowable length is very large ( 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ζ 0 than .
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ζ 0 . For analytical simplicity, we shall consider the special case of C 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 δ = λ 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 2002): Invoking the Debye-Huckel linearization forζ 0 1, the potential distribution takes the form φ =ζ 0 e −y/δ . Accordingly, it may be shown that the velocity field takes the form (y = 0 on the flat plate) De 2 β 3 δ −2ζ 3 0 1 − e −3y/δ + . . . .

(4.22)
Our analytical solutions for the Oldroyd-B model assume that viscoelasticity inside the EDL remains subdominant so that the velocity field scales similarly to that of a Newtonian fluid. However, we observe from (4.22) that the viscoelastic contribution to the velocity scales as δ −2 when shear thinning effects are present and hence this component may easily dominate in the thin EDL limit (δ 1), whenζ 0 is larger than a critical value. As a result, the analysis of flow within the EDL, as outlined in § 3, remains valid for the FENE-P model as well so long as the viscoelastic contributions are negligible at the leading order of ζ 0 . From (4.22), we note that such a paradigm is realized whenζ 2 0 δ 2 6 ( 2 − 3) −2 . By inserting the values of δ and used in figure 3(a), we observe that this condition is satisfied whenζ 0 0.07, which is in line with the conclusions drawn at the end of the previous paragraph. We would like to clarify that the above constraint on the value ofζ 0 only applies when the FENE-P or another similar nonlinear viscoelastic model is being considered. As far as the Oldroyd-B model is concerned,ζ 0 1 is sufficient for the validity of the analysis in § 3.4. Further, we note from (4.21) that the scaling τ xx ∼ O(τ 2 xy ) equally holds for the FENE-P model as well, despite its distinctive velocity field, which underlines the ability of the Oldroyd-B model in capturing many of the critical flow features within the EDL.
Second, it is interesting to note the unique curvature dependence of the electrophoretic velocity of the particle, as evident from (4.11). Recall that, for a Newtonian fluid, the electrophoretic velocity is independent of its size, as evident from the expression for U N mentioned above (as also reported previously in the literature, see Ohshima 2006). However, notice from (4.11a) that the particle's radius a explicitly appears in the expression for the velocity and is always associated with the characteristic relaxation/retardation time scale λ 0 -a viscoelastic property of the fluid. Further, observe that particles of larger sizes tend to cause smaller deviations from the Newtonian behaviour -this non-intuitive coupling between the particle's curvature and the fluid rheology is unique to the complex fluid under consideration here. While this coupling is present for a uniformly charged particle, it is augmented by non-uniform surface charge, as evident from the λ 0 /a appearing in the O(ζ 2 0 ) term in the particle velocity. Another point to note in this regard is that the particle size and the non-uniformity in charge density are intricately coupled in controlling the net change in the velocity (whether positive or negative) as compared with the case of a Newtonian fluid medium -more discussion on this is included later. Third, it is to be noted that the electrophoretic velocity is not proportional to the applied external field β, as is expected for Newtonian fluids, in the presence of thin EDLs. The O(ζ 2 0 ) and O(ζ 3 0 ) corrections introduce nonlinear dependence on the external field through the β 2 and β 3 factors respectively. These corrections are valid irrespective of the strength of the surface charge. However, for Newtonian fluids, the nonlinear relation between U and β may feature only for strongly charged surfaces, when the effect of surface conduction becomes important (Schnitzer & Yariv 2012;Schnitzer et al. 2013).
Fourth, it is well established that, for flow past spherical or cylindrical objects in polymeric liquids, when De ∼ O(1) or more, the polymers just behind the rear stagnation point witness strong extension rates. This leads to the formation of the so-called 'birefringent strands', wherein the extensional viscosity of the solution is greatly enhanced, even in dilute solutions -see the work of Harlen and coworkers (Harlen, Rallison & Chilcott 1990;Harlen, Hinch & Rallison 1992;Becherer, van Saarloos & Morozov 2009) for further details. This can profoundly alter the force on the particle and hence change its mobility significantly. Interestingly, formation of these strands has been found in both FENE-P (Harlen et al. 1990) as well as UCM type constitutive models (Becherer et al. 2009). However, in the present study, despite De ∼ O(1), the effective Deborah number is actually De eff ∼ O(ζ 0 ) 1, on account of the weakly charged particle (see § 3.5). As a result, the stretching of polymers remains limited and they are likely to retain their random coil shapes (Harlen et al. 1992) which hinders the formation of birefringent strands. Hence, it is expected that formation of such polymeric strands will not play a key role in the scenarios described within the scope of the present model. For completeness, in the supplementary material (see § S4 therein), we have shown an example of the variation in the polymeric stresses in the outer layer, around the particle.
Finally, it is intriguing to note the effect of non-uniform charge (i.e. a 1 ) density in combination with the medium's viscoelasticity on the particle's mobility. For a Newtonian fluid, non-zero a 1 always slows down (see the expression for U N ) a particle's velocity, which indicates that a non-uniformly charged particle would move slower than a uniformly charged one, at least when the surface charge is weak. However, in a viscoelastic medium, the sign of a 1 dictates whether the particle would speed up or slow down. For an Oldroyd-B fluid, λ 1 > λ 2 , as required for elongational viscosity (Bird et al. 1987). Therefore, from (4.9b), we note that for a 1 < 0, a non-uniformly charged particle would speed up in a viscoelastic medium; the reverse is true when a 1 > 0. On the other hand, the O(ζ 3 0 ) contribution always slows down the particle, regardless of the sign of a 1 . Perhaps the most striking feature of non-uniform surface charge and viscoelasticity is the appearance of the O(ζ 2 0 ) correction to the velocity, as evident from (4.11a). Notice that this term is absent when the particle is uniformly charged or the medium is Newtonian. Therefore, it is a unique outcome of the confluence of a non-uniform surface charge and the medium's viscoelasticity. One of its major consequences is that the presence of non-uniform charge augments the influence of the medium's viscoelasticity on the particle's electrophoretic velocity. Since this correction is proportional toζ 2 0 , it naturally leads to the breaking of fore-aft symmetry in electrophoresis of axisymmetric particles. More discussion is included in § 4.6.
Another promising way to probe the effects of non-uniform surface charge on the particle's mobility is to explore the influence of multipole moments of the surface charge, following the lead of Anderson and coworkers (Anderson 1985;Fair & Anderson 1989). They showed that (see (4.14)) in a Newtonian medium, subject to a uniform external electric field, a particle's velocity would depend on the zeroth (i.e. the average) and the quadrupole (i.e. second) moments of the 'zeta potential'. On the other hand, the dipole (i.e. the first) moment of φ P governs the angular velocity of the particle. These results were exact (at least in the thin EDL limit), while the particle velocity was essentially a linear combination of contributions from the zeroth and quadrupole moments of the potential. However, in a viscoelastic medium, because of the nonlinear constitutive relation, we expect that the various multipole moments will interact with each other to influence the particle velocity. At the same time, moments other than the quadrupole and the zeroth moment may also contribute to the particle's mobility. Here, we shall only consider contributions from up to the quadrupole moments of the surface potential.
Next, we shall try to express the velocity reported in (4.11a) using the imposed electric field, E ∞ = E 0êz and the multipole moments mentioned above, in tensor form. We further define, ν E = εψ c /η as the electrophoretic mobility. Then, the structure of the electrophoretic velocity in (4.11a) suggests that it may be rewritten order-wise in the following general form: where, Admittedly, only an axisymmetric analysis is not enough to obtain all of the coefficients b 11 , b 31 , . . . and so on. That said, based on Anderson's work (Anderson 1985) and (4.23d) we can conclude that b 11 = 1 2 . Further, comparing (4.23c) and (4.11a), we note that b 21 = − 9 5 (λ 0 /a)(λ 1 − λ 2 ). Also, taking the example of a uniformly charged particle, whence H 2 = H 1 = 0, it follows from (4.12) that b 33 = 24(λ 0 /a) 2 (λ 1 − λ 2 )( 20 819 5720 λ 1 − 23 8 λ 2 ). The coefficients b 31 and b 32 are also functions of λ 1 , λ 2 , a etc.; however, it is not possible to comment further on these coefficients based on the present analysis alone.
There are two key features of the equation (4.23), which are distinct from what is observed in a Newtonian medium. First, notice from (4.23c) that the dipole moment also contributes to the overall mobility in a viscoelastic fluid -this is what ultimately leads to symmetry breaking, as discussed in § 4.6. Second, it is quite obvious that, in a complex fluid, the electrophoretic velocity also has contributions from the interactions between the various multipole moments of the surface potential. The coefficients b 21 , b 31 , . . . etc. bear the signature of the medium's viscoelasticity and the associated interactions between various moments on the mobility of the particle. It is naturally expected that the higher-order moments may also influence the particle velocity, although they have not been considered here. Hence, the coefficients (such as b 21 ) are expected to have errors of O(||H 3 ||). Figure 4(a,b) demonstrates the contours of U/U N as a function of a 0 and a 1 , while values of other relevant parameters are given in the caption. In panel (a), contours for a 1 < 0 are shown, whereas, in (b), contours for a 1 > 0 are depicted. The assertions made in § 4.5 are clearly observed in these two panels. We note that, for a 1 < 0, the maximum increment in the velocity compared with a Newtonian medium is obtained when a 0 is small and a 1 is relatively large. It is further important to note that, for a 0 = 0, U = 0this is expected because a 0 = 0 indicates that the particle does not carry any net charge and hence its net velocity is constrained to be zero. When a 1 = 0, U < U De=0 always; this reduction occurs owing to the O(ζ 3 0 ) correction to U, which does not vanish even when the particle is uniformly charged. The O(ζ 3 0 ) correction underscores the effects of streamwise stretching of polymers inside the EDL due to the particle's curvature on its mobility -see § 3.2 for more discussion. Finally, it is interesting to note that a change in the sign of a 1 dictates whether the particle's velocity increases or decreases in a viscoelastic medium - this is in stark contrast to Newtonian fluids. This indicates that, in a viscoelastic fluid, a particle's mobility strongly depends on how the charge is distributed around its surface. It is to be noted that changing the sign of a 1 does not change U De=0 and hence the relative change in U between panels (a,b) indicates fore-aft symmetry breaking in a viscoelastic medium. Recall that this symmetry breaking occurs owing to the contributions from the dipole moment to the particle velocity, which comes at O(ζ 2 0 ). This symmetry breaking is a unique feature of the rheological paradigm considered here, since symmetry always prevails in Newtonian media in the Stokes flow regime (Leal 2007). In the next section, we further show that a similar feature also prevails in case of pure particle rotation in a viscoelastic medium, subject to an external electric field, provided that the surface charge is non-axisymmetric.

Results for electrophoretic velocity
5. Case study-II: electrophoretic rotation of a non-uniformly charged particle As a second example, we consider pure electrophoretic rotation of a particle, carrying non-uniform surface charge. In an effort to keep things simple and focused, we shall only consider the particle's angular velocity at a given instant, when the charge distribution on its surface is known. Further, in this section, we shall only report up to the first correction to the angular velocity occurring at O(ζ 2 0 ), caused by viscoelasticity. Any arbitrary non-axisymmetric surface chargeζ(μ, ϕ) may be represented using surface harmonics as (Happel & Brenner 2012)ζ(μ, ϕ) = ∞ n=0 S n (μ, ϕ), where S n (μ, ϕ) = n m=0 (A n cos mϕ +Â n sin mϕ)P m n (μ) is the surface harmonic of order n and P m n (x) is the associated Legendre polynomial (Happel & Brenner 2012) of orders n and m. In order to keep the algebra simple, we chooseζ(μ, ϕ) = sin θ sin ϕ + cos θ = P 0 1 (μ) − P 1 1 (μ) sin ϕ, which indicates that the particle does not carry any net charge and hence U = 0. As we will establish shortly, the particle will rotate with angular velocity Ω, which has an asymptotic expansion of the form Ω =ζ 0 Ω (1) +ζ 2 0 Ω (2) + . . . . Finally, we note that the fluid velocity vanishes far away from the particle (v → 0 as r → ∞).

The reciprocal theorem
The total stresses in a viscoelastic fluid may be expressed at any order ofζ 0 as (say, at kth order ofζ 0 ) τ , where σ k = −p k δ + 2D k is the Newtonian contribution to the stresses and τ (exc) k is the excess polymeric stress, which may be written as τ (k) ), see (4.5b). As a consequence, the momentum conservation equation in the outer layer at the kth order becomes ∇ · τ (tot) k = 0 (see (3.2c)), along with the continuity equation, ∇ · v k = 0. This may be rewritten as may be treated as the body force caused by the polymeric stresses. An auxiliary Newtonian flow with stressesσ and body forceb would satisfy (inertia is neglected) ∇ ·σ +b = 0 and ∇ ·v = 0,v being the auxiliary velocity field. Assuming that both viscoelastic flow and the auxiliary flow vanish sufficiently far away from the particle, the generalized reciprocal theorem may be written as (Masoud & Stone 2019;Li & Koch 2020) (5.1) In the above equation, the integration on the left-hand side is carried out over a spherical surface just outside the EDL, as mentioned in (2.4a,b) and the ones on the right-hand side are carried out over the entire outer region. Therefore, on S -see § 3.3. As the auxiliary flow, we choose the velocity field generated by an identical sphere rotating with angular velocityê (a constant unit vector) in a Newtonian fluid. As such,b = 0 andv = ∇ × (r −3ê · rr), which yields (Pozrikidis & Jankowski 1997)ê r ·σ = −3ê ×ê r andv =ê ×ê r on S p .

5.2.
The O(ζ 0 ) flow field Forζ(θ, ϕ) = cos θ + sin θ sin ϕ, at O(ζ 0 ), the velocity on the particle surface (r = 1) may be expressed as (using (3.14)) v 1 The resulting velocity field may be determined using Lamb's general solution (Happel & Brenner 2012), whereby v 1 may be written as where Λ −(n+1) , Θ −(n+1) and p −(n+1) are solid harmonics and have the following expressions: A m n cos mϕ +Â m n sin mϕ P m n (μ). (5.5) The rest of the solid harmonics also have analogous expressions and are not explicitly shown here for brevity. The velocity profile in (5.4) trivially satisfies the condition that the fluid velocity vanishes far away from the particle. Following the procedure outlined by Brenner (1964), we compute the quantities,ê r · v (1) S , −r∇ · v (1) S and r · ∇ × v (1) S , from which the spherical harmonics can be uniquely determined. It may be easily verified that The net torque on the particle at O(ζ 0 ) is given by (Happel & Brenner 2012) . Noting that N 1 = 0 (see (2.3a)), it follows that the leading-order (Newtonian) contribution to the angular velocity of the particle becomes It may be easily verified that (5.2), derived using the reciprocal theorem, also yields the same result, albeit without actually requiring the O(ζ 0 ) velocity field. Nevertheless, this velocity field will be essential for calculating the O(ζ 2 0 ) contribution to Ω, as we show next.
5.3. The O(ζ 2 0 ) correction The O(ζ 2 0 ) correction to the angular velocity of the particle may be directly evaluated using the reciprocal theorem, as outlined in (5.3a), since we have already determined the complete O(ζ 0 ) velocity field. The leading-order shear strain components (D 1 ) may be easily evaluated from the expression for v 1 in (5.7). Subsequently, its convected derivatives ( ∇ D 1 ) may be evaluated, by noting their components in spherical coordinate -see the book by Bird et al. (1987)  (5.9b) From the leading-order angular velocity, it may be noted that χ 1 = − 3 4 βμ cos ϕ and Γ 1 = − 3 4 β sin ϕ. This allows us to compute ω 2 and ω 3 , which are required for evaluating the O(ζ 2 0 ) slip velocity. It then follows from (3.18) that Notice that, at O(ζ 2 0 ), there is an azimuthal component of the slip velocity (v s,ϕ ), purely because of the viscoelasticity of the fluid and rotation at the leading order, which makes ω 3 / = 0. In other words, this component of the slip velocity actually originates from the interactions between particle rotation at O(ζ 0 ) and the excess polymeric stresses inside the EDL. This serves as an example of how the particle's angular velocity influences the Smoluchowski slip, as discussed in § 3.5. Expressions in (5.9) and the O(ζ 2 0 ) slip velocity 924 A41-36 in (5.10) can be inserted in (5.3a) to deduce the O(ζ 2 0 ) correction to the angular velocity of the particle, which reads The total angular velocity of the particle, up to O(ζ 2 0 ) is thus given by (5.12) 5.4. General discussion on particle rotation There are several interesting points to note from the angular velocity reported above. First, just like the translational velocity as discussed in § 4.2, the first viscoelastic contribution to the angular velocity also appears at O(ζ 2 0 ), when the particle bears a non-uniform surface charge. Interestingly, this contribution does not change upon reversing the direction of the electric field, which essentially indicates that the angular velocity may either increase or decrease depending on the direction of the imposed electric field. This is another example of symmetry breaking in viscoelastic media, which naturally comes about for non-uniformly charged surfaces, as already discussed in § § 4.5 and 4.6. Further, the sign of |Ω (2) | also depends on the choice of the constitutive model, given the charge distribution. For the Oldroyd-B model (where λ 1 > λ 2 ) or, the UCM model (λ 2 = 0), the O(ζ 2 0 ) viscoelastic correction increases Ω, while for a second-order model with vanishing second normal stress coefficient (λ 1 = 0), the rotation actually slows down.
In a Newtonian medium Ω (2) = 0 identically, even if the particle carries non-uniform and non-axisymmetric surface charge. The next correction in such cases would come at O(ζ 3 0 ). Interestingly, it is straightforward to verify (using the procedure outlined in § 5.3) that, if one choosesζ(θ, ϕ) = sin θ sin ϕ, Ω (2) = 0 identically and then the first viscoelastic correction comes at O(ζ 3 0 ). Therefore, the dominant viscoelastic correction to particle's rotation strongly depends on how the charge is distributed across its surface. Finally, we note that particle rotation will alter its surface charge distribution with respect to a non-rotating frame and the mathematical form ofζ(θ, ϕ) would continually change and thus the Ω reported in (5.12) is only valid at the instant whenζ(θ, ϕ) = sin θ sin ϕ + cos θ. Although we do not actively pursue how the surface charge evolves and how that change in turn alters the particle's motion here, the rate of change ofζ may be determined from equations (S2) in the supplementary material.

Experimental perspectives
In this section, we briefly discuss a potential experimental set-up that can be used to test some of the key predictions of our analysis. In this regard, it should be noted that experimental studies on electrophoresis in complex fluids are scarce (Lu et al. 2014(Lu et al. , 2015Malekanfard et al. 2019), while any focused investigation on the effects of viscoelasticity on single particle electrophoresis is completely missing, to the best of our knowledge. Figure 5 depicts a possible experimental set-up, which is very similar to those used by Malekanfard et al. (2019) as well as Liang et al. (2010). The set-up consists of two reservoirs connected by a conduit, containing a viscoelastic liquid. Several polymer solutions such as aqueous poly-ethylene oxide (PEO) solution with added electrolytes (such as NaCl/HCl) may be used as the viscoelastic medium (Ardekani et al. 2009). . Buffer solutions may be added to the viscoelastic medium to make the particle concentration lower, so that the conditions approach that of a single particle electrophoresis. A DC electric field of O(10 kV m −1 ) can be applied by placing electrodes inside the reservoirs and an inverted microscope may be used to capture particle motion (Liang et al. 2010). The channel should be large enough so that its walls do not influence the particle velocity (see Liang et al. 2010;Malekanfard et al. 2019).
We will end this section by noting the specific predictions which may be compared with experimental results. First, it is expected that a uniformly charged particle would move slower as compared with a Newtonian fluid with comparable viscosity. Second, particles with different sizes should move at different velocities; in particular, for a uniformly charged particle, larger particles will move faster than smaller ones. For weakly and non-uniformly charged particles, the velocity may either increase or decrease in a viscoelastic fluid, depending on how the charge is distributed around the surface; however, this change is expected to be far more pronounced in a viscoelastic medium. Finally, our analysis predicts breaking of fore-aft symmetry for non-uniformly charged particles. This can be tested by noting whether the particle's mobility changes upon changing the direction of the imposed DC field in the channel.

Conclusions
The present study delves into the electrokinetics around a weakly charged spherical particle in a viscoelastic fluid, in the limit of thin EDLs. Analytically, we apply singular perturbation to shed light on the dynamics within the EDL and subsequently use a regular power series expansion in the characteristic potentialζ 0 , which is assumed to be small and hence effectively renders the resulting flow weakly viscoelastic. We develop a generic framework to evaluate the modified Smoluchowski slip velocity for an arbitrarily charged particle, which might undergo arbitrary motion. We further assess the accuracy of the results derived using the Oldroyd-B model, by comparing them with numerically computed results for the FENE-P model, which is one of the simplest nonlinear viscoelastic models that remains valid in the presence of strong elongational stresses. Our analysis revealed that the modified slip velocity shows a strong dependence on the variability of the surface charge, as well as the particle's curvature and angular velocity. These features are unique to the complex fluid considered here and are not witnessed in Newtonian fluids; most of them occur because the normal stress components within the EDL become asymptotically large as compared with the shear stresses and hence they play a key role in governing the momentum transport. To demonstrate an application of the modified Smoluchowski slip thus derived, we separately considered two specific cases of electrophoretic translation and rotation of a particle carrying non-uniform surface charge.
Our analysis reveals that the electrophoretic velocity in a viscoelastic fluid depends on the particle's curvature, even in the regime of thin EDLs and weak surface chargesin sharp contrast to Newtonian fluids. As a result, particles of different sizes move with different velocities, when acted upon by the same external electric field. We further demonstrated that the particle's radius, the medium's viscoelasticity as well as the extent of non-uniformity of the surface charge together dictate the particle's velocity through the fluid. We subsequently explore how the multipole moments of the particle's potential (φ P ) influence its movement. To this end, we establish that, because of the nonlinear nature of the constitutive model of viscoelastic fluids, the multipole moments interact between each other to alter the particle's velocity. In particular, the dipole moment results in an O(ζ 2 0 ) contribution to the mobility, which is absent in a Newtonian medium. As a result, the particle's velocity may either increase or decrease, depending on the nature of surface charge distribution; the same is not true for a Newtonian fluid. At the same time, the fore-aft symmetry, which prevails in Stokes flows, is broken in a viscoelastic medium, because of the O(ζ 2 0 ) contribution from the dipole moment of φ P . As a consequence, the particle's radius plays a more important role in determining its velocity, when it carries a non-uniform surface charge. The mobility of a uniformly charged particle may be derived as a special case of the analysis and recovers the result that a viscoelastic medium will trivially slow down the particle under such conditions.
Comparison with numerical simulations of the FENE-P model reveals that the predictions using the Oldroyd-B model are reasonably accurate, so long as the surface potential is sufficiently weak and the maximum permissible spring length in a bead-spring model of polymers is large. This ensures that shear thinning behaviour remains subdominant, which enhances the applicability of the Oldroyd-B model. The comparison between the two models also works well in the limit of low polymer viscosity (C 1). Beyond a critical value ofζ 0 , however, the Oldroyd-B model underpredicts the electrophoretic mobility of the particle, as the shear thinning behaviour starts to dominate. The accuracy of the Oldroyd-B model, when compared with the FENE-P model is, however, far more sensitive to the strength of the surface charge, as compared withrepresentative of the maximum permissible length of the polymers.
The second case study on pure electrophoretic rotation of a non-axisymmetrically charged particle was carried out using a combination of Lamb's general solution and the generalized reciprocal theorem. It was shown that the medium's viscoelasticity may alter the angular velocity of the particle in a very similar way as observed for the translational velocity. For instance, the angular velocity may both increase or decrease because of the medium's viscoelasticity, depending on the direction of the applied electric field as well as the fluid properties. This brings out another example of fore-aft symmetry breaking in polymeric liquids, driven by non-uniformity of the surface charge on the particle.
In summary, the present analysis provides a platform to improve the state of the art understanding of the transport of charged species in biologically relevant fluidic media, offering a new design paradigm for in vitro analytics. Further studies, however, are deemed essential to focus on some more exclusive aspects of the fluid rheology that do not fall in the purview of the constitutive model that has been considered in the present work, as well as the implications of temporal variations in the surface charge due to particle rotation.