Numerical simulation is used to calculate accurately the electrophoretic mobility of a charged spherical nanoparticle confined in a nanochannel, under a weakly applied electric field. Classic models for electrophoretic mobility are valid only in the linear regime of small particle zeta potential, and for an unbounded fluid domain. However, these models fail to predict the electrophoretic mobility measured experimentally in bounded nanochannels. We adopt asymptotically expanded formulations and solve the fully coupled equations on a 3D finite element domain. Factors affecting particle mobility include electrolyte concentration, channel size, and zeta potentials on both the particle surface and channel walls. Specifically, spherical particles are examined with diameters 2a = 10 and 50 nm, in a 100 nm high channel. The non-dimensional electric double layers were varied between 0.1 < κa < 100. The results indicate that the mobility of a particle located at the nanochannel centerline agrees to within 1% of the average mobility of a particle distributed transversely throughout the nanochannel. Furthermore, confinement by the nanochannel walls was found to affect greatly nanoparticle mobility. As a result, it is feasible to use nanochannels to separate two different size nanoparticles, even when the particles have equal zeta potentials. Finally, a new method is proposed to estimate accurately particle and wall zeta potentials by contrasting the observed differences in mobility between observed in two different height channels.
I. INTRODUCTION
Electrokinetic phenomena in microchannels and nanochannels contribute to many practical applications in colloidal and biomedical sciences.1–3 For example, the electrophoretic mobility of a charged particle can be used to characterize its zeta potential,4 which is the key parameter for determining the stability of colloid dispersions. In addition, the separation of different colloidal and analyte species in a microfluidic device can be applied to directly identify biomolecules and particles for a variety of different applications, for example, disease diagnosis.5,6 A particle behaves differently in a channel compared to an unbounded domain, especially when the thickness of electric double layer (EDL) is on the order of channel size. The confinement effect not only induces additional hydrodynamic drag on the particle but also affects the electric field near the particle, thereby altering the resulting electrophoretic force. Therefore, a comprehensive model is required to capture the underlying coupled physics.
Electrokinetic phenomena in EDL have been studied extensively,7–18 including specifically the electrophoretic mobility of particles.19–28 Smoluchowksi22 and Huckel23 studied limiting cases of thin and thick EDLs on both particles and surfaces, and obtained a simple linear relation between electrophoretic mobility and particle zeta potential. Henry24 derived an expression for electrophoretic mobility of a spherical particle with a finite EDL thickness. However, these studies only focused on the linear regime, i.e., using small zeta potentials. Wiersema et al.25 calculated mobility of a particle for high zeta potentials using a Gouy-Chapman model for the EDL. O’Brien and White26 later directly solved the linearized equations for a particle of high zeta potentials from low to high electrolyte concentration. They found that the linear relation between mobility and zeta potential is only valid for low zeta potential cases. As the zeta potential increases, the electrophoretic mobility reaches a maximum and then decreases in a thin double layer system. This is because the drag force increases faster with increasing zeta potential than the driving electric force. Ohshima et al.27 derived an analytical expression for the mobility of a spherical particle in a symmetric electrolyte solution. This expression agrees well with the results of O’Brien and White. Khair and Squires28 showed the slip enhances particle mobility for thick double layers by reducing viscous drag. In the case of thin double layers, however, mobility increases from low to moderate values (0–50 mV) of zeta potential. Further increases in zeta potential leads to decreased mobility, and approach a limiting value independent of slip length, due to nonuniform surface conduction.
In terms of examining the mobility of particles in confined systems, the boundary effects for spheres located near walls have been widely studied.29–34 Keh and Anderson29 studied a non-conducting sphere near a boundary for a very thin EDL on the particle. A charged boundary generates electro-osmotic flow, which affects both electrical and velocity fields. Keh and Chen30 derived an analytical expression of electrophoretic mobility of a charged spherical particle near a charged plane wall in the thin EDL limit. When the gap between the particle and the wall is small, the mobility is enhanced up to 23%, due the squeezed electrical field lines. Ennis and Anderson31 used method of reflections to study a spherical particle near a charged wall and in a cylindrical pore. Later Shugai and Carnie32 studied similar problems with different methods. They found that the linear relation between electrophoretic mobility and zeta potential does not hold due to complex EDL interactions. Hsu et al.33 analyzed electrophoretic mobility of a spherical particle in a confined cylindrical pore with linearized equations assuming a low particle zeta potential. In the case of both charged particle and pore, they found that the magnitude of mobility has a local maximum as electrolyte concentration increases, which is due to the interaction of the double layer between the particle and the pore. Hsu and Chen34 then extended their model to include the effects of double layer polarization and electro-osmotic flow caused by the charged pore. For a positively charged particle in a positively charged pore with low surface potential, the mobility has a minimum and the direction of the motion may change twice as double layer thickness changes. To summarize, in all cases to date, electrophoresis of a particle is only considered under certain restrictions (low zeta potential, very thin or thick EDL) in an unbounded or idealized domain (spherical/cylindrical pore). However, these studies can never accurately capture nanochannel experimental data, where the zeta potential is not low, the EDL is of variable thickness (not limited to thick or thin) and the domain is a fixed rectangular system.
In the current study, numerical simulation is used to examine the mobility of a charged spherical particle of any size EDL driven by a weak electric field in a rectangular channel, where the magnitude of the applied electric field is much smaller than the electric field generated inside the EDL. We use asymptotically expanded formulations to simplify the governing equations, and solve the fully coupled equations on a 3D domain. The results are first validated in an unbounded domain with existing theory. Next, the effect of particle location in the nanochannel is examined. The results indicate that the mobility of a particle distributed transversely throughout the nanochannel agrees to within 1% of a particle located at the nanochannel centerline. Using the centerline for particle location, we investigated particle mobilities for varying zeta potentials and electrolyte concentrations, which showed a large dependence on EDL thickness and channel height. We next investigate the ability to improve separation efficacy of particles in a nanochannel. Finally, a new method is proposed to characterize accurately particle and wall zeta potentials, by comparing observed particle mobilities measured in microchannels to those measured in nanochannels.
II. THEORY
A. Governing equations
Consider a spherical particle in a channel filled with electrolyte with valences of z1 and z2 of cations and anions respectively as shown in Figure 1. The flow is assumed to be incompressible due to assumed low Reynolds number flow in nanoconfined channels. To simplify the system, we further assume that the electrolyte is symmetric and monovalent (KCl for example), that is, z1 = −z2 = 1. Thus the bulk ion concentrations and diffusion coefficients are simplified as n1∞ = n2∞ = n∞, D1 = D2 = D. Following the work of Khair and Squires,28 for a weakly applied electric field, a regular perturbation expansion can be used to simplify the governing equations. Suppose particle zeta potential is a constant value ζp, the scale of electric field in the EDL is ζp/λD, where λD = 1/κ = (ɛkBT/2n∞e2)1/2 is Debye length. Next, a small parameter δ is defined as the ratio of applied electric field to the electrostatic field, resulting from the surface charge on the nanoparticle, orδ ≡ E∞λD/ζp, where E∞ is the applied electric field. In this study we assume the following constants: E∞ ∼ 1000 V/m, λD ∼ 10−8 m, ζp ∼ 10−2 V, so δ is O(10−3).
The dimensionless governing equations can then be rewritten as electrostatic and perturbed equations as follows:
where dependent variables |$\hat \phi$|, |$\hat n_i$|, |${\bf \hat u}$|, and |$\hat p$| are the dimensionless electric potential, ion number concentration of species i, flow velocity, and pressure, respectively. The subscript “0” and “1” correspond to electrostatic and perturbed states, respectively. The dimensionless charge density is defined as |$\hat \rho = \sum {z_i \hat n_i } = \hat n_1 - \hat n_2$|. Here we use ϕc = kBT/e, the particle radius a, and the bulk concentration n∞ to scale electric potential, length, ion number concentration, respectively. The velocity is scaled by the electrophoretic velocity |$U_c = \varepsilon \phi _c^2 /\eta a$|. The pressure is scaled by pc = ηUc/a, and the electric force is scaled Fc = ρca3(ϕc/a) = n∞kBTa2. In addition, we define the Peclet number as α = 2Uca/D. The parameters ɛ, e, kB, T, and η correspond to permittivity, electron charge, Boltzmann constant, temperature, and viscosity, respectively. Proper boundary conditions for both electrostatic and perturbed states are required to solve the complete system.
B. Boundary conditions
In our system, the surfaces of the particle and the channel walls are nonconductive and impermeable with a no-slip condition. The coordinates are fixed on the particle, with the origin at the center as shown in Figure 1. Zeta potentials are specified at the surfaces of the particle and walls, and we apply constant voltages at both ends of the channel to produce a constant applied electric field along the channel. We fix the coordinates to the particle so that the channel walls are assigned the negative particle velocity. A no stress condition is used at the inlet and outlet far from the particle, where the concentration should reach the equilibrium value. Thus we have the following boundary conditions:
where n is the unit normal vector. The applied electric field is specified by E∞ = (Vin − Vout)/Lx, where Lx is the distance between the inlet and the outlet.
C. Calculation of electric and hydrodynamic forces
Equations (1)–(5) are solved along with the appropriate boundary conditions to find the electrophoretic velocity of the particle under weak applied electric field. Since the unknown nanoparticle velocity is coupled directly to the wall boundary conditions, we vary the particle velocity until the drag force on the particle balances the electric force. For a nonconductive particle surface with specified zeta potential, we apply the simplified electric force expression:35
In addition, the hydrodynamic force (drag force) on the particle can be obtained by
where |${\bm {\hat \sigma }}$| is the dimensionless stress tensor. Note that a factor 2/(κa)2 appears in (8) due to the scales used to make the velocity and force nondimensional. After we obtain the equilibrium particle velocity Up, the dimensionless mobility can be calculated by
D. Numerical simulation
The equations and associated boundary conditions are solved in a 3D finite element domain using COMSOL V4.3a (COMSOL, Inc., Stockholm, Se). Quadratic elements are used for electric potential and ionic potential, and linear elements are used for velocity and pressure fields. The mesh is refined adaptively in the Debye length near the surfaces of the particle and walls. Mesh independence is checked in all cases, and the relative tolerance is chosen to be 0.001.
Our computational domain is a rectangular channel (length Lx, width Ly, and height Lz) with a spherical particle as shown in Figure 1. Symmetry is used to reduce the computational domain to one-half of the channel. To simulate a practical nanochannel, we consider the confinement effect only in the channel height direction (Lz = 2h), while the channel width is computationally chosen to be large, Ly = 10Lz, and has negligible influence. The ratio of the channel half height to particle radius is defined as |$\hat h \equiv h/a$|. An electric field of 100 V/m is applied along the channel, and KCl is used as the background electrolyte. Zeta potentials are specified at the surface of the particle and walls, where zero corresponds to uncharged condition. We use viscosity and permittivity of water at 20 °C and assume diffusion coefficient of both K+ and Cl− to be 1.96 × 10−9m2s−1.36
III. RESULTS AND DISCUSSIONS
A. Numerical model verification
To validate our numerical model, we compared our results of electrophoretic mobility with literature values26 in an unbounded domain for both low and high particle zeta potentials. Figure 2 shows the dimensionless mobility |$\hat \mu$| as a function of |$\hat \zeta _p$| for different values of κa. Here we use a 50 nm-diameter particle in a 2.5 μm-wide uncharged channel to minimize influence of channel walls and to simulate an unbounded environment. Solid lines represent the current simulation, while dashed lines are adopted from O’Brien and White.26 Note that according to the definition in the literature, the mobility of thick EDL approaches 1 instead of 2/3. Therefore the mobility data from Ref. 26 need to be multiplied by a factor 2/3 to be consistent with our definition of dimensionless mobility. The good agreement indicates that our model correctly captures interactions between electric and flow fields over a wide range of both |$\hat \zeta _p$| and κa.
B. Effect of particle location
In order to obtain useful results from simulations with the particle fixed at the center of the channel, it is important to make sure that the effect of the particle being located at various positions across the channel is minimal. For low electrolyte concentrations, we show that the electric double layers are large and confine the particle to the center of the channel. For high electrolyte concentrations, the electric double layers are thinner and the particle position is distributed across the channel. However, even in this case, we can show that the mobility of the particle within this distribution of transverse channel positions varies only by ∼1%.
Specifically, in our case the particle and channel walls are both negatively charged, if a particle deviates from the midplane of the channel, it will be repelled due to the electrostatic force. However, this force may be screened at least partially by the EDLs, no matter how thick or thin. Therefore, we calculated the electric force in the transverse direction for different particle locations with different electrolyte concentrations. Figure 3 shows the relative electric force as a function of particle location zp for a 50 nm-diameter and a 10 nm-diameter particle in a 100 nm channel, where zp = 0 corresponds to the nanochannel midplane. The relative electrostatic repulsive force is defined as |$\tilde F_{e,z} = F_{e,z} /F_{e,x}$|, which is a comparison of the cross-channel to the driving electric force. When |$\tilde F_{e,z} = 1$| the magnitude of the repulsive force is the same as the axial driving force, which means the particle has strong tendency to move towards the midplane. For low electrolyte concentrations (∼1 mM) the repulsive force is significant, thereby confining the particle very close to the midplane due to the thick EDLs. In contrast, for high electrolyte concentrations (∼3 M), the thin EDLs effectively screen the repulsive force until the particle is close to the wall. Therefore for high electrolyte concentration (thin EDL) the particle location is distributed across the channel.
Subsequently, we focused our attention on the highest electrolyte concentration, 3 M, as a worst case scenario where the widest particle distribution occurs. The electrostatic repulsive force distributionis integrated it to obtain the energy barrier Ue, where Ue ≡ 0 at the midplane due to zero repulsive force. In addition, the mobility of the particle is calculated as a function of transverse location. The probability distribution of particle position is related to the energy barrier, fY(y) ∼ exp [ − Ue(y)/kBT]. The zeta potentials are chosen as ζp = −1 mV and ζw = −2 mV. The average dimensionless mobility for a 50 nm-diameter and a 10 nm-diameter particle is 0.0396 and 0.043, respectively. The mobilities of the two particles at the centerline are 0.0399 and 0.0431, which indicates the deviation is within 1% for both cases. Since the velocity field generated by electro-osmosis behaves like plug flow, the mobility remains relatively unchanged unless the particle is located near the wall. However, the repulsive force confines the particle near the centerline, which decreases substantially the probability of a particle being located near the wall. Based on these results, we are motivated to restrict our attention to analysis of particle mobility at the channel centerline.
C. Particle mobility
In particle mobility experiments, particle total mobility can be separated into electrophoretic and electro-osmotic components, which are usually treated independently.37,38 However, EDLs generated by the particle and by channel walls interact with each other, which affects mobility especially when channel height is comparable to thickness of EDL. Therefore, particle mobility should in fact deviate from that predicted by superposition of electrophoretic and electro-osmotic mobility. Further restricting our attention to a 50 nm-diameter particle, Figure 4 shows particle total mobility from numerical simulation as a function of κa for |$\hat h$| = 50 and 2. This result is compared to total mobility obtained by superimposing electrophoretic and electro-osmotic mobility (|$\hat \mu _{total} = \hat \mu _{EP} - \hat \mu _{EOF}$| as in classical models), using the classic expression |$\hat \mu _{EOF} = \hat \zeta _w$| to calculate electro-osmotic flow. The particle and walls are assumed to be negatively charged to simulate practical conditions,22,23 and values of zeta potential are chosen as |$\hat \zeta _p = \hat \zeta _w = - 0.04$| to avoid any nonlinear effects caused by double layer polarization. For thin EDL systems in a large channel (|$\hat h$| = 50), the results from superposition are valid, since no interaction occurs between the EDLs of the channel walls and the particle. However, as κa decreases, the EDLs overlap, and lead to a moderate deviation at thick EDL region.
When the channel height further reduces (|$\hat h$| = 2), deviation in mobility exists even for thin EDLs. This indicates that charged walls influence the confinement effect by increasing hydrodynamic drag. Near the thick EDL limit, on the other hand, a significant deviation can be observed. The speed of background flow (electro-osmosis) is strongly suppressed by overlapping EDLs, which provides a lower magnitude of |$\hat \mu _{EOF}$| compared to the expression (|$\hat \mu _{EOF} = \hat \zeta _w$|) used in the superposition method. Therefore, the total mobility calculated by superposition overpredicts mobility. Given the influence of overlapping EDLs on electrophoresis and electro-osmosis, we conclude that incorporating both effects simultaneously is required in order to accurately predict the mobility of a particle in a nanochannel.
Figure 5 shows total mobility verses κa for two different size channels (solid lines for |$\hat h$| = 50, dashed lines for |$\hat h$| = 2). The value of wall zeta potential |$\hat \zeta _w$| is fixed at −0.08, and particle zeta potential |$\hat \zeta _p$| is chosen as −0.04, −0.08, and −0.12. First, in the thin EDL regime (κa > 100), there is no interaction between EDLs from the particle and walls. The Helmholtz-Smoluchowski model is valid and total mobility can be approximated by superposition of electrophoretic and electro-osmotic mobility, that is, |$\hat \mu _{total} \approx \hat \mu _{EP} - \hat \mu _{EOF} = ( {\hat \zeta _p - \hat \zeta _w } )$|. Note that both the particle and walls are negatively charged, and consequently electrophoresis is in the opposite direction to electro-osmosis. For example, if |$\hat \zeta _p$| = −0.12, electrophoretic mobility dominates and the particle moves against direction of the background EOF flow (|$\hat \mu _{total}$|< 0).
The deviation in mobility between the two channel heights depends on the value of |$( {\hat \zeta _p - \hat \zeta _w } )$| at the thin EDL limit. In the case of |$\hat \zeta _p = \hat \zeta _w$| = −0.04, mobility approaches zero for both channel heights, and the confinement effect is not observed. This is due to the charged condition at the wall surfaces. Specifically, electro-osmotic flow generated in thin EDL systems can be interpreted as a shear flow, which superimposes shear drag on the particle. In this case, the velocity near the walls has the same magnitude and direction as velocity near the particle. Consequently, no shear stress from the walls can be perceived by the particle, which effectively eliminates the confinement effect. Similarly, if flow near the walls is faster, for example in the |$\hat \zeta _p$| = −0.04 case, the shear stress acts on the particle along the direction of electrophoresis and increases electrophoretic mobility, which results in a lower total mobility. Since shear stress depends on the velocity gradient, this effect is predominant for small channel heights. The results suggest that particle mobility can be altered by using different materials or surface coatings for nanochannels.
As κa decreases (λD increases), |$\hat \mu _{total}$| increases due to decreased |$\hat \mu _{EP}$|. As the EDL thickness increases, the wall EDLs will increasingly interact with the particle EDLs. This reduces surface charge density on the particle, which reduces |$\hat \mu _{EP}$|. In addition, the magnitude of |$\hat \mu _{EOF}$| remains constant, since the EDLs from the top and bottom walls do not overlap. Subsequently, total mobility continues to increase with decreasing κa. As κa decreases even further, the EDLs from the top and bottom walls start to interact, thereby reducing the electro-osmotic flow, as observed in Figure 5. Overlapping can occur for large EDL thicknesses, even in relatively large channels. Consequently, the maximum mobility appears at lower κa for |$\hat h$| = 50.
As κa approaches zero (thick EDL limit), electro-osmotic mobility approaches zero due to extremely low charge density in flow. In this case, particle total mobility is dominated by electrophoretic mobility. According to Huckel-Onsager's model, the drag force on the particle can be predicted by Stokes’ law (with modification to include channel walls). In addition, since κa approaches zero, the Poisson-Boltzmann equation (1) reduces to Laplace's equation. With specified zeta potential at the particle and walls, surface charge density of the particle should be proportional to the difference in the zeta potentials. The electric force as well as total particle mobility should proportional to (|$\hat \zeta _p - \hat \zeta _w$|) near the thick EDL limit.
At low κa, the effect of EDL overlapping leads to significant differences in mobility. The mobility of a nanoparticle in a nanochannel is greater than that in a microchannel at moderate κa. Finally, when magnitude of particle zeta potential is higher (|$\hat \zeta _p$| = −0.12), mobility changes from negative to positive value, and then back to negative again. The results indicate that particle motion can be manipulated by changing electrolyte concentration in a nanochannel.
D. Application to nanochannels
In the above sections we showed that particle mobility is affected by particle zeta potential, electrolyte concentration, channel size, and wall zeta potential. However, particle zeta potential and electrolyte concentration cannot be arbitrarily specified in real applications. Consequently, particle mobility can be manipulated by choosing channels with different heights or surface charges.39 For example, if the value of wall zeta potential is chosen between the values of particle zeta potential of two types of particles, the two particle types will migrate in the opposite directions under an applied electric field. The results shown in Figure 5 can also be viewed as particle mobilities with different sizes in the same charged channel. For example, consider two particle types: 50 nm-diameter and 10 nm-diameter. Figure 6 shows particle mobility as a function of electrolyte concentration in a microchannel and a nanochannel. The zeta potentials are chosen as |$\hat \zeta _p$| = −0.04, |$\hat \zeta _w$| = −0.08. Figure 6 indicates how nanochannels (dashed lines) can be used to separate 10 and 50 nm diameter particles using 1 μM–1 M electrolyte concentrations, with a significant improvement being observed at the lower and higher electrolyte concentration regions. The microchannels (solid lines) can separate the two types of particles, but only for higher electrolyte concentrations. This indicates that one can exploit nanochannels to separate particles efficiently over a large range of electrolyte concentrations.
The conventional method to estimate particle zeta potential incorporates Helmholtz-Smoluchowski's model and electrophoretic mobility, which is calculated by subtracting particle total mobility from the electro-osmotic mobility. As discussed above, this approach fails for cases with thick EDLs or with high particle zeta potentials. In addition, the difference between total mobility and electro-osmotic mobility is relatively small and can introduce significant error. Subsequently, we propose here a new method for estimating particle and wall zeta potentials by using only the total particle mobility as measured in two different height channels. Figure 7 shows contours of total particle mobility as a function of |$\hat \zeta _p$| and |$\hat \zeta _w$| with 1 mM electrolyte concentration (κa = 2.61) for |$\hat h$| = 50 (solid lines) and |$\hat h$| = 2 (dashed lines). Numbers on contour lines correspond to dimensionless total particle mobility. Once the total particle mobility has been measured in the two different height channels, the intersection of corresponding contours leads to definitive values of particle and wall zeta potentials.
For example, if we want to determine the zeta potential of a 50 nm-diameter particle, we can measure its mobility in a 2.5 μm-high channel (|$\hat h = 50$|) and a 100 nm-high channel (|$\hat h = 2$|) filled with 1 mM KCl solution (κa = 2.61). For demonstration purposes, let's assume that the measured total particle mobilities are, say, |$\mu _{\hat h = 50} = 1$| and |$\mu _{\hat h = 2} = 1$|. Then by comparing the measured mobilities with Figure 7, we can estimate that the particle and channel zeta potentials are 19.8 mV (|$\hat \zeta _p$| = 0.792) and 39.8 mV (|$\hat \zeta _w$| = 1.59), respectively. If the error from mobility measurement is, say, 3%, the error for zeta potential would also be about 3%. Our approach can avoid the errors generated from measuring electro-osmotic mobility, and works for a wide range of κa and |$\hat \zeta _p$|. Note that Figure 7 is restricted to the value κa = 2.61. This approach can be readily extended for other values of κa.
IV. CONCLUSION
Numerical simulation was used to investigate the electrophoretic mobility of a spherical particle in a confined nanochannel. The numerical model was validated for a wide range of zeta potentials, electrolyte concentrations, and channel sizes. The results indicate that, for the cases under current consideration, the mobility of a particle located at the centerline of the nanochannel agrees to within 1% of the average mobility for a particle distributed transversely throughout the nanochannel.
When a nanoparticle is confined in a nanochannel, overlapping EDLs between the charged particles and nanochannels walls can be important. Charged walls not only induce a background flow (electro-osmotic flow), but also affect the particle's hydrodynamic drag and surface charge density. Particle mobility in nanochannels can be greater than that in microchannels, if the electrolyte concentration is chosen properly.
The numerically simulated results indicate that different size nanoparticles may be electrophoretically separated using nanochannels, even if the particles have similar zeta potentials. For example, a 100 nm-wide nanochannel could be used to separate 50 nm and 10 nm-diameter particles for a wide range of electrolyte concentrations. Finally, a new method is proposed for determining zeta potentials of the particle and channel walls by measuring the mobility of a particle using two different height channels. This method can avoid errors generated from measuring electro-osmotic velocity and it is applicable over a wide range of zeta potentials.
ACKNOWLEDGMENTS
This work is supported in part by the Institute for Collaborative Biotechnologies through grant W911NF-09-0001 from the U.S. Army Research Office. The content of the information does not necessarily reflect the position or the policy of the Government and no official endorsement should be inferred.