Transport coefficients for electrolytes in arbitrarily shaped nano and micro-fluidic channels

We consider laminar flow of incompressible electrolytes in long, straight channels driven by pressure and electro-osmosis. We use a Hilbert space eigenfunction expansion to address the general problem of an arbitrary cross section and obtain general results in linear-response theory for the hydraulic and electrical transport coefficients which satisfy Onsager relations. In the limit of non-overlapping Debye layers the transport coefficients are simply expressed in terms of parameters of the electrolyte as well as the geometrical correction factor for the Hagen-Poiseuille part of the problem. In particular, we consider the limits of thin non-overlapping as well as strongly overlapping Debye layers, respectively, and calculate the corrections to the hydraulic resistance due to electro-hydrodynamic interactions.


Introduction
Laminar Hagen-Poiseuille and electro-osmotic flow is important to a variety of labon-a-chip applications and microfluidics [1,2,3] and the rapid development of micro and nano fabrication techniques during the past decade has put even more emphasis on flow in channels with a variety of shapes depending on the fabrication technique in use. The list of examples includes rectangular channels obtained by hot embossing in polymer wafers, semi-circular channels in isotropically etched surfaces, triangular channels in KOH-etched silicon crystals, Gaussian-shaped channels in laser-ablated polymer films, and elliptic channels in stretched PDMS devices [4]. While general results for the shape-dependence of the hydraulic resistance in the case of a nonconducting fluid were reported recently [5] there has, according to our knowledge, been no analogous detailed study of the shape-dependence of flow of electrolytes in the presence of a zeta potential which is a scenario of key importance to lab-on-achip applications involving biological liquids/samples in both microfluidic [6,7,8] and nanofluidic channels [9,10,11,12,13,14].
In this work we address the general steady-state flow problem in Fig. 1 where pressure gradients and electro-osmosis (EO) are playing in concert [15]. We consider a long, straight channel of length L having a constant cross section Ω of area A and boundary ∂Ω of length P. The channel contains an incompressible electrolyte, which we for simplicity assume to be binary and symmetric, i.e., containing ions of charge +Ze and −Ze and equal diffusivities D. The electrolyte has the Debye screening length λ D , bulk conductivity σ o , viscosity η, permittivity ǫ and at the boundary ∂Ω it has a zeta potential ζ. The laminar, steady-state flow is driven by a linear pressure drop ∆p and a linear voltage drop ∆V . With these definitions flow will be in the positive x direction. In the linear-response regime the corresponding volume flow rate Q and charge current I are related to the driving fields by where, according to Onsager relations [16], G is a symmetric, G 12 = G 21 , two-by-two conductance matrix. The upper diagonal element is the hydraulic conductance or Figure 1. A translation invariant channel of arbitrary cross section Ω of area A containing an electrolyte driven by a pressure gradient −∆p/L and by electroosmosis through the potential gradient −∆V /L. The channel wall ∂Ω has the electrical potential ζ, which induces a thin, charged Debye layer (dark gray) that surrounds the charge neutral bulk (light gray).
inverse hydraulic resistance given by where α is the dimensionless geometrical correction factor, shown in Ref. [5] to be a linear function of the dimensionless compactness parameter C = P 2 /A. While there is no intrinsic length scale influencing G 11 , the other elements of G depend on the Debye screening length λ D . This length can be comparable to and even exceed the transverse dimensions in nano-channels [9,10,11], in which case the off-diagonal elements may depend strongly on the actual cross-sectional geometry. However, for thin Debye layers with a vanishing overlap all four matrix elements in G are independent of the details of the geometry. For a free electro-osmotic flow, a constant velocity field v eo = (ǫζ/η)∆V /L is established throughout the channel, except for in the thin Debye layer of vanishing width. Hence Q = v eo A and From Ohm's law I = (σ o A/L)∆V it follows that For strongly overlapping Debye layers we find, see Sections 3.2.2 and 4.2 that We emphasize that the above results are generally valid for symmetric electrolytes, even beyond the Debye-Hückel approximation. In the Debye-Hückel limit Zeζ ≪ k B T they also hold for asymmetric electrolytes. We also note that in the Debye-Hückel limit the expressions agree fully with the corresponding limits for a circular cross section and the infinite parallel plate system, were explicit solutions exist in the Debye-Hückel limit in terms of Bessel functions [17,18] and cosine hyperbolic functions [18], respectively. From the corresponding resistance matrix R = G −1 we get the hydraulic resistance where β ≡ G 12 G 21 /(G 11 G 22 ) is the Debye-layer correction factor to the hydraulic resistance. In the two limits we have For ζ going to zero β vanishes and we recover the result in Ref. [5] for the hydraulic resistance.

Governing equations
For the system illustrated in Fig. 1, an external pressure gradient ∇p = −(∆p/L) e x and an external electrical field E = E e x = (∆V /L) e x is applied. There is full translation invariance along the x axis, from which it follows that the velocity field is of the form v( r) = v( r ⊥ ) e x where r ⊥ = y e y + z e z . For the equilibrium potential and the corresponding charge density we have φ eq ( r) = φ eq ( r ⊥ ) and ρ e eq ( r) = ρ e eq ( r ⊥ ), respectively. We follow our related recent work [19] and use the Dirac bra-ket notation [20,21], where functions f ( r ⊥ ) in Ω are written as f with inner products defined by the cross-section integral From the Navier-Stokes equation it follows that the velocity is governed by [22,23] where ∇ 2 ⊥ = ∂ 2 y + ∂ 2 z is the 2D Laplacian. The equilibrium potential φ eq and the charge density ρ e eq are related by the Poisson equation The velocity v is subject to a no-slip boundary condition on ∂Ω while the equilibrium potential φ eq equals the zeta potential ζ on ∂Ω. Obviously, we also need a statistical model for the electrolyte, and in the subsequent sections we will use the Boltzmann model where the equilibrium potential φ eq is governed by the Poisson-Boltzmann equation. However, before turning to a specific model we will first derive general results which are independent of the description of the electrolyte. We first note that because Eq. (10) is linear we can decompose the velocity as v = v p + v eo , where v p is the Hagen-Poiseuille pressure driven velocity governed by and v eo is the electro-osmotic velocity given by The latter result is obtained by substituting Eq. (11) for ρ e eq in Eq. (10). The upper diagonal element in G is given by G 11 = 1 v p /∆p which may be parameterized according to Eq. (2). The upper off-diagonal element is given by G 12 = 1 v eo /∆V and combined with the Onsager relation we get where we have used that 1 1 = A and introduced the average potentialφ eq = φ eq 1 / 1 1 . There are two contributions to the lower diagonal element G 22 ; one from migration, G mig 22 = 1 σ /L, and one from electro-osmotic convection of charge, G conv

22
= ρ e eq v eo /∆V , so that where the electrical conductivity σ( r ⊥ ) depends on the particular model for the electrolyte. For thin non-overlapping Debye layers we note thatφ eq ≃ 0 so that Eq. (14) reduces to Eq. (3) and, similarly since the induced charge density is low, Eq. (15) reduces to Eq. (4). For strongly overlapping Debye layers the weak screening means that φ eq approaches ζ so that the off-diagonal elements G 12 = G 21 and the G conv 22 part of G 22 vanish entirely. In the following we consider a particular model for the electrolyte and calculate the asymptotic suppression as a function of the Debye screening length λ D .

Debye-Hückel approximation
In the Debye-Hückel approximation the equilibrium potential φ eq is governed by the linearized Poisson-Boltzmann equation [3] where λ D is the Debye screening length. The validity of this model will be discussed in more detail in Sec. 4.

Hilbert space formulation
In order to solve Eqs. (10), (11), and (16) we will take advantage of the Hilbert space formulation [24], often employed in quantum mechanics [21], and recently employed by us on the problem of an accelerating Poiseuille flow [19]. The Hilbert space of real functions on Ω is defined by the inner product in Eq. (9) and a complete, countable set ψ n of orthonormal basis functions, i.e., where δ nm is the Kronecker delta. We choose the eigenfunctions ψ n of the Helmholtz equation (with a zero Dirichlet boundary condition on ∂Ω) as our basis functions, With this complete basis any function in the Hilbert space can be written as a linear combination of basis functions. In the following we write the fields as v = ∞ n=1 a n ψ n , Inserting Eqs. (18) and (20) into Eq. (16), and multiplying by ψ m , yields b n = ζ ψ n 1 1 + (κ n λ D ) 2 , n = 1, 2, 3, . . .

Transport coefficients
The flow rate and the electrical current are conveniently written as where the second relation is the linearized Nernst-Planck equation with the first term being the convection/streaming current while the second is the ohmic current. Substituting Eqs. (19) and (21) into these expressions we identify the transport coefficients as where is the effective area of the eigenfunction ψ n . The fraction A n /A is consequently a measure of the relative area occupied by ψ n satisfying the sum-rule ∞ n=1 A n = A [19]. We note that as expected G obeys the Onsager relation G 12 = G 21 . Furthermore, using that we get the following bound between the off-diagonal elements G 12 = G 21 and the lower diagonal element G 22 , In the context of the geometrical correction factor α studied in detail in [5] we note that the first diagonal element may be written as is a characteristic hydraulic conductance and the geometrical correction factor α can be expressed as [19] α where k n = κ n A/P is a dimensionless eigenvalue. In passing we furthermore note that this formal result is a convenient starting point for perturbative analysis of the correction due to small changes in the boundary ∂Ω [25].
3.2.1. Non-overlapping, thin Debye layers. For the off-diagonal elements of G we use that . In Section 5 we numerically justify that the smallest dimensionless eigenvalue k 2 1 is of the order unity, so we may approximate the sum by a factor of unity, see Table 1, whereby we arrive at Eq. (3) for λ D ≪ A/P. These results for the off-diagonal elements are fully equivalent to the Helmholtz-Smoluchowski result [18]. For G 22 we use that (κ n λ D ) 2 , thus we may neglect the second term, whereby we arrive at Eq. (4).

Strongly overlapping Debye layers.
In the case of κ 1 λ D ≫ 1 we may use the This is the Debye-Hückel limit of Eq. (5) where I n is the nth modified Bessel function of the first kind, and were we have explicitly introduced the variable A/αλ 2 D to emphasize the asymptotic dependence in Eq. (5) for strongly overlapping Debye layers. We note that we recover the limits in Eqs. (3) and (5) for λ D ≪ A/P and λ D ≫ A/P, respectively.

Beyond the Debye-Hückel approximation
In order to go beyond the Debye-Hückel approximation we consider, for simplicity, a symmetric binary (Z:Z) electrolyte. Next, we neglect strong correlations between the ions so that the equilibrium properties of the electrolyte are governed by Boltzmann statistics [3], i.e., the concentrations of the two type of ions are given by This is equivalent to assuming equilibrium with bulk reservoirs at the ends of the channel in which the potential φ eq tends to zero and both concentrations c ± eq to c o . Substituting the charge density ρ e eq = Ze(c + eq − c − eq ) into the Poisson equation (11) we arrive at the nonlinear Poisson-Boltzmann equation [3,18] where the Debye screening length is given by The conductivity σ of the electrolyte depends on the local ionic concentrations assuming equal diffusivities D for the two ionic species. In the Debye-Hückel limit, We calculate the off-diagonal elements from G 12 = G 21 = ρ e eq v p /∆p and find Similarly, Eq. (15) for the two components in the electrical conduction G 22 = G mig 22 + G conv 22 we get where we have used that σ o = ǫD/λ 2 D and introduced the dimensionless quantity m, which indicates the importance of electro-osmosis relative to electro-migration.

Non-overlapping, thin Debye layers
In the limit of thin Debye layers we have already discussed how Eq. (14) in general leads to Eq. (3) because the screening is good and φ eq is nonzero only on a negligible part of Ω. This property is more implicit when G 12 or G 21 is written in the form of Eq. (42), which is a more appropriate starting point for analyzing the limit of strongly overlapping Debye layers. For G 22 the calculations are more involved; we assume that the channel wall is sufficiently smooth on the Debye-length scale so that we can everywhere use the Gouy-Chapman (GC) solution for a semi-infinite planar geometry [3,18], where r n denotes the normal distance to the channel wall. Substituting this into Eq. (43) and (44) the integrals can be carried out analytically resulting in where Du is the Dukhin number defined as the ratio of the surface conductivity in the charged Debye layers to the bulk conductivity σ o times the geometrical length scale A/P (see Ref. [28] and references therein). Clearly, when the Debye layer becomes very thin, surface conduction is negligible and we recover the simple result in Eq. (4).

Strongly overlapping Debye layers
When the Debye layers are strongly overlapping the screening is weak and φ eq ≈ ζ throughout the cross section. Hence we can pull the integrand sinh(Zeφ eq /k B T ) outside the bra-ket in Eq. (42) and we arrive at Eq. (5). Here, we have used that 1|v p /∆p = G 11 and introduced the parameterization in Eq. (2). Similarly, from Eqs. (43) and (44) we obtain Eq. (6) where we have used Eqs. (14) and (5) to eliminate 1|ζ − φ eq . We note that due to shifts in free energies, the zeta potential inside a narrow channel with significant Debye-layer overlap is generally not the same as in a macroscopic channel with no overlap, see e.g. [11,29] for a discussion.

The Helmholtz basis
Only few geometries allow analytical solutions of both the Helmholtz equation and the Poisson equation. The circle is of course among the most well-known solutions and the equilateral triangle is another example. However, in general the equations have to be solved numerically, and for this purpose we have used the commercially available finite-element software Femlab [26]. The first eigenstate of the Helmholtz equation is in general non-degenerate and numbers for a selection of geometries are tabulated in Table 1. Note how the different numbers converge when going through the regular polygons starting from the equilateral triangle through the square, the regular pentagon, and the regular hexagon to the circle. In general, k 2 1 is of the order unity, and for relevant high-order modes (those with a nonzero A n ) the eigenvalue is typically much larger. Similarly, for the effective area we find that A 1 /A ≤ 4/γ 2 1 ≃ 0.69 and consequently we have A n /A < 1 − 4/γ 2 1 ≃ 0.31 for n ≥ 2. The transport coefficients in Eqs. (27) to (30) are thus strongly influenced by the first eigenmode which may be used for approximations and estimates of the transport coefficients. As an example the column for α/C is well approximated by only including the first eigenvalue in the summation in Eq. (34).

Transport coefficients
Our analytical results predict that when going to either of the limits of thin nonoverlapping or strongly overlapping Debye layers, the transport coefficients only   depend on the channel geometry through the cross sectional area A and the correction factor α. Therefore, when plotted against the rescaled Debye length α/A λ D , all our results should collapse on the same asymptotes in the two limits.
In Fig. 2 we show the results for the off-diagonal coefficients obtained from finiteelement simulations in the Debye-Hückel limit Zeζ ≪ k B T for three different channel cross sections, namely two parabola shaped channels of aspect ratio 1:1 and 1:5, respectively, and a rectangular channel of aspect ratio 1:5. In all cases we find excellent agreement between the numerics and the asymptotic expressions. For the comparison we have also included exact results, Eq. (36), for the circular cross section as well as results based on only the first eigenvalue in Eq. (28). Even though Eq. (36) is derived for a circular geometry we find that it also accounts remarkably well for even highly non-circular geometries in the intermediate regime of weakly overlapping Debye layers.
In Fig. 3 we show numerical results for the off-diagonal transport coefficients beyond the Debye-Hückel approximation. At large zeta potentials the Debye layer is strongly compressed and the effective screening length reduced. Therefore the suppression of the electro-osmotic flow/streaming current at strong Debye-layer overlap is shifted to larger values of λ D as compared to the Debye-Hückel limit. For the comparison we have also included the exact result for a circular cross section in the Debye-Hückel approximation as well as the asymptotic expression for non-overlapping and strongly overlapping Debye layers, Eq. (5). As seen the asymptotic expressions account well for the full numerical solutions independently of the geometry.
In Fig. 4 we show numerical results for the electrical conductance beyond the Debye-Hückel approximation. Open symbols show the electro-migration part G mig 22 , subtracted the trivial bulk contribution σ o A/L, whereas solid symbols show the  (44), respectively. Again, we find that the numerics are in excellent agreement with our asymptotic results. For the λ D ≪ P/A regime we note that the Dukhin number, Eq. (48), is proportional to λ D P/A = C/A λ D and not α/A λ D . Therefore, strictly one would not in general expect data to collapse on the same asymptote for λ D ≪ P/A. Looking carefully at this part of the figure, small variations can be seen from geometry to geometry -the reason why the variations are still so small is that α/C ∼ 2 independently of geometry, see Table 1.

Conclusion
We have analyzed the flow of incompressible electrolytes in long, straight channels driven by pressure and electro-osmosis.
By using a powerful Hilbert space eigenfunction expansion we have been able to address the general problem of an arbitrary cross section and obtained general results for the hydraulic and electrical transport coefficients. Results for strongly overlapping and thin, non-overlapping Debye layers are particular simple, and from these analytical results we have calculated the corrections to the hydraulic resistance due to electro-hydrodynamic interactions. These analytical results reveal that the geometry dependence only appears through the area A and the correction factor α, as the expressions only depend on the rescaled Debye length α/A λ D . Our numerical analysis based on finite-element simulations indicates that these conclusions are generally valid also for intermediate values of λ D . Combined with recent detailed work on the geometrical correction factor [5] the present results constitute an important step toward circuit analysis [15] of complicated micro and nanofluidic networks incorporating complicated cross-sectional channel geometries.