2 Experiments

In our experiments, we use a 2 cm long Polydimethylsiloxane (PDMS) microchannel with a rectangular cross-section of $66\times 54~\unicode[STIX]{x03BC}\text{m}^{2}$. The electric fields were generated using a high-voltage direct current power supply using a wire connection to electrodes at the channel inlet and outlet. We dissolved polyethylene oxide (PEO) powder (molecular weight $=2\times 10^{6}~\text{g}~\text{mol}^{-1}$, Sigma-Aldrich) into a buffer solution (0.01 mM phosphate buffer mixed with 0.5 % Tween 20, Fisher Scientific) at a concentration of 250 ppm. A dilute suspension (<0.1 % in volume fraction) of polystyrene spheres of 2.2 μm diameter is prepared in the PEO solution and is pumped at $25~\unicode[STIX]{x03BC}\text{l}~\text{h}^{-1}$. The particle zeta potential was measured to be $-83~\text{mV}$. Particle motion was visualized at the outlet of the microchannel through an inverted microscope (Nikon Eclipse TE2000U) equipped with a CCD camera (Nikon DS-Qi1MC). Digital videos were recorded at a rate of 15 frames per second, from which the superimposed images were obtained and further processed in the Nikon imaging software (NIS-Elements AR 3.22).

Figure 1(a) shows the particles at the inlet and outlet of the channel when no electric field is applied, i.e. the focusing is purely viscoelastic (Karnis & Mason Reference Karnis and Mason1966; Gauthier et al. Reference Gauthier, Goldsmith and Mason1971; Ho & Leal Reference Ho and Leal1976). The particles remain almost uniformly distributed at the channel exit. This weak focusing is depicted in figure 1(d), which shows the change in particle stream width (lateral width of the particle distribution) along the channel length. In figure 1(b) an application of an electric field parallel to the flow shows that a significant fraction of the particles are focused at the centreline. Figure 1(d) shows that the introduction of an electric field results in substantial reduction in particle stream width. Upon reversing the direction of the electric field (see figure 1c), the particles focus near the walls. These results demonstrate the potential to efficiently separate particles and cells based on their intrinsic electrical properties. To extrapolate these observations towards potential applications, it is pivotal to understand the role of electrical properties and flow profile characteristics.

Figure 1.

Superimposed experimental images showing the focusing of polystyrene particles (negative zeta potential), at the channel inlet and outlet, for: (a) no external electric field; (b) electric field in the direction of the flow. (c) Particle focusing when electric field is applied opposite to the flow. (d) Particle stream width for (a,b). Scale bar represents 100 μm. Parameters: ${\mathcal{E}}_{\infty }^{\ast }=300~\text{V}~\text{cm}^{-1}$; $C_{PEO}^{\ast }=250~\text{ppm}$; $a^{\ast }=1.1~\unicode[STIX]{x03BC}\text{m}$; $l^{\ast }=66~\unicode[STIX]{x03BC}\text{m}$; $U_{max}^{\ast }=2.9~\text{mm}~\text{s}^{-1}$; $\unicode[STIX]{x1D707}^{\ast }=1.3~\text{m}~\text{Pa}~\text{s}$; $\unicode[STIX]{x1D701}_{p}^{\ast }=-83~\text{mV}$; $\unicode[STIX]{x1D701}_{w}^{\ast }=-120~\text{mV}$; $Re_{p}=1.6\times 10^{-4}$. Experimental movies corresponding to the above superimposed images can be found as movies 1–6 in the supplementary material, available at https://doi.org/10.1017/jfm.2020.397.

In what follows, we formulate a mathematical framework to predict, quantify and draw physical insights from the migration behaviour observed in the experiments. Since the particle suspension in the experimental study is dilute, we analyse the migration of a single particle. We model the fluid medium using a second-order fluid (SOF) model (Bird et al. Reference Bird, Armstrong and Hassager1987), as we are only interested in capturing the effect of small deviations from the Newtonian behaviour. Although the SOF model is only valid for slow and slowly varying (thus nearly Newtonian) flows, it is an asymptotic approximation for most viscoelastic flows (Ho & Leal Reference Ho and Leal1976). Furthermore, it has been shown in the past that dilute PEO solutions ($C_{PEO}^{\ast }<500~\text{ppm}$) exhibit pure elastic behaviour (i.e. negligible shear thinning) (Rodd et al. Reference Rodd, Scott, Boger, Cooper-White and McKinley2005; Yang et al. Reference Yang, Kim, Lee, Lee and Kim2011; Kim et al. Reference Kim, Ahn, Lee and Kim2012). Hence employing the SOF model can help us capture the dynamics in slow Boger fluid flows, both qualitatively and quantitatively.

3 Mathematical model

Figure 2 shows the schematic of a neutrally buoyant electrophoretic sphere suspended in a pressure driven flow of SOF, subjected to a parallel electric field (${\mathcal{E}}_{\infty }^{\ast }$). The fluid moves with a maximum velocity of $U_{max}^{\ast }$ at the centreline. The frame of reference is placed at the centre of the particle. Typical PDMS microchannel walls and polystyrene particles usually acquire a surface charge (formally quantified as zeta-potential $\unicode[STIX]{x1D701}^{\ast }$) when in contact with an electrolytic solution (Sze et al. Reference Sze, Erickson, Ren and Li2003), which gives rise to an electrical double layer (EDL) (Anderson Reference Anderson1989). In comparison to the channel and particle size, the EDL can be assumed to be asymptotically thin, i.e. $\unicode[STIX]{x1D706}_{D}^{\ast }/a^{\ast }\ll 1$ (where $\unicode[STIX]{x1D706}_{D}^{\ast }$ is the EDL thickness). Under the influence of an external electric field, the ions inside the EDL slip in the tangential direction, generating an electro-osmotic slip. For low zeta potentials ($\unicode[STIX]{x1D701}_{p}^{\ast }\sim 25~\text{mV}$), the slip at the particle surface can be represented by Smoluchowski’s relation (Keh & Anderson Reference Keh and Anderson1985; Anderson Reference Anderson1989). Recently, Ghosh, Chaudhury & Chakraborty (Reference Ghosh, Chaudhury and Chakraborty2016) showed that the surface slip in the SOF does not get modified for a surface with constant zeta potential. Hence we employ the Smoluchowski slip in this work.

Figure 2.

Schematic of the problem formulation. The origin of the coordinate system is at the centre of the particle. Here, $y^{\ast }$ points in the direction of vorticity of the pressure driven flow. The dashed layer around the particle and walls depicts the outer edge of the double layer. Here, * denotes the dimensional variables.

3.1 Governing equations

We represent the system in terms of dimensionless variables, where the coordinates, electrostatic potential, velocity and pressure are non-dimensionalized using $a^{\ast }$, $|{\mathcal{E}}_{\infty }^{\ast }|a^{\ast }$, $\unicode[STIX]{x1D705}U_{max}^{\ast }$ and $\unicode[STIX]{x1D707}^{\ast }\unicode[STIX]{x1D705}U_{max}^{\ast }/a^{\ast }$, respectively ($\unicode[STIX]{x1D707}^{\ast }$ being the dynamic viscosity). Here $\unicode[STIX]{x1D705}$ is the ratio of particle radius ($a^{\ast }$) to the channel width ($l^{\ast }$). To obtain analytical results, we will later restrict our analysis to small particles ($\unicode[STIX]{x1D705}\ll 1$), which is consistent with the experimental conditions.

The potential ($\unicode[STIX]{x1D719}$) in the electroneutral bulk region is governed by the Laplace equation

(3.1)$$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}=0.\end{eqnarray}$$

The boundary conditions at the particle and wall surface (i.e. slip plane of the EDL) are governed by the no-flux condition

(3.2a,b)$$\begin{eqnarray}\boldsymbol{e}_{r}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=0\quad \text{at}~r=1\quad \text{and}\quad \boldsymbol{e}_{z}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=0\quad \text{at}~z=-s/\unicode[STIX]{x1D705}~\text{and}~z=(1-s)/\unicode[STIX]{x1D705},\end{eqnarray}$$

where $s\equiv d^{\ast }/l^{\ast }$ (as defined in figure 2). The potential is undisturbed ($\unicode[STIX]{x1D719}_{\infty }\sim -x$) far away from the particle

(3.3)$$\begin{eqnarray}\unicode[STIX]{x1D719}\rightarrow \unicode[STIX]{x1D719}_{\infty }\quad \text{as}~r\rightarrow \infty .\end{eqnarray}$$

The hydrodynamics is governed by the continuity and the Navier–Stokes equations

(3.4a,b)$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{U}=0\quad \text{and}\quad Re_{p}(\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{U})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64F},\end{eqnarray}$$

where $Re_{p}$ is particle Reynolds number defined as $\unicode[STIX]{x1D70C}^{\ast }U_{max}^{\ast }\unicode[STIX]{x1D705}a^{\ast }/\unicode[STIX]{x1D707}^{\ast }$ ($\unicode[STIX]{x1D70C}^{\ast }$ being the fluid density). We model the hydrodynamic stress using the second-order-fluid model (Bird et al. Reference Bird, Armstrong and Hassager1987)

(3.5)$$\begin{eqnarray}\unicode[STIX]{x1D64F}=-P\unicode[STIX]{x1D644}+\unicode[STIX]{x1D640}^{(1)}+De\unicode[STIX]{x1D64E}.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D640}^{(1)}$ represents the rate of strain tensor ($\unicode[STIX]{x1D735}\boldsymbol{U}+\unicode[STIX]{x1D735}\boldsymbol{U}^{\dagger }$), $\dagger$ denotes the transpose, $\unicode[STIX]{x1D64E}$ represents the polymeric stress

(3.6)$$\begin{eqnarray}\unicode[STIX]{x1D64E}=\unicode[STIX]{x1D640}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D640}^{(1)}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D640}^{(2)}.\end{eqnarray}$$

$\unicode[STIX]{x1D640}^{(2)}$ is the steady-state covariant convected derivative of $\unicode[STIX]{x1D640}^{(1)}$ (also known as the Rivlin–Ericksen tensor)

(3.7)$$\begin{eqnarray}\unicode[STIX]{x1D640}^{(2)}=\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D640}^{(1)}+\unicode[STIX]{x1D640}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{U}^{\dagger }+\unicode[STIX]{x1D735}\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D640}^{(1)}.\end{eqnarray}$$

Here, $De$ (Deborah number) is the ratio of viscoelastic time scale ($(\unicode[STIX]{x1D6F9}_{1}^{\ast }+\unicode[STIX]{x1D6F9}_{2}^{\ast })/\unicode[STIX]{x1D707}^{\ast }$) to the time scale corresponding to the average shear in the background flow ($l^{\ast }/U_{max}^{\ast }$), $\unicode[STIX]{x1D6FF}$ is a viscometric parameter, $\unicode[STIX]{x1D6F9}_{1}^{\ast }$ and $\unicode[STIX]{x1D6F9}_{2}^{\ast }$ are the steady shear viscometric coefficients (Bird et al. Reference Bird, Armstrong and Hassager1987)

(3.8a-d)$$\begin{eqnarray}De=\frac{(\unicode[STIX]{x1D6F9}_{1}^{\ast }+\unicode[STIX]{x1D6F9}_{2}^{\ast })U_{max}^{\ast }}{\unicode[STIX]{x1D707}^{\ast }l^{\ast }},\quad \unicode[STIX]{x1D6FF}=\frac{-\unicode[STIX]{x1D6F9}_{1}^{\ast }}{2(\unicode[STIX]{x1D6F9}_{1}^{\ast }+\unicode[STIX]{x1D6F9}_{2}^{\ast })},\quad \unicode[STIX]{x1D6F9}_{1}^{\ast }=\frac{T_{xx}^{\ast }-T_{zz}^{\ast }}{(\unicode[STIX]{x2202}U_{x}^{\ast }/\unicode[STIX]{x2202}z^{\ast })^{2}},\quad \unicode[STIX]{x1D6F9}_{2}^{\ast }=\frac{T_{yy}^{\ast }-T_{zz}^{\ast }}{(\unicode[STIX]{x2202}U_{x}^{\ast }/\unicode[STIX]{x2202}z^{\ast })^{2}}.\end{eqnarray}$$

Since the migration time scale is much larger than the convective time scale ($a^{\ast }/\unicode[STIX]{x1D705}U_{max}^{\ast }$), we analyse the migration in quasi-steady-state (Becker, McKinley & Stone Reference Becker, McKinley and Stone1996; Choudhary et al. Reference Choudhary, Renganathan and Pushpavanam2019). Also, we neglect inertial effects (i.e. $Re_{p}\ll De$) which leaves the polymeric stress as the sole nonlinearity in the system.

The sphere moves with a translational ($\boldsymbol{U}_{s}$) and rotational ($\unicode[STIX]{x1D734}_{s}$) velocity, which is to be determined from force and torque balance. Since the domain is infinite in the $y$-direction, we can consider $\boldsymbol{U}_{s}=U_{sx}\boldsymbol{e}_{x}+U_{sz}\boldsymbol{e}_{z}$ and $\unicode[STIX]{x1D734}_{s}=\unicode[STIX]{x1D6FA}_{sy}\boldsymbol{e}_{y}$. To account for the electrokinetic effects in the EDL, we employ the slip-boundary condition (Smoluchowski Reference Smoluchowski1903) and impose at the particle surface

(3.9)$$\begin{eqnarray}\boldsymbol{U}=\unicode[STIX]{x1D734}_{s}\times \boldsymbol{r}+Ha\unicode[STIX]{x1D701}_{p}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\quad \text{at}~r=1.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D701}_{p}$ is the dimensionless particle surface zeta potential ($\unicode[STIX]{x1D701}_{p}^{\ast }/|{\mathcal{E}}_{\infty }^{\ast }|a^{\ast }$) and $Ha$ is a Hartmann-type number, $Ha=\unicode[STIX]{x1D700}^{\ast }{{\mathcal{E}}^{\ast }}_{\infty }^{2}a^{\ast }/4\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}^{\ast }\unicode[STIX]{x1D705}U_{max}^{\ast }$, which denotes the ratio of electrical energy density to shear stress. Here, $\unicode[STIX]{x1D700}^{\ast }$ is the electrical permittivity of the medium. Similarly, the flow satisfies the Smoluchowski slip condition at the walls

(3.10)$$\begin{eqnarray}\boldsymbol{U}=Ha\unicode[STIX]{x1D701}_{w}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}-\boldsymbol{U}_{s}\quad \text{at walls}.\end{eqnarray}$$

As the particle is freely suspended, the translational ($\boldsymbol{U}_{s}$) and rotational velocity ($\unicode[STIX]{x1D734}_{s}$) of the particle are obtained using force-free and torque-free conditions: the total force and torque on the particle must be zero. Apart from the hydrodynamic stress, we also have to account for the Maxwell stress in the calculation of total force and torque. Maxwell stress arises due to a lateral asymmetry in the electrostatic potential distribution around the particle due to the presence of the wall. The Maxwell stress is defined as

(3.11)$$\begin{eqnarray}\unicode[STIX]{x1D648}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}-{\textstyle \frac{1}{2}}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719})\unicode[STIX]{x1D644}.\end{eqnarray}$$

It should be noted that there is no contribution from Maxwell stress to the fluid momentum equation (3.4) because $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D648}=0$ in the electroneutral fluid (Yariv Reference Yariv2006). The total force-free and torque-free conditions for particle are given by

(3.12a)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}=\boldsymbol{F}_{H}+\boldsymbol{F}_{M}=\int _{S_{p}}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}\,\text{d}S+4\unicode[STIX]{x03C0}Ha\int _{S_{p}}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D648}\,\text{d}S=\mathbf{0}, & \displaystyle\end{eqnarray}$$

(3.12b)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{L}=\boldsymbol{L}_{H}+\boldsymbol{L}_{M}=\int _{S_{p}}\boldsymbol{n}\times (\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D64F})\,\text{d}S+4\unicode[STIX]{x03C0}Ha\int _{S_{p}}\boldsymbol{n}\times (\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D648})\,\text{d}S=\mathbf{0}. & \displaystyle\end{eqnarray}$$

Here, subscript $H$ and $M$ represent contributions from hydrodynamic and Maxwell stress, respectively.

3.1.1 Undisturbed flow

In the absence of particle, the flow profile for a fully developed pressure driven second-order fluid flow (in the frame of reference of the particle moving with $\boldsymbol{U}_{s}$) is

(3.13)$$\begin{eqnarray}\boldsymbol{U}_{\infty }=(\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FD}z+\unicode[STIX]{x1D6FE}z^{2})\boldsymbol{e}_{x}+Ha\unicode[STIX]{x1D701}_{w}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{\infty }-\boldsymbol{U}_{s}.\end{eqnarray}$$

This is identical to that obtained for a Newtonian fluid (Choudhary et al. Reference Choudhary, Renganathan and Pushpavanam2019). However, the pressure exhibits a variation in the lateral direction for SOF

(3.14)$$\begin{eqnarray}P_{\infty }=2\unicode[STIX]{x1D6FE}x+4De(1+2\unicode[STIX]{x1D6FF})(\unicode[STIX]{x1D6FD}z+\unicode[STIX]{x1D6FE}z^{2})+{\mathcal{I}},\end{eqnarray}$$

where ${\mathcal{I}}$ is an integration constant. The constants $\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FE}$ are

(3.15a-c)$$\begin{eqnarray}\unicode[STIX]{x1D6FC}=4s(1-s)/\unicode[STIX]{x1D705},\quad \unicode[STIX]{x1D6FD}=4(1-2s),\quad \unicode[STIX]{x1D6FE}=-4\unicode[STIX]{x1D705},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FE}$ represent the shear and curvature of the background flow.

3.2 Equations governing the disturbance

We analyse the system in terms of disturbance variables (i.e. deviation from the undisturbed flow). For the different variables, these are defined as: potential ($\unicode[STIX]{x1D713}=\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{\infty }$), velocity ($\boldsymbol{v}=\boldsymbol{U}-\boldsymbol{U}_{\infty }$) and pressure ($p=P-P_{\infty }$). Since the bulk fluid medium is electroneutral, the disturbance potential is also governed by the Laplace equation

(3.16)$$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}=0.\end{eqnarray}$$

Far away from the particle, $\unicode[STIX]{x1D713}$ decays to zero. Using the definition of disturbance potential in (3.2),

(3.17a,b)$$\begin{eqnarray}\boldsymbol{e}_{\boldsymbol{r}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}=-\boldsymbol{e}_{\boldsymbol{r}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{\infty }\quad \text{at}~r=1\quad \text{and}\quad \boldsymbol{e}_{\boldsymbol{z}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}=0\quad \text{at walls}.\end{eqnarray}$$

The disturbance flow is governed by continuity and Cauchy’s momentum equation

(3.18a,b)$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0\quad \text{and}\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}_{H}=\mathbf{0}.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D748}_{H}$ is the hydrodynamic disturbance stress tensor which contains the Newtonian and polymeric contributions ($\unicode[STIX]{x1D748}_{H}=-p\unicode[STIX]{x1D644}+\unicode[STIX]{x1D65A}^{(1)}+De\,\unicode[STIX]{x1D668}$). Here, $\unicode[STIX]{x1D644}$ is the identity tensor, $\unicode[STIX]{x1D65A}^{(1)}$ is the rate of strain tensor for the disturbance flow ($\unicode[STIX]{x1D735}\boldsymbol{v}+\unicode[STIX]{x1D735}\boldsymbol{v}^{\dagger }$) and $\unicode[STIX]{x1D668}$ is the disturbance polymeric stress tensor

(3.19)$$\begin{eqnarray}\unicode[STIX]{x1D668}=\unicode[STIX]{x1D65A}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D65A}^{(1)}+\unicode[STIX]{x1D66C}^{(1)}+\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D65A}^{(2)}+\unicode[STIX]{x1D66C}^{(2)}),\end{eqnarray}$$

(3.20)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D66C}^{(1)}=\unicode[STIX]{x1D640}_{\infty }^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D65A}^{(1)}+\unicode[STIX]{x1D65A}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D640}_{\infty }^{(1)},\quad \unicode[STIX]{x1D65A}^{(2)}=\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D65A}^{(1)}+\unicode[STIX]{x1D65A}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}^{\dagger }+\unicode[STIX]{x1D735}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D65A}^{(1)},\\ \displaystyle \unicode[STIX]{x1D66C}^{(2)}=\boldsymbol{U}_{\infty }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D65A}^{(1)}+\unicode[STIX]{x1D65A}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\boldsymbol{U}_{\infty }}^{\dagger }+\unicode[STIX]{x1D735}\boldsymbol{U}_{\infty }\boldsymbol{\cdot }\unicode[STIX]{x1D65A}^{(1)}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D640}_{\infty }^{(1)}+\unicode[STIX]{x1D640}_{\infty }^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}^{\dagger }+\unicode[STIX]{x1D735}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D640}_{\infty }^{(1)}.\end{array}\right\}\end{eqnarray}$$

Here, $\unicode[STIX]{x1D640}_{\infty }^{(1)}$ is the rate of strain tensor for the undisturbed flow, $\unicode[STIX]{x1D65A}^{(2)}$ is the disturbance Rivlin–Eriksen (RE) tensor, $\unicode[STIX]{x1D66C}^{(1)}$ is the ‘interaction tensor’ (arising from the interaction between background flow and disturbance flow) and $\unicode[STIX]{x1D66C}^{(2)}$ is the RE interaction tensor. Using (3.9) and (3.10), we find that the disturbance flow is subject to the following boundary conditions:

(3.21a)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}=\unicode[STIX]{x1D734}_{s}\times \boldsymbol{r}+Ha\unicode[STIX]{x1D701}_{p}\unicode[STIX]{x1D735}(\unicode[STIX]{x1D713}+\unicode[STIX]{x1D719}_{\infty })-\boldsymbol{U}_{\infty }\quad \text{at}~r=1, & \displaystyle\end{eqnarray}$$

(3.21b)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}=Ha\unicode[STIX]{x1D701}_{w}\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\quad \text{at walls}, & \displaystyle\end{eqnarray}$$

(3.21c)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}\rightarrow \mathbf{0}\quad \text{as}~r\rightarrow \infty . & \displaystyle\end{eqnarray}$$

The total Maxwell stress can be decomposed into four components as

(3.22)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D648} & = & \displaystyle \left(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{\infty }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{\infty }-{\textstyle \frac{1}{2}}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{\infty }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{\infty })\unicode[STIX]{x1D644}\right)+\left(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{\infty }-{\textstyle \frac{1}{2}}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{\infty })\unicode[STIX]{x1D644}\right)\nonumber\\ \displaystyle & & \displaystyle +\left(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{\infty }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}-{\textstyle \frac{1}{2}}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{\infty }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713})\unicode[STIX]{x1D644}\right)+\left(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}-{\textstyle \frac{1}{2}}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713})\unicode[STIX]{x1D644}\right).\end{eqnarray}$$

The first term arises from the undisturbed component, the second and third terms denote the contribution from the interaction between applied ($\unicode[STIX]{x1D719}_{\infty }$) and disturbed potential field ($\unicode[STIX]{x1D713}$) and the fourth term denotes the disturbance Maxwell stress. Yariv (Reference Yariv2006) showed that the first three components do not result in a net force or torque and the leading-order contribution arises from the fourth term. Since $\unicode[STIX]{x1D648}$ is divergence free in the fluid domain, only the terms (in (3.22)) decaying as $O(r^{-2})$ contribute to the force (Guazzelli & Morris Reference Guazzelli and Morris2011). Only the disturbance Maxwell stress generates $O(r^{-2})$ terms and provides the leading-order contribution (Yariv Reference Yariv2006, p. 3). Thus, the disturbance Maxwell stress is denoted as

(3.23)$$\begin{eqnarray}\unicode[STIX]{x1D748}_{M}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}-{\textstyle \frac{1}{2}}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713})\unicode[STIX]{x1D644}.\end{eqnarray}$$

Similarly, the undisturbed hydrodynamic field cannot generate a net force or torque on the particle. Hence, the force-free and torque-free conditions are given by

(3.24a)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}=\boldsymbol{F}_{H}+\boldsymbol{F}_{M}=\int _{S_{p}}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D748}_{H}\,\text{d}S+4\unicode[STIX]{x03C0}Ha\int _{S_{p}}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D748}_{M}\,\text{d}S=\mathbf{0}, & \displaystyle\end{eqnarray}$$

(3.24b)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{L}=\boldsymbol{L}_{H}+\boldsymbol{L}_{M}=\int _{S_{p}}\boldsymbol{n}\times (\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D748}_{H})\,\text{d}S+4\unicode[STIX]{x03C0}Ha\int _{S_{p}}\boldsymbol{n}\times (\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D748}_{M})\,\text{d}S=\mathbf{0}. & \displaystyle\end{eqnarray}$$

3.3 Migration due to Maxwell stress

Here we evaluate the force and lateral migration arising from the Maxwell stress (3.23) alone. Since the potential distribution is not affected by the fluid rheology, the Maxwell force on the particle in a second-order fluid is identical to that in a Newtonian fluid (Yariv Reference Yariv2006). Following our previous analysis (Choudhary et al. Reference Choudhary, Renganathan and Pushpavanam2019), we use the method of reflections to determine the potential distribution in the particle–wall configuration (Brenner Reference Brenner1962). Since our aim is to obtain qualitative insights and analytical results, we restrict the analysis to small particles i.e. $\unicode[STIX]{x1D705}\ll 1$. Assuming the small particle is not too close to the walls (i.e. $s\gg \unicode[STIX]{x1D705}$), the disturbance potential is sought as successive reflections, $\unicode[STIX]{x1D713}=_{(1)}\unicode[STIX]{x1D713}+_{(2)}\unicode[STIX]{x1D713}+\cdots$. Here, $_{(i)}\unicode[STIX]{x1D713}$ represents the $i\text{th}$ reflection, where the odd reflections satisfy boundary conditions at the particle surface and the even reflections satisfy the wall boundary conditions. Accounting for each successive pair of reflections increases the accuracy by $O(\unicode[STIX]{x1D705})$ (Happel & Brenner Reference Happel and Brenner2012, p. 286). Since $\unicode[STIX]{x1D705}\ll 1$, we take only the first two reflections of disturbance potential into account. Evaluation of $\boldsymbol{F}_{M}$ and $\boldsymbol{L}_{M}$ (detailed in § A.2) yields

(3.25a,b)$$\begin{eqnarray}\boldsymbol{F}_{M}=4\unicode[STIX]{x03C0}Ha\left(\frac{3\unicode[STIX]{x03C0}}{16}\unicode[STIX]{x1D705}^{4}(\mathfrak{Z}(4,s)-\mathfrak{Z}(4,1-s))\right)\boldsymbol{e}_{z}+O(\unicode[STIX]{x1D705}^{7})\quad \text{and}\quad \boldsymbol{L}_{M}=\mathbf{0}.\end{eqnarray}$$

The Maxwell force acts along the cross-stream direction alone and is wall-repulsive in nature; the Maxwell torque is zero. The contribution of Maxwell stresses to the migration velocity is

(3.26)$$\begin{eqnarray}U_{mig}^{M}\approx \frac{4\unicode[STIX]{x03C0}Ha}{6\unicode[STIX]{x03C0}(1+O(\unicode[STIX]{x1D705}))}\left(\frac{3\unicode[STIX]{x03C0}}{16}\unicode[STIX]{x1D705}^{4}(\mathfrak{Z}(4,s)-\mathfrak{Z}(4,1-s))\right).\end{eqnarray}$$

Here, superscript $M$ shows that the contribution is from Maxwell stress, $\mathfrak{Z}$ is the generalized Riemann zeta function and $6\unicode[STIX]{x03C0}(1+O(\unicode[STIX]{x1D705}))$ is the viscous resistance which includes wall correction at the first order (Brenner Reference Brenner1961, p. 246) (the details of evaluation through the method of reflections are provided in the supplementary material). Later we show that $U_{mig}^{M}$ decays rapidly away from the walls and is negligible in the bulk of the channel, where migration due to hydrodynamic stress dominates.

3.4 Migration due to hydrodynamic stress

3.4.1 Perturbation expansion

We find the migration at $O(De)$ using a regular perturbation expansion. For small values of $De$, the disturbance field variables are expanded as

(3.27)$$\begin{eqnarray}\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{(0)}+De\,\unicode[STIX]{x1D709}_{(1)}+\cdots .\end{eqnarray}$$

Here, $\unicode[STIX]{x1D709}$ is a generic field variable which represents velocity ($\boldsymbol{v}$), pressure ($p$), translational ($\boldsymbol{U}_{s}$) and angular velocity ($\unicode[STIX]{x1D734}_{s}$). Since the electrostatic potential is decoupled from the hydrodynamics, the variable $\unicode[STIX]{x1D713}$ is not expanded as in (3.27). We substitute (3.27) in the equations governing the hydrodynamics (3.18)–(3.21). We obtain the problem at $O(1)$ as

(3.28)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}_{(0)}=0,\\ \unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{(0)}-\unicode[STIX]{x1D735}p_{(0)}=\mathbf{0},\\ \boldsymbol{v}_{(0)}={\unicode[STIX]{x1D734}_{s}}_{(0)}\times \boldsymbol{r}+Ha\unicode[STIX]{x1D701}_{p}\unicode[STIX]{x1D735}(\unicode[STIX]{x1D713}+\unicode[STIX]{x1D719}_{\infty })-{\boldsymbol{U}_{\infty }}_{(0)}\quad \text{at}~r=1,\\ \boldsymbol{v}_{(0)}=Ha\unicode[STIX]{x1D701}_{w}\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\quad \text{at walls},\\ \boldsymbol{v}_{(0)}\rightarrow \mathbf{0}\quad \text{as}~r\rightarrow \infty ,\end{array}\right\}\end{eqnarray}$$

and at $O(De)$ as

(3.29)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}_{(1)}=0,\\ \unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{(1)}-\unicode[STIX]{x1D735}p_{(1)}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D668}_{(0)},\\ \boldsymbol{v}_{(1)}={\boldsymbol{U}_{s}}_{(1)}+{\unicode[STIX]{x1D734}_{s}}_{(1)}\times \boldsymbol{r}\quad \text{at}~r=1,\\ \boldsymbol{v}_{(1)}=\mathbf{0}\quad \text{at walls},\\ \boldsymbol{v}_{(1)}\rightarrow \mathbf{0}\quad \text{as}~r\rightarrow \infty .\end{array}\right\}\end{eqnarray}$$

Here, $\unicode[STIX]{x1D668}_{(0)}$ is the disturbance polymeric stress of $O(De^{0})$ field. The angular ($\unicode[STIX]{x1D734}_{s}$) and translational velocity ($\boldsymbol{U}_{\boldsymbol{s}}$) are evaluated by using (3.24), at each order.

3.4.2 Reciprocal theorem

The symmetry of Stokes flow implies that the $O(1)$ field cannot produce a lateral lift. Therefore, the lift must arise from the $O(De)$ field. The reciprocal theorem (Ho & Leal Reference Ho and Leal1976) allows us to determine the migration associated with $\boldsymbol{v}_{(1)}$, without having to solve for the $O(De)$ field. The reciprocal theorem relates the properties of an unknown Stokes flow to a known test flow field, provided both flow fields correspond to the same geometry. The test field ($\boldsymbol{u}^{t},p^{t}$) is chosen to be that generated by a moving sphere in the positive $z$-direction (towards the upper wall) with unit velocity in a quiescent Newtonian fluid (refer to appendix A for the details of test field). Following Ho & Leal (Reference Ho and Leal1976, p. 789), we obtain the lateral migration velocity at $O(De)$ as

(3.30)$$\begin{eqnarray}De(U_{s(1)})_{z}=U_{mig}^{H}=-\frac{1}{6\unicode[STIX]{x03C0}(1+O(\unicode[STIX]{x1D705}))}De\int _{V_{f}}\unicode[STIX]{x1D668}_{(0)}\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{u}^{t}\,\text{d}V.\end{eqnarray}$$

Here, superscript $H$ denotes the hydrodynamic contribution and $6\unicode[STIX]{x03C0}(1+O(\unicode[STIX]{x1D705}))$ is the viscous resistance including wall correction at the first order (the details of derivation are shown in the supplementary material).

3.4.3 Evaluation of $\pmb{v}_{(0)}$

Following the method used for the evaluation of potential ($\unicode[STIX]{x1D713}$), we use the method of reflections to find the velocity field which satisfies the boundary conditions at the particle and wall surface with sufficient accuracy (Brenner Reference Brenner1962). As in § 3.3, we restrict the analysis to small particles i.e. $\unicode[STIX]{x1D705}\ll 1$. This condition is in addition to $De\ll 1$, which was used to expand the field variables in (3.27). We seek the disturbance ($\boldsymbol{v}_{(0)},p_{(0)}$) as

(3.31a,b)$$\begin{eqnarray}\boldsymbol{v}_{(0)}=_{(1)}\boldsymbol{v}_{(0)}+_{(2)}\boldsymbol{v}_{(0)}+\cdots ,\quad p_{(0)}=_{(1)}p_{(0)}+_{(2)}p_{(0)}+\cdots .\end{eqnarray}$$

Here $_{(i)}\boldsymbol{v}_{(0)}$ represents the $i\text{th}$ reflection of the creeping flow disturbance velocity. Since $\unicode[STIX]{x1D705}\ll 1$, we take the first two reflections into account. We substitute the above expansion in the equations governing the creeping flow hydrodynamics (3.28) and obtain the following set of problems:

(3.32)$$\begin{eqnarray}\left.\begin{array}{@{}rcl@{}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }_{(1)}\boldsymbol{v}_{(0)}\ & =\ & 0,\quad {\unicode[STIX]{x1D6FB}^{2}}_{(1)}\boldsymbol{v}_{(0)}-\unicode[STIX]{x1D735}_{(1)}p_{(0)}=\mathbf{0},\\ _{(1)}\boldsymbol{v}_{(0)}\ & =\ & U_{sx(0)}\boldsymbol{e}_{x}+\unicode[STIX]{x1D6FA}_{sy(0)}\boldsymbol{e}_{y}\times \boldsymbol{r}+Ha\unicode[STIX]{x1D701}_{p}(\unicode[STIX]{x1D735}_{(1)}\unicode[STIX]{x1D713}+\unicode[STIX]{x1D735}_{(2)}\unicode[STIX]{x1D713})\\ \ & \ & -(\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FD}z+\unicode[STIX]{x1D6FE}z^{2}+Ha(\unicode[STIX]{x1D701}_{p}-\unicode[STIX]{x1D701}_{w}))\boldsymbol{e}_{x}\quad \text{at}~r=1,\\ _{(1)}\boldsymbol{v}_{(0)}\ & \rightarrow \ & \mathbf{0}\quad \text{at}~r\rightarrow \infty ;\end{array}\right\}\end{eqnarray}$$

(3.33)$$\begin{eqnarray}\left.\begin{array}{@{}rcl@{}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }_{(2)}\boldsymbol{v}_{(0)}\ & =\ & 0,\quad {\unicode[STIX]{x1D6FB}^{2}}_{(2)}\boldsymbol{v}_{(0)}-\unicode[STIX]{x1D735}_{(2)}p_{(0)}=\mathbf{0},\\ _{(2)}\boldsymbol{v}_{(0)}\ & =\ & Ha\unicode[STIX]{x1D701}_{w}(\unicode[STIX]{x1D735}_{(1)}\unicode[STIX]{x1D713}+\unicode[STIX]{x1D735}_{(2)}\unicode[STIX]{x1D713})-_{(1)}\boldsymbol{v}_{(0)}\quad \text{at the walls},\\ _{(2)}\boldsymbol{v}_{(0)}\ & \rightarrow \ & \mathbf{0}\quad \text{at}~r\rightarrow \infty .\end{array}\right\}\end{eqnarray}$$

The first reflection ($_{(1)}\boldsymbol{v}_{(0)},_{(1)}p_{(0)}$) is found by employing Lamb’s general solution (Lamb Reference Lamb1975; Happel & Brenner Reference Happel and Brenner2012)

(3.34)$$\begin{eqnarray}\displaystyle \text{}_{(1)}\boldsymbol{v}_{(0)} & = & \displaystyle A_{1}\left(\boldsymbol{e}_{x}+\frac{x\boldsymbol{r}}{r^{2}}\right)\frac{1}{r}+B_{1}\left(-\boldsymbol{e}_{x}+\frac{3x\boldsymbol{r}}{r^{2}}\right)\frac{1}{r^{3}}+C_{1}\left(\frac{z\boldsymbol{e}_{x}}{r^{3}}-\frac{x\boldsymbol{e}_{z}}{r^{3}}\right)+D_{1}\frac{zx\boldsymbol{r}}{r^{5}}\nonumber\\ \displaystyle & & \displaystyle +\,E_{1}\left(z\boldsymbol{e}_{x}+x\boldsymbol{e}_{z}-\frac{5xz\boldsymbol{r}}{r^{2}}\right)\frac{1}{r^{5}}+F_{1}\left(\boldsymbol{e}_{x}-\frac{2z^{2}\boldsymbol{e}_{x}+x\boldsymbol{r}}{r^{2}}+\frac{2xz\boldsymbol{e}_{z}}{r^{2}}\right)\frac{1}{r^{3}}\nonumber\\ \displaystyle & & \displaystyle +\,G_{1}\left(\boldsymbol{e}_{x}-\frac{5z^{2}\boldsymbol{e}_{x}+10xz\boldsymbol{e}_{z}+13x\boldsymbol{r}}{r^{2}}+\frac{75xz^{2}\boldsymbol{r}}{r^{4}}\right)\frac{1}{r^{3}}\nonumber\\ \displaystyle & & \displaystyle +\,H_{1}\left(\boldsymbol{e}_{x}-\frac{5z^{2}\boldsymbol{e}_{x}+10xz\boldsymbol{e}_{z}+5x\boldsymbol{r}}{r^{2}}+\frac{35xz^{2}\boldsymbol{r}}{r^{4}}\right)\frac{1}{r^{5}}.\end{eqnarray}$$

(3.35)$$\begin{eqnarray}\displaystyle \text{}_{(1)}p_{(0)}=A_{1}\frac{3x}{2r^{3}}+D_{1}\frac{2xz}{r^{5}}+G_{1}\left(\frac{-30x}{r^{5}}+\frac{150xz^{2}}{r^{7}}\right). & & \displaystyle\end{eqnarray}$$

Here, the coefficients are defined as

(3.36)$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle A_{1}=\frac{3}{4}\left(U_{sx(0)}-\unicode[STIX]{x1D6FC}-\frac{\unicode[STIX]{x1D6FE}}{3}-Ha(\unicode[STIX]{x1D701}_{p}-\unicode[STIX]{x1D701}_{w})\right),\\ \displaystyle B_{1}=-\frac{1}{4}\left(U_{sx(0)}-\unicode[STIX]{x1D6FC}-\frac{3\unicode[STIX]{x1D6FE}}{5}-Ha(\unicode[STIX]{x1D701}_{p}-\unicode[STIX]{x1D701}_{w})\right)+\frac{Ha\unicode[STIX]{x1D701}_{p}}{2},\\ \displaystyle C_{1}=\unicode[STIX]{x1D6FA}_{sy(0)}-\frac{\unicode[STIX]{x1D6FD}}{2},\quad D_{1}=-\frac{5\unicode[STIX]{x1D6FD}}{2},\quad E_{1}=-\frac{\unicode[STIX]{x1D6FD}}{2},\\ \displaystyle F_{1}=\frac{\unicode[STIX]{x1D6FE}}{3},\quad G_{1}=-\frac{7\unicode[STIX]{x1D6FE}}{120},\quad H_{1}=\frac{\unicode[STIX]{x1D6FE}}{8}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The terms multiplying the coefficients $A_{1}$, $B_{1}$, $C_{1}$, $D_{1}$, $E_{1}$ represent the force-monopole (stokeslet), source-dipole, antisymmetric dipole (rotlet), symmetric dipole (stresslet) and octupole singularities, respectively (Guazzelli & Morris Reference Guazzelli and Morris2011; Kim & Karrila Reference Kim and Karrila2013). The other disturbances (terms multiplying $F_{1},G_{1},H_{1}$) are further singularities in the multipole expansion which arise due to the curvature ($\unicode[STIX]{x1D6FE}$) in the background flow field.

The evaluation of higher-order reflections can be found using Faxén transformation (Faxén Reference Faxén1922) and is identical to our previous analysis (Choudhary et al. Reference Choudhary, Renganathan and Pushpavanam2019); the details of evaluation are provided in the supplementary material. We shall later show (in § 3.4.4) that the higher reflections of the velocity field (i.e. $_{(2)}\boldsymbol{v}_{(0)}$, $_{(3)}\boldsymbol{v}_{(0)},\ldots$) are dispensable for the current problem.

The unknown translational (${U_{sx}}_{(0)}$) and rotational (${\unicode[STIX]{x1D6FA}_{sy}}_{(0)}$) velocities can be evaluated using the force-free and torque-free conditions (3.24) for first three reflections (since the second reflection does not contribute to the hydrodynamic drag and torque (Happel & Brenner Reference Happel and Brenner2012; Kim & Karrila Reference Kim and Karrila2013), we evaluate the third reflection to estimate the wall corrections)

(3.37a)$$\begin{eqnarray}\displaystyle & \displaystyle U_{sx(0)}=\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FE}/3+Ha(\unicode[STIX]{x1D701}_{p}-\unicode[STIX]{x1D701}_{w})+O(\unicode[STIX]{x1D705}^{2}), & \displaystyle\end{eqnarray}$$

(3.37b)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{sy(0)}=\unicode[STIX]{x1D6FD}/2+O(\unicode[STIX]{x1D705}^{3}). & \displaystyle\end{eqnarray}$$

The $O(\unicode[STIX]{x1D705}^{2})$ and $O(\unicode[STIX]{x1D705}^{3})$ corrections are contributions from the wall reflections of the disturbance velocity, which arise when the third reflection is taken into account. The details of the derivation of (3.37) can be found in the supplementary material. Substitution of (3.37) in (3.36) yields $A_{1}\sim O(\unicode[STIX]{x1D705}^{2})$, $B_{1}\sim Ha\unicode[STIX]{x1D701}_{p}O(1)$ and $C_{1}\sim O(\unicode[STIX]{x1D705}^{3})$. Since $\unicode[STIX]{x1D705}\ll 1$, we neglect the $O(\unicode[STIX]{x1D705}^{2})$ and $O(\unicode[STIX]{x1D705}^{3})$ wall corrections, and thus use

(3.38a-c)$$\begin{eqnarray}A_{1}=C_{1}=0,\quad U_{sx(0)}=\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FE}/3+Ha(\unicode[STIX]{x1D701}_{p}-\unicode[STIX]{x1D701}_{w}),\quad \unicode[STIX]{x1D6FA}_{sy(0)}=\unicode[STIX]{x1D6FD}/2.\end{eqnarray}$$

This yields the background flow velocity in the frame of reference of $\boldsymbol{U}_{s(0)}$ (3.13),

(3.39)$$\begin{eqnarray}\boldsymbol{U}_{\infty (0)}=\left(\unicode[STIX]{x1D6FD}z+\unicode[STIX]{x1D6FE}z^{2}-Ha\unicode[STIX]{x1D701}_{p}-\frac{\unicode[STIX]{x1D6FE}}{3}\right)\boldsymbol{e}_{x}.\end{eqnarray}$$

The above estimates of $A_{1}$, $C_{1}$ and $U_{sx(0)}$ in (3.38) when substituted in (3.34) provide an insight into leading-order velocity disturbances: (i) stresslet (${\sim}O(1/r^{2})$ arising from the resistance to strain); (ii) source-dipole field (${\sim}O(1/r^{3})$ arising from the electrophoretic slip); and (iii) octupole field (${\sim}O(1/r^{4})$ also arising from the resistance to strain). The other disturbances (multiplying $F_{1},G_{1},H_{1}$), arising due to curvature in the background flow are of $O(\unicode[STIX]{x1D705})$ and are accounted for in the evaluation of (3.30).

3.4.4 Evaluation of volume integral

To estimate the migration velocity using (3.30), we substitute the velocity field $\boldsymbol{v}_{(0)}$ into the integrand. We analyse the integral for computational convenience and divide the fluid domain ($V_{f}$) into two asymptotic subdomains (Ho & Leal Reference Ho and Leal1974, Reference Ho and Leal1976)

(3.40a,b)$$\begin{eqnarray}V_{1}=\{\boldsymbol{r}\boldsymbol{ : }1\leqslant |r|\leqslant \unicode[STIX]{x1D70C}\},\quad V_{2}=\{\tilde{\boldsymbol{r}}\boldsymbol{ : }\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\leqslant |\tilde{r}|<\infty ,-s\leqslant \tilde{z}\leqslant 1-s\}.\end{eqnarray}$$

Here, the intermediate radius $\unicode[STIX]{x1D70C}$ satisfies $1\ll \unicode[STIX]{x1D70C}\ll 1/\unicode[STIX]{x1D705}$, $\tilde{r}$ represents the outer coordinates ($\tilde{r}=r\unicode[STIX]{x1D705}$). The characteristic length scale corresponding to the inner domain ($V_{1}$) is particle radius $a^{\ast }$. Whereas, channel width $l^{\ast }$ is the characteristic length corresponding to the outer domain ($V_{2}$). Below, we inspect the nature of the integrand in each subdomain and show that the integration in $V_{2}$ is not required to attain the leading-order migration velocity.

Integration in subdomain $V_{1}$: substituting (3.34) and (3.39) in (3.19) and performing the integration (3.30) in $V_{1}$ (which is essentially a spherical shell, provided the particle is not too close to the wall, i.e. $s\gg \unicode[STIX]{x1D705}$, where $\unicode[STIX]{x1D705}\ll 1$) (the analytical integration of the volume integral is performed using Mathematica 12; the code is provided in the supplementary material) we obtain

(3.41)$$\begin{eqnarray}\displaystyle \frac{6\unicode[STIX]{x03C0}(1+O(\unicode[STIX]{x1D705}))U_{mig}^{H}}{De} & = & \displaystyle -\int _{V_{1}:r=1}^{r=\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D668}_{(0)}\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{u}^{t}\,\text{d}V=-\int _{V_{1}:r=1}^{r=\infty }\unicode[STIX]{x1D668}_{(0)}\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{u}^{t}\,\text{d}V+o(\unicode[STIX]{x1D705})\nonumber\\ \displaystyle & {\approx} & \displaystyle \frac{10\unicode[STIX]{x03C0}}{3}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}(1+3\unicode[STIX]{x1D6FF})-\frac{3\unicode[STIX]{x03C0}}{2}\unicode[STIX]{x1D6FD}Ha\unicode[STIX]{x1D701}_{p}(1+\unicode[STIX]{x1D6FF}).\end{eqnarray}$$

The details of $o(\unicode[STIX]{x1D705})$ error incurred while approximating the upper limit of the integral to $r=\infty$ is detailed in appendix B.

In view of the definition of $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FE}$, the first term in (3.41) is $O(\unicode[STIX]{x1D705})$, and is identical to that obtained by Ho & Leal (Reference Ho and Leal1976), who reported the viscoelastic migration of a particle in the absence of electrokinetic effects. The second term is $O(1)$ and is the dominant contribution to the migration in the current study. It is interesting to note that the first term arises purely from the interaction of disturbance field and the bulk flow (i.e. $\unicode[STIX]{x1D66C}_{(0)}^{(1)}$ and $\unicode[STIX]{x1D66C}_{(0)}^{(2)}$). However, the second term arises from the interactions among various disturbances (i.e. $\unicode[STIX]{x1D65A}_{(0)}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D65A}_{(0)}^{(1)}$ and $\unicode[STIX]{x1D65A}_{(0)}^{(2)}$). This insight and the proportionality to $\unicode[STIX]{x1D6FD}Ha\unicode[STIX]{x1D701}_{p}$ (second term) suggests that the $O(1)$ modification to migration arises from the interaction of stresslet and octupole disturbances (both proportional to $\unicode[STIX]{x1D6FD}$) with the electrophoresis induced source dipole disturbance (proportional to $Ha\unicode[STIX]{x1D701}_{p}$).

Integration in subdomain$V_{2}$: we now find the contribution of the integral in the outer subdomain $V_{2}$. Since the characteristic length in the outer domain is $l^{\ast }$, the disturbances are stretched into the outer coordinates ($\tilde{r}=r\unicode[STIX]{x1D705}$). In the outer coordinate representation, the order of magnitude of the leading-order disturbance and test field velocity is

(3.42a)$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\boldsymbol{v}}_{(0)}\sim Ha\unicode[STIX]{x1D701}_{p}O(\unicode[STIX]{x1D705}^{3})+\unicode[STIX]{x1D6FD}O(\unicode[STIX]{x1D705}^{2})+\unicode[STIX]{x1D6FD}O(\unicode[STIX]{x1D705}^{4})+\unicode[STIX]{x1D6FE}O(\unicode[STIX]{x1D705}^{3})+\unicode[STIX]{x1D6FE}O(\unicode[STIX]{x1D705}^{5}), & \displaystyle\end{eqnarray}$$

(3.42b)$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\boldsymbol{u}}^{t}\sim O(\unicode[STIX]{x1D705})+O(\unicode[STIX]{x1D705}^{3}). & \displaystyle\end{eqnarray}$$

In the outer coordinates, the undisturbed background flow velocity is

(3.43)$$\begin{eqnarray}\tilde{\boldsymbol{U}}_{\infty (0)}\sim \unicode[STIX]{x1D6FD}O(1/\unicode[STIX]{x1D705})+\unicode[STIX]{x1D6FE}O(1/\unicode[STIX]{x1D705}^{2})+Ha\unicode[STIX]{x1D701}_{p}O(1)+\unicode[STIX]{x1D6FE}O(1).\end{eqnarray}$$

We substitute (3.42a) and (3.43) in (3.19) and obtain $\tilde{\unicode[STIX]{x1D668}}_{(0)}$.

The volume integral (3.30), in $V_{2}$ subdomain, is represented in outer coordinates as

(3.44)$$\begin{eqnarray}\frac{6\unicode[STIX]{x03C0}(1+O(\unicode[STIX]{x1D705}))U_{mig}^{H}}{De}=-\int _{V_{2}}\tilde{\unicode[STIX]{x1D668}}_{(0)}\boldsymbol{ : }(\unicode[STIX]{x1D705}\tilde{\unicode[STIX]{x1D735}}\tilde{\boldsymbol{u}}^{t})\unicode[STIX]{x1D705}^{-3}\,\text{d}\tilde{V}.\end{eqnarray}$$

Upon simplification, we find that the integrand is $O(\unicode[STIX]{x1D705}^{2})$. Since $\unicode[STIX]{x1D705}\ll 1$, the dominant contribution can be concluded to arise from the inner integral (i.e. $O(1)$) and the contribution from the outer integral can be neglected. In § D.1, the expression for leading-order migration velocity (3.41) is verified using an independent calculation following the theorem of Giesekus (Reference Giesekus1963) for the special case of polymeric fluids ($\unicode[STIX]{x1D6F9}_{2}^{\ast }=-0.5\unicode[STIX]{x1D6F9}_{1}^{\ast }$).

4 Results and discussion

The expression for the migration velocity (3.41) shows that the direction of lateral migration is determined by the magnitude and sign of the viscometric parameter $\unicode[STIX]{x1D6FF}$ and $De$. For most viscoelastic fluids $De$ is positive, and $\unicode[STIX]{x1D6FF}$ is reported to be between $-0.5$ and $-0.7$ (Caswell & Schwarz Reference Caswell and Schwarz1962; Leal Reference Leal1975; Bird et al. Reference Bird, Armstrong and Hassager1987; Koch & Subramanian Reference Koch and Subramanian2006). Since $\unicode[STIX]{x1D6FE}$ is negative, the first term in (3.41) is strictly positive below the centreline ($\unicode[STIX]{x1D6FD}>0$) and negative above it ($\unicode[STIX]{x1D6FD}<0$), i.e. the pure elastic stresses push the particle towards the channel axis at $O(\unicode[STIX]{x1D705})$. The second term is $O(1)$; its proportionality to $Ha\unicode[STIX]{x1D701}_{p}$ indicates that this dominant contribution originates due to the electrophoretic slip, i.e. incorporation of electrophoresis enhances the migration by $O(\unicode[STIX]{x1D705}^{-1})$. For $-0.7<\unicode[STIX]{x1D6FF}<-0.5$, a particle with positive $\unicode[STIX]{x1D701}_{p}$ will migrate towards the walls (at $O(1)$) when the electric field is applied in the direction of the flow. The reverse would hold for a particle with negative $\unicode[STIX]{x1D701}_{p}$ or if the electric field is applied in the opposite direction. These findings qualitatively agree with the experimental observations depicted in figure 1.

The enhancement demonstrates that an addition of electrokinetic slip on the particle surface can drastically alter the normal stresses which consequently enhances the particle migration. The slip creates a source-dipole disturbance in the fluid which decays rapidly away from the particle (${\sim}O(1/r^{3})$). This suggests that the migration occurs from short-ranged interactions, which are unaffected by the boundaries unless the particle is very close to the walls. Furthermore, equation (3.41) reveals that the parameters such as wall zeta potential ($\unicode[STIX]{x1D701}_{w}$) and curvature of the background flow ($\unicode[STIX]{x1D6FE}$) do not affect the migration at the present order of approximation.

Figure 3.

Particle stream width, PSW, ($\unicode[STIX]{x03BC}\text{m}$) of negative zeta potential particles suspended in a Poiseuille flow of second-order fluid. Parameters: $De=0.12$; $\unicode[STIX]{x1D6FF}=-0.5$; $E_{\infty }^{\ast }=300~\text{V}~\text{cm}^{-1}$; $\unicode[STIX]{x1D701}_{p}=-2.52$ ($\unicode[STIX]{x1D701}_{p}^{\ast }=-83~\text{mV}$); $l^{\ast }=66~\unicode[STIX]{x03BC}\text{m}$; $\unicode[STIX]{x1D705}=0.017$.

Since the experimental parameters are in the regime of the theoretical analysis ($De<1$, $De\gg Re_{p}$, $Ha\unicode[STIX]{x1D701}_{p}\sim O(1)$, $\unicode[STIX]{x1D705}\ll 1$), we now quantitatively compare the two. Theoretically, the trajectories of a single particle are calculated using (3.41) and are converted to an equivalent particle stream width using $|l^{\ast }-2(l^{\ast }-d^{\ast })|$, for particle focusing at the centreline ($l^{\ast }$ being the channel height and $d^{\ast }$ is the particle to wall distance). (In § D.2, the theoretical calculation of trajectory is detailed and validated with that reported by Ho & Leal (Reference Ho and Leal1976) for $Ha=0$.) We use the relaxation time scale ($(\unicode[STIX]{x1D6F9}_{1}^{\ast }+\unicode[STIX]{x1D6F9}_{2}^{\ast })/\unicode[STIX]{x1D707}^{\ast }$) corresponding to $C_{PEO}^{\ast }=250~\text{ppm}$ as 2.75 ms (Lu & Xuan Reference Lu and Xuan2015). Figure 3 shows the comparison of particle stream widths for two cases: (i) no electric field (i.e. viscoelastic focusing) and (ii) parallel electric field (i.e. electroviscoelastic focusing). For case (i), the theoretical prediction shows modest particle migration. The experiments show negligible focusing as some particles get trapped near the channel corners, which is one of the limitations of viscoelastic focusing at low flow rates (Lu et al. Reference Lu, Liu, Hu and Xuan2017; Tian et al. Reference Tian, Feng, Chen, Liu, Li and Sun2019). Since our model is valid for two-dimensional flows, it does not account for the trapping of particles at three-dimensional corners. For case (ii), the theoretical predictions and experimental results confirm an order of magnitude enhancement in focusing. The trapping of particles at the corners is not seen here as the repulsive Maxwell lift is strongest near the corners (Liang, Qian & Xuan Reference Liang, Qian and Xuan2010). Thus, the agreement between the experiments and theoretical analysis is improved.

Figure 4.

Comparison of various components of the migration velocity ($\text{m}~\text{s}^{-1}$) for $\unicode[STIX]{x1D705}=0.017$; Maxwell migration, MM; viscoelastic migration, VM; and electroviscoelastic migration, EVM. Here EVM corresponds to figure 1(b). The inset shows the comparison of MM and VM. Parameters: $De=0.12$; $\unicode[STIX]{x1D6FF}=-0.5$; $E_{\infty }^{\ast }=300~\text{V}~\text{cm}^{-1}$; $\unicode[STIX]{x1D701}_{p}=-2.52$ ($\unicode[STIX]{x1D701}_{p}^{\ast }=-83~\text{mV}$); $l^{\ast }=66~\unicode[STIX]{x03BC}\text{m}$; $\unicode[STIX]{x1D705}=0.017$.

Our mathematical framework reveals that the migration in figure 1(b,c) arises from: (i) Maxwell stresses (3.26); (ii) pure viscoelastic stresses (3.41-term 1); and (iii) electrophoresis enhanced viscoelastic stresses (3.41-term 2). A comparison of these migration velocities, namely, Maxwell migration (MM), viscoelastic migration (VM) and electroviscoelastic migration (EVM) is shown in figure 4. The curves correspond to the case in figure 1(b): a particle with negative zeta potential and electric field applied parallel to the flow. The migration due to decoupled Maxwell stresses (i.e. MM) is dominant near the walls and decays rapidly away from it. Although this wall-repulsive force is negligible in the bulk of the channel (for $\unicode[STIX]{x1D705}\ll 1$), it helps eliminate the corner equilibrium positions (Liang et al. Reference Liang, Qian and Xuan2010). The other components VM and EVM act towards the centre as depicted by the negative gradient at the central equilibrium position. For conditions corresponding to figure 1(c), where flow and electric field are opposite to each other ($Ha<0$), the EVM curve would reverse, as we can see from (3.41). As a consequence, the central equilibrium position would become unstable, resulting in particle migration towards the walls.