Nonlinear electro-osmosis in dilute non-adsorbing polymer solutions with low ionic strength

Nonlinear behavior of electro-osmosis in dilute non-adsorbing polymer solutions with low salinity is investigated with Brownian dynamics simulations and a kinetic theory. In the Brownian simulations, the hydrodynamic interaction between the polymers and a no-slip wall is considered with Rotne-Prager approximation of Blake tensor. In a plug flow under a sufficiently strong applied electric field, the polymer migrates toward the bulk, forming a depletion layer thicker than the equilibrium one. Consequently, the electro-osmotic mobility increases nonlinearly with the electric field and gets saturated. This nonlinear mobility qualitatively does not depend on the details of rheological properties of the polymer solution. Analytical calculation of the kinetic theory for the same system reproduces quantitatively well the results of the Brownian dynamics simulation.


Introduction
Electro-osmosis is observed widely in many systems such as colloids, porous materials, and biomembranes. It characterises the properties of interfaces between solids and electrolyte solutions. 1,2 Interest in the applications of electro-osmosis has been growing recently. For instance, it is used to pump fluids in microfluidic devices, as it is more easily implemented than pressure-driven flow. 3 Its application to electrical power conversion in chemical engineering is also very fascinating. 4,5 When the electrokinetic properties of a surface are characterised by a zeta potential, the Smoluchowski equation is often employed in conjunction with measurements of the electro-osmotic or electrophoretic mobilities. However, the validity of this equation should be considered more seriously. It is derived from the Poisson-Boltzmann equation and Newton's constitutive equation for viscous fluids. The zeta potential is defined as the electrostatic potential at the plane where a no-slip boundary condition is assumed. When these equations cannot be validated, the Smoluchowski equation is also questionable. The Poisson-Boltzmann equation, the simple hydrodynamic equations, or both sometimes do not work well for a strong coupling double layer, 6 inhomogeneity of viscosity and dielectric constant near the interface, 7,8 and non-Newtonian fluids, [9][10][11][12] for example.
To control the electrokinetic properties of charged capillaries, the structures of liquid interfaces in contact with charged surfaces are modified by grafting or adding polymers. 13 In capillary electrophoresis, for example, the electro-osmotic flow is reduced by grafted polymers on the interfaces. Several studies of surfaces with end-grafted charged and uncharged polymers have also been reported. [14][15][16][17][18][19] Under a weak applied electric field, the grafted polymer remains in the equilibrium configuration, and the resultant electro-osmotic velocity behaves linearly with respect to the electric field. To measure the mobility of such a surface, the hydrodynamic screening and anomalous charge distributions due to the grafted polymers are important. [14][15][16][17][18][19] When a sufficiently strong electric field is applied, the polymers are deformed by the flow and the electric field; thus, the electroosmotic velocity becomes nonlinear. 16 Note that the end-grafted polymers cannot migrate toward the bulk, as one of the ends is fixed on the surfaces.
When we add polymers to solutions, a depletion or adsorption layer is often formed near a solid wall, as well as diffusive layers of ions in equilibrium states. The interaction between the polymers and the wall determines whether the polymers are depleted from or adhere to the surfaces. The thickness of the depletion or adsorption layer is of the same order as the gyration length of the polymers. When polymers adhere to the wall, the viscosity near the wall becomes large, so the electro-osmotic mobility is strongly suppressed. 20 Moreover, it is known that an adsorption layer of charged polymers can change the sign of mobility. 7,[21][22][23][24] The curvature of the surface also modulates the surface charge density and even increases the mobility beyond the suppression caused by viscosity enhancement. 20 Electro-osmosis of a non-adsorbing polymer solution has been analysed in terms of two length scales: the equilibrium depletion length d 0 and the Debye length l. 7,25 In the depletion layer, the viscosity is estimated to be approximately the same as that of the pure solvent, and it is smaller than the solution viscosity in the bulk. When the Debye length is smaller than the depletion length, the electro-osmotic mobility is larger than that estimated from the bulk value of viscosity. Typically, for 10 mM electrolyte solutions, one has l E 3 nm and d 0 E 100 nm. In such a case, [26][27][28] an electro-osmotic flow with a high shear rate is localised at a distance l from the wall. Thus, the electro-osmotic flow profile and resultant electro-osmotic mobility are almost independent of the polymers. Such behaviours are experimentally observed in solutions of carboxymethyl cellulose with urea. 27 On the other hand, in the solutions of small polymers with low salinity (typically for 0.1 mM electrolyte solutions, l E 30 nm and d 0 E 5 nm), the electro-osmotic mobility is suppressed by the polymeric stress. 25 When a sufficiently strong electric field is applied, the electroosmosis of a polymer solution shows nonlinear behaviours. 27,29 These nonlinearities are theoretically analysed by models of uniform non-Newtonian shear-thinning fluids. [9][10][11][12] Assuming that polymers remain in the interfacial layers and viscosity depends on the local shear rate, as in power-law fluids, their phenomenological parameters differ from those in the bulk, as the concentration in the interfacial layers differs from the bulk concentration. 27 Thus, the understanding of nonlinear electroosmosis remains phenomenological. Furthermore, when shear flow is applied to polymer solutions near a wall, it is experimentally and theoretically confirmed that cross-stream migration is induced toward the bulk. [30][31][32] The concentration profiles of the polymer near the wall have been calculated, and the depletion length dynamically grows tenfold larger than the gyration radius. 31 However, these hydrodynamic effects in the electrokinetics have not been studied to date, to the best of our knowledge.
In this context, this paper discusses another origin of nonlinearity, which is induced by hydrodynamic interaction between the polymer and the wall, considering mainly situations where d 0 { l. For this purpose, this paper is organised as follows. Section 2 presents a toy model of nonlinear electro-osmosis of dilute polymer solutions. Section 3 describes the Brownian dynamics simulation, and Section 4 presents the results of the simulation. In Section 5, we discuss an analytical approach to nonlinear electroosmosis by using the kinetic theory of cross-stream migration. 31 Section 6 outlines the main conclusions.

Toy model
First, we propose a toy model for electro-osmosis of polymer solutions. A dilute solution of non-adsorbing polymers is considered. The viscosity of the solution is given by where Z 0 is the viscosity of the pure solvent, and Z sp is the specific viscosity of the solution. The gyration length of the polymers is defined as d 0 , which is of the same order as the equilibrium depletion length. It is assumed that the polymers have d 0 E 100 nm. Ions are also dissolved in the solution with the Debye length l. When well-deionised water is considered, the Debye length is on the order of l E 10 3 nm, although such salt-free water is rarely realised owing to the spontaneous dissolution of carbon dioxides. The interfacial structure near a charged surface is characterised by l and d 0 . When an external electric field is applied, a shear flow is imposed locally within a distance l from the wall, and the resultant shear rate is where m 0 is the electro-osmotic mobility of the pure solvent and is typically estimated as m 0 E 10 À8 m 2 (V s) À1 . According to studies of cross-stream migration in uniform shear flow, 31 the depletion layer thickness depends on the shear rate, Fig. 1(a) schematically shows the thickness of the depletion layer as a function of the electric field strength. Fig. 1(b) shows the nonlinear electro-osmotic mobility. The mobility increases and is saturated with increasing electric field. The threshold electric field E 0 is typically 10 3 V m À1 , which is experimentally accessible.

Model for simulation
In this section, our method of Brownian dynamics simulation is described. As shown in Fig. 2(a), a dumbbell is simulated in an electrolyte solution with a no-slip boundary at z = 0. The dumbbell behaves like a dilute solution. The solvent is described as a continuum fluid with the viscosity Z 0 and fills the upper half of the space (z 4 0). Electrolytes are also treated implicitly with the Debye length l = k À1 . The dumbbell has two beads whose hydrodynamic radii are a, and each bead consists of many monomeric units of the polymer [see Then we solve the overdamped Langevin equations 33 given by where X na is the a component of vector X n . Furthermore, u 0 (z) is the external plug flow, r na = q/qx na , G is the mobility tensor, F n = Àr n U is the force exerted on the nth bead, and U is the interaction energy given as a function of x n . n n is the thermal noise which satisfies the fluctuation-dissipation relation as To include the effects of the no-slip boundary, the Rotne-Prager approximation for the Blake tensor 34 is used as the mobility tensor for distinct particles (n a m), 35,36 although it is valid only for particles separated by a large distance. In this study, we neglect lubrication corrections for nearby particles. 37 The Blake tensor for the velocity at x 2 induced by a point force at x 1 with the no-slip boundary at z = 0 is given by the Oseen tensor and the coupling fluid-wall tensor as 34 where q = x 2 À x 1 , R = x 2 À % x 1 , and % x 1 is the mirror image of x 1 with respect to the plane z = 0 [see Fig. 2 The Oseen tensor is given by where q is the magnitude of q. The second term in eqn (10) is where and r qg = q/qq g . The Rotne-Prager approximation of the Blake tensor is given by 6,35-38 The mobility tensor for the self-part is given by 6,35-38 where Finally we obtain the mobility tensor as The non-uniform mobility term in eqn (8) can be simplified within the Rotne-Prager approximation of the Blake tensor because the following relation holds: X b¼x;y;z Thus, the non-uniform mobility term is rewritten as X The interaction energy includes the spring and bead-wall interaction given by where U s is the spring energy, where a FENE dumbbell is a finitely extensible nonlinear elastic dumbbell, and the parameter b = Hq 0 2 /k B T is defined for convenience. U w is the bead-wall interaction, 39 which is purely repulsive: Eqn (8) is solved numerically. A reflection boundary condition is set at z = D. When the centre of the dumbbell crosses the boundary, the z coordinates of each bead are transformed from z to 2D À z. In the lateral directions, periodic boundary conditions are imposed. The size of the lateral directions is L Â L.

Results of simulation
The concentration and velocity profiles are calculated as and where d(z) is the delta function, du(z) = u(z) À u 0 (z) is the increase in velocity due to polymeric stress, and hÁ Á Ái is the statistical average in the steady state. Eqn (25) is derived in Appendix A. For a surface with a small zeta potential compared to k B T/e, where e is the elementary charge, the imposed electro-osmotic flow u 0 (z) is given by where m 0 is the electro-osmotic mobility in the pure solvent, and E is the applied electric field. 2 Eqn (8) is rewritten in a dimensionless form with the length scale d 0 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi k B T=H p and time scale t = 6pZ 0 a/4H. Different types of dumbbells are simulated using the parameters listed in Table 1.
Note that the simulated systems are treated as dilute systems, and the linearity with respect to the bulk polymer concentration is preserved. After sample averaging, we obtain the concentration at the upper boundary c(D), which deviates slightly from (L 2 D) À1 because of the inhomogeneity near the surface. Hereafter, we define the normalised concentration as The velocity increment du(z), as well as the concentration profile, is linear with respect to c(D). For convenience, we set a characteristic concentration c b = 0.1d 0

À3
, and the nonlinear electro-osmotic mobility is defined as The top boundary is placed at D = 100d 0 , the lateral size is set to L = 1000d 0 , and the Debye length is set to l = k À1 = 10d 0 . We also set w = 3k B T and define a hydrodynamic parameter h Ã as 31 Fig. 3 shows the steady-state profiles of the Hookian dumbbell concentration as functions of the distance from the wall. In the equilibrium state of E = 0, the profile shows a depletion layer whose width is of the same order as the gyration length d 0 .
When the applied electric field is increased further, the depletion layer becomes larger than the equilibrium one, and a peak is formed. The inset in Fig. 3 shows the depletion length as a function of the applied field. The depletion length is defined by the position of the concentration peak. It shows power-law behaviour with an exponent of 0.22, which is much smaller than that of 2.0 for a uniform shear flow. 31 The value of concentration at the peak also increases as the electric field increases. The results given above are for the Hookian dumbbell which is infinitely extensible with shear deformation. To consider more realistic polymers, the finitely extensible nonlinear elastic (FENE) dumbbell is simulated. Fig. 4(a) shows the concentration profiles at E = 1000E 0 . Interestingly, one-peak behaviour is also observed in the FENE dumbbells. For the Hookian dumbbell, the concentration near the surface remains finite. On the other hand, for the FENE dumbbells, the concentrations near the surface are negligibly small. Fig. 4(b) plots the electro-osmotic mobilities with respect to the applied electric field. The resultant electro-osmosis clearly increases nonlinearly with respect to the applied electric field. When the applied field becomes stronger, the mobility increases and is saturated, similar to that in the toy model. The two types of dumbbells have different rheological properties in the bulk, 40-42 so this nonlinearity is not due to the rheological properties of the dumbbells. On the other hand, the mobility is almost constant for E t 10E 0 , and this threshold of linearity is larger than E 0 , which is predicted by the toy model. Likewise, saturation is observed when E E 10 4 E 0 , which is larger than E 1 . To clarify the difference in the profiles near the surface, hq 2 i and hq z 2 /q 2 i are plotted with respect to the distance from the surface. Fig. 5(a) shows the profiles of hq z 2 /q 2 i. In the bulk, they approach 1/3, which means that the dumbbells are distributed isotropically. Near the surface, the polymers are inclined by the shear flow. The Hookian dumbbell has the largest angle between the z axis and the dumbbell direction. Fig. 5(b) plots the profiles of hq 2 i. In the bulk, they approach 3d 0 , which is their equilibrium value. Near the surface, they become larger because the polymers are elongated by the shear flow. For the FENE dumbbells, saturation of the dumbbell length is observed. These behaviours differ greatly from the minor differences in the concentration profiles.

Kinetic theory
In this section, a kinetic theory for a dumbbell is developed on the basis of the Ma-Graham theory. 31 The probability function C(x 1 ,x 2 ,t) obeys the continuity equation where : x n is the flux velocity given by 43 In the kinetic model, the beads are treated as point-like particles. Thus, the mobility tensor is obtained by using G B instead of G RPB for both the self and distinct parts. The continuity equation can be rewritten using q and r as @C @t ¼ Àr r Á ð_ rCÞ À r q Á ð _ qCÞ; where r = (x 1 + x 2 )/2 is the centre of the mass of the dumbbell. We also define r 1 and r 2 as Then, the probability function is also regarded as a function of r and q. Here we neglect the interaction between the wall and the beads. The flux velocities for r and q are obtained as   where F s = Àr 1 U s is the spring force, and D K is the Kirkwood diffusion tensor, which characterises the diffusivity of macromolecules and is given by % G and Ĝ are variants of the mobility tensors defined as The concentration field c(r,t) can be obtained by integrating the probability function with respect to the spring coordinate. It is given by We also define the probability function only for the spring coordinate asĈ ðq; t; rÞ ¼ Cðr; q; tÞ cðr; tÞ : These new fields satisfy the continuity conditions, such that where hÁ Á Ái q represents the average with the spring coordinate: For the limit of q { r, % G and D K can be expanded using r. Keeping only the leading term, we obtain where Note that the approximation is more accurate than that in a previous study 31 as that study considered only w E 1, which is not satisfied near the surface. With the approximation, eqn (36) is averaged byĈ, and finally we obtain the concentration flux for the z direction as where Eqn (48) indicates two opposite fluxes of the polymers due to the external flow field. One is the migration flux from the wall toward the bulk and originates from the hydrodynamic interaction between the wall and the force dipoles. 31 The other is the diffusion flux from the bulk to the surface wall and is not found for polymers in uniform shear flows. 31 Note that the second flux includes not only the ordinary diffusion flux hD K zz i q r rz c, but also the diffusion flux due to q inhomogeneity, cr rz hD K zz i q . When the external shear flow is uniform, the second flux vanishes, and the depletion length is proportional to the square of the shear rate, as the migration velocity is proportional to the normal stress difference. 31 In a plug flow, the diffusion flux suppresses the growth of the depletion layer, and this behaviour may explain why the exponent of the depletion length is much smaller than 2.0 in a uniform shear flow. In the steady states of electro-osmosis, the total flux in eqn (48) becomes zero; thus, This equation shows that the migration flux and the diffusion flux are balanced at the peak of the concentration profiles. Finally the concentration profile is calculated as The resultant flow profile can be calculated as where r p is the polymeric part of the stress tensor: To obtain explicit expressions for c and du, it is necessary to estimate u mig , hD K zz i q , and r p . For this purpose, eqn (43) should be analysed. However, eqn (43) is highly complicated. Even without the wall effects, it cannot be solved exactly, so several approximation methods have been proposed. 44 For simplicity, all the hydrodynamic interactions are ignored; thus, the continuity equation is given as For the Hookian dumbbell, eqn (54) can be solved for the second moment of q, and for the FENE dumbbell, a pre-averaged approximation 41,42 is employed. The curved lines in Fig. 5 are calculated using these approximations, and they exhibit good quantitative agreement with the simulation results. Appendix B gives approximate expressions for these quantities of the Hookian and FENE dumbbells. Fig. 6(a), (c), and (e) show the concentration profiles for the applied field E = 1000E 0 . The points are obtained by the Brownian dynamics simulation and the curved lines are obtained by the kinetic theory. The theoretical calculations quantitatively cover the simulations well. Moreover, they reproduce the differences in the concentration near the surface between the Hookian and FENE dumbbells, as the migration velocity can be approximately proportional to hwi q (see Appendix B), and it is greatly suppressed in the Hookian dumbbells. Fig. 6(b), (d), and (f) show the nonlinear electro-osmotic mobilities with respect to the applied field. The theoretical curved lines also exhibit acceptable agreement with the simulation results. However, they are not as consistent with the simulation results under weak applied electric fields, as the equilibrium depletion layer is not considered in the kinetic theory.

Summary and remarks
Brownian dynamics simulations are used to study nonlinear electro-osmotic behaviour of dilute polymer solutions. The simulation results agree with a toy model and analytical calculations using a kinetic theory. The main results are summarised below.
(i) Under an external plug flow, the polymer migrates toward the bulk. The concentration profile of the polymer shows a depletion layer and a single peak. The thickness of the depletion layer depends on the electric field. At the peak, the migration flux is balanced by the diffusion flux.
(ii) The growth of the depletion layer leads to an increase and saturation of the electro-osmotic mobility. This behaviour does not depend qualitatively on the rheological properties of the dumbbells.
(iii) The results of analytical calculation of the concentration and the nonlinear mobility using the kinetic theory agree with the results of the Brownian dynamics simulation. The threshold of the electric field for nonlinear growth and saturation of the mobility is much larger than that predicted by the toy model, as the diffusive flux suppresses the migration toward the bulk due to the inhomogeneous shear flow.
We conclude this study with the following remarks.
(1) Nonlinear electro-osmosis with l { d 0 has already been observed experimentally. 26,27 These studies reported that the mobility increased with increasing electric field. However, nonlinear electro-osmosis with l c d 0 has not been reported experimentally; therefore, experimental verification of our findings is highly desired.
(2) It remains a future problem to determine whether the hydrodynamic interaction between the polymers and the surface plays an important role in the electro-osmosis of polymer solutions with l { d 0 . In this case, the elongation of the polymers is strongly inhomogeneous under plug flow with a short Debye length; thus, more realistic chain models should be considered.
(3) Addition of charged polymers to solutions can change the direction of the linear electro-osmotic flow. 7,23 When a sufficiently strong electric field is applied to this system, the flow might recover its original direction. This needs to be investigated theoretically and experimentally.

A Derivation of the velocity equation for Brownian dynamics simulation
In this appendix, the derivation of eqn (25) is explained. The velocity field induced by the polymer is given by and the polymeric part of the stress tensor is obtained by averaging the microscopic expressions in the lateral directions as Here the microscopic expression of the stress tensor is given bŷ where F nm is the force exerted on the n-th bead by the m-th bead, and d s nm (x) is the symmetrised delta function given by d s nm ðxÞ ¼ The symmetrised delta function is integrated in the lateral directions as where y(z n À z) = 1 À y(z À z n ). Then we obtain Þy z À z n ð ÞÀ z À z m ð Þy z À z m ð Þ z m À z n ¼ min z; z n ð ÞÀmin z; z m ð Þ z n À z m ; where min(z,z n ) = zy(z) À (z À z n )y(z À z n ). Finally, the velocity increment is expressed as where Therefore, we have and the polymeric stress tensor is The Kirkwood diffusion constant can be estimated as where the second term is split into second-order moments; thus, we obtain It is differentiated by z as The migration velocity can be estimated using the splitting approximation of the averages as u mig ðzÞ ¼ 3k B T 64pZ 0 z 2 w q x 2 þ q y 2 À Á À 2w q % 3k B T 32pZ 0 z 2 hwi q q x 2 þ q y 2 À 2 q ¼ 3k B T 32pZ 0 z 2 hwi q f 2 ; (69) where

B.2 FENE dumbbell
The second moment of the spring coordinate for a FENE dumbbell can be obtained by pre-averaged closures of a p-FENE model. 41,42 It is given by The polymer stress tensor is The Kirkwood diffusion constant is D K