A nonlinear equation for ionic diffusion in a strong binary electrolyte

A nonlinear equation for ionic diffusion in a strong binary electrolyte

Abstract

The problem of the one dimensional electro-diffusion of ions in a strong binary electrolyte is considered. In such a system the solute dissociates completely into two species of ions with unlike charges. The mathematical description consists of a diffusion equation for each species augmented by transport due to a self consistent electrostatic field determined by the Poisson equation. This mathematical framework also describes other important problems in physics such as electron and hole diffusion across semi-conductor junctions and the diffusion of ions in plasmas. If concentrations do not vary appreciably over distances of the order of the Debye length, the Poisson equation can be replaced by the condition of local charge neutrality first introduced by Planck. It can then be shown that both species diffuse at the same rate with a common diffusivity that is intermediate between that of the slow and fast species (ambipolar diffusion). Here we derive a more general theory by exploiting the ratio of Debye length to a characteristic length scale as a small asymptotic parameter. It is shown that the concentration of either species may be described by a nonlinear integro-differential equation which replaces the classical linear equation for ambipolar diffusion but reduces to it in the appropriate limit. Through numerical integration of the full set of equations it is shown that this nonlinear equation provides a better approximation to the exact solution than the linear equation it replaces.

pacs:
82.45.Gj,82.45.-h,82.70.y,47.57.eb,47.57.jd,72.20.Dp

The one dimensional electro-diffusion equations describing the evolution of ionic concentrations , in a fully dissociated binary electrolyte with a self consistent electrostatic potential may be put in the form (1)

 ∂tn−∂x(n∂xϕ) = ∂xxn (1) u−1∂tp−z∂x(p∂xϕ) = ∂xxp (2) ∂xxϕ = −n−zp (3)

where is the ratio of mobilities of species p and n and is the ratio of their valences. We will take n as the species with lesser mobility, so that . Since the two species are oppositely charged, is negative. The above equations have been rendered dimensionless by scaling the concentrations by the concentration of n at infinitely far points () so that , the co-ordinate by a length scale related to Debye length, the time by a diffusion time and the potential by the thermal energy . In the above, and are the diffusivities of the two species which are related to the mobilities through the Einstein relation , where is the Boltzmann constant and the absolute temperature. The length scale referred to above is defined by the expression in terms of the electronic charge (), permittivity of vacuum () and dielectric constant (). It is related to the Debye length in a homogeneous charge neutral electrolyte where the less mobile species has a dimensionless concentration of as .

Eq.(1)-(3) and its modifications also describe other transport problems of interest in various areas of physics. For example, if where is an external field and is the perturbation in the potential and if a term representing the density of free charges is added to the right hand side of Eq. (3), then these equations become the Van Roosbroeck equations describing the migration of charge carriers in solid state physics, where , and represent the concentrations of electrons, holes and dopants (2). If is interpreted as the density of ion exchange sites within a membrane, then the same set of equations also describes various filtration processes such as electrodialysis (1). In plasma physics (3) one often encounters situations where the electrons and positive ions may be regarded as two interpenetrating charged fluids of different diffusivities coupled by a self-consistent electric field. It is then described mathematically by Eq. (1)-(3). The generalization of these equations to three or more species in the presence of an applied electric field describe the separation of charged macromolecules in capillary electrophoresis (4) and other separation processes based on electric charge.

One of the simplest experimental realizations of Eq. (1)-(3) is the “Liquid Junction” problem, where a barrier initially separates two semi-infinite regions of a binary electrolyte (such as sodium chloride in water) of different concentrations, and subsequently, at time , the barrier is removed. Since generally the positive and negative components would diffuse at different rates a polarization vector and consequently a measurable “Liquid Junction Potential” (LJP) appears across the interface. Planck (5) derived an expression for the LJP on the basis of the local charge neutrality assumption, an idea that has since become a central paradigm in many problems involving electro-diffusion. The physical basis of the approximation is that the electrostatic force that would arise if a charge separation was realized is so strong that in practice it precludes any departure from electro-neutrality anywhere in space.

Mathematically, the local electro-neutrality assumption amounts to neglecting the term on the left hand side of Eq. (3), so that it reduces to . This relation can now be used to eliminate the terms involving in Eq. (1) and (2). It then follows that (as well as ) obeys the diffusion equation

 ∂tn=D∂xxn (4)

where , a result due to Henderson (6); (7). Since , clearly . The coulomb attraction between the two species enhances the rate of spreading of the slower species and reduces that of the faster species so that both diffuse at an equal, ambipolar (that is, the same for either charge) rate. In the analysis outlined above, though the term in is dropped in Eq. (3), it must be retained in Eq. (1)-(2), a situation that appeared self contradictory and led to some controversy until Hickman (8); (9) provided a proper interpretation within the framework of an asymptotic theory based on expansion in the small parameter , where is a characteristic wave number and is the Debye length. Since the Debye length is normally a very small quantity (about nm in a M solution of a common salt like sodium or potassium chloride) in most applications the charge neutrality assumption is very accurate. However, it could be violated in many modern applications such as in nanochannels where one or more geometric dimensions may be of the order of nanometers. Such departures from charge neutrality may give rise to new effects. It would therefore seem worthwhile to first study these effects in the context of the simple model system represented by Eq. (1)-(3).

In this paper we use the method of multiple scales (10) to reduce Eq. (1)-(3) to a nonlinear one dimensional system in the limit of long (compared to the Debye length) length scales and slow (compared to the diffusion time over a Debye length) time scales. We show that at the lowest order the equation for ambipolar diffusion is recovered. If the asymptotic theory is continued to the next order, a nonlinear integro-differential equation for the concentration emerges. Numerical solutions of this reduced equation is seen to agree better with that of the full system of Eq. (1)-(3) and lead to effects not captured in the lowest order theory based on local charge neutrality.

We would like to study the behavior of Eq. (1)-(3) at long length and time scales, under the boundary conditions that all approach constant values and respect local charge neutrality as . Thus, we are considering passive diffusion in the absence of imposed electric fields and currents. Following the usual procedure (10), we introduce slow variables and and suppose that depend solely on and . Then Eq. (1)-(3) reduce to

 ∂τn−∂ξ(n∂ξϕ) = ∂ξξn (5) u−1∂τp−z∂ξ(p∂ξϕ) = ∂ξξp (6) ϵ∂ξξϕ = −n−zp (7)

where is a small parameter. Expanding all dependent variables in asymptotic series in , such as , substituting in Eq. (5)-(7) and equating the coefficients of on both sides, we get a series of equations starting with the lowest order ones:

 ∂τn0−∂ξ(n0∂ξϕ0) = ∂ξξn0 (8) u−1∂τp0−z∂ξ(p0∂ξϕ0) = ∂ξξp0 (9) −n0−zp0 = 0 (10)

If the last equation is used to eliminate the second terms from the first two, we get

 ∂τn0 = D∂ξξn0 (11) p0 = −n0/z (12)

with , representing ambipolar diffusion and if the time derivative terms are eliminated instead and the result integrated, we get

 ϕ0=u−11−uzlnn0. (13)

If the above equation is evaluated at , we get Planck’s formula (5) for the potential drop across a liquid junction. At the next order in , we have

 ∂τn1−∂ξ(n1∂ξϕ0)−∂ξ(n0∂ξϕ1) = ∂ξξn1 (14) u−1∂τp1−z∂ξ(p1∂ξϕ0)−z∂ξ(p0∂ξϕ1) = ∂ξξp1 (15) ∂ξξϕ0=−n1−zp1 (16)

If we add Eq. (14) and (15) and use Eq.(10), the term in is eliminated. We then obtain after substituting from Eq. (13) and after some algebra;

 ∂τn1 = D∂ξξn1+α∂ξξξξ(lnn0)−β∂ξξ(∂ξlnn0)2 (17) p1 = −n1z−u−1z(1−uz)∂ξξ(lnn0), (18)

where and are positive constants defined by

 α = −uz(u−1)2(1−uz)3 (19) β = u(u−1)(2−z−uz)2(1−uz)3 (20)

If we add Eq. (11) to times Eq. (17) and transform back to our original variables and , we get

 ∂tn=D∂xxn+α∂xxxx(lnn)−β∂xx(∂xlnn)2 (21)

with an error which is higher than order . In arriving at Eq. (21) we have replaced the terms by and by , which is justified because doing so introduces an error in Eq. (21) that is higher than order . Similarly, combining Eq. (18) with Eq. (12) we get

 p=−nz−u−1z(1−uz)∂xx(lnn) (22)

with an error of order . The last term in Eq.(22) is due to the departure from charge neutrality.

Let us now consider some of the consequences of Eq. (21) and (22). First, we consider weak perturbations from the constant values at infinity: and , where . This gives a hyper-diffusion equation for :

 ∂tn′=D∂xxn′+α∂xxxxn′. (23)

Therefore, solutions of the form have growth rates . Since , high wavenumber modes are unstable. We will show that this instability is entirely spurious. Indeed, Eq. (21) is not even valid for modes since these solutions violate the premise . The instability arises as a consequence of truncating the asymptotic series, as can be seen by considering the linearized version of Eq. (1)-(3) with :

 ∂tn′−∂xxϕ = ∂xxn′, (24) ∂xxϕ = ∂xxp′, (25) ∂xxϕ = −n′−zp′. (26)

If Eq. (25) is integrated and the result substituted in Eq. (26) an exact expression for may be found in terms of :

 ϕ = −1z(1+z−1∂xx)−1n′, (27) = −∫+∞−∞G(|x−y|;√−z)n′(y)dy, (28) = −1z[1−z−1∂xx+z−2∂xxxx+⋯]n′ (29)

where is the Green’s function of the Helmholtz operator . If the exact inversion, that is Eq. (28), is substituted in Eq. (24) one gets an integro-differential equation which in Fourier-space yields a growth rate

 a=−k2[1+1k2−z]∼−k2(1−1z)+k4z2+⋯ (30)

The exact formula for indicated by the sign above shows that there is no high wavenumber instability if the exact inversion, Eq. (28), is employed. The second part, indicated by the sign, shows the result of the approximate inversion, Eq.(29), based on the low wavenumber approximation. In this case there is a high wavenumber instability if the asymptotic series is truncated after an even number of terms. That is, the root cause of the instability is the oscillatory approach to the limit of the series indicated on the right hand side of Eq. (30).

The above analysis suggests a simple way of modifying Eq. (21) without violating its asymptotic validity, but at the same time eliminating the spurious high wavenumber instability. We recognize that the term in Eq. (21) most likely resulted from a truncated expansion of the Helmholtz operator: and simply replace in Eq. (21) by . Thus, we get the nonlinear integro-differential equation:

 ∂tn=∂xxF[n] (31)

where is the functional:

 F[n] = Dn−β[∂x(lnn)]2 (32) + 12√α∫+∞−∞e−|x−y|/√αln{n(y)n(x)}dy.

Eq. (31) is asymptotically equivalent to Eq. (21) since they differ only by terms of order higher than but it is free from the spurious high wavenumber instability. By virtue of its construction it also has the property that when , the linearized version of Eq.(31) is the exact solution of Eq.(24)-(26). Thus, when is large and amplitudes are small, the solutions of Eq.(31) match closely the true solution, irrespective of the validity of the assumption of long length scales and slow time scales exploited in the asymptotic theory. The formal condition for the validity of Eq. (31) is , where the bars on top of the variables indicate that they are in dimensional form. In dimensionless notation, the condition of validity reduces to . It should be noted that, Eq. (31) as well as the lower order effective diffusivity approximation fails in the limit of very low concentrations (). This is because the Debye length becomes infinitely large in this limit invalidating our premise of long length scales compared to the local value of the Debye length. Eq. (31) together with Eq. (22) to calculate the concentration of the faster diffusing species are our principal results.

A numerical method employing a fourth order finite difference algorithm for spatial derivatives coupled with a fourth order Runge-Kutta time stepping scheme with fixed grid and step size was used to solve a time dependent electrodiffusion problem. An initial condition was used with and the other parameters were set as and . The problem was solved at three levels of approximation:
(I) Using the ambipolar diffusion equation, Eq.(4).
(II) Using the nonlinear equation, Eq.(31).
(III) Using the full electro-diffusion model Eq. (1)-(3).
It is clear that with (I) the variance would increase linearly with time, but not so in the case of (II). Further, a distribution that is initially Gaussian remains so at future times under (I) but not under (II). Therefore, the excess kurtosis is expected to remain zero at all times in the case of (I) but not in the case of (II). In Figure 1 the quantities and are plotted as a function of time () for (II) and (III). In case of (I), at all and this case is not shown for clarity. It is seen that at sufficiently long times (I), (II) and (III) all approach a common asymptotic limit. However, at shorter times, (II) is in better accord with (III) than (I) is.

The qualitative nature of the time variation of and may also be understood on the basis of Eq. (31). To see this, we use Eq.(29) to expand the integral operator and put Eq. (31) in the form

 ∂tn=D∂xxn+α∂xxxx(lnn)−β∂xx(∂xlnn)2+⋯ (33)

where indicate terms involving sixth, eighth, tenth, …. order derivatives. If we further restrict ourselves to small amplitudes, , with , then , so that the Eq.(33) becomes

 ∂tn′=D∂xxn′+α∂xxxxn′−α2∂xxxx(n′)2−β∂xx(∂xn′)2+⋯ (34)

Here the now include the nonlinear terms of higher order. If we multiply Eq. (34) by and integrate, we can generate ordinary differential equations for the moments , starting with . The neglected sixth and higher derivative terms do not contribute for . Without loss of generality, we can assume so that here and . By combining the equations for the second and fourth moments, we derive

 dσ2dt=2D−2βm0∫+∞−∞(∂xn′)2dx (35)

and

 dγdt=−4γDσ2+12ασ4+⋯ (36)

where indicate the nonlinear terms that we suppress. Now, if , then either (if it is zero initially) or else monotonically decays to zero. When deviations from the ambipolar diffusion limit is considered (), we find that initially increases (when is small enough for the second term to dominate), but as becomes larger, the first term eventually dominates resulting in a decrease in towards zero. Further, Eq.(35) shows, that since the rate of increase of the variance is slightly less than but the deficit in the growth rate becomes progressively smaller at large times. This is indeed what is observed in Figure 1.

In summary, a prototypical problem in electro-diffusion was considered in the limit where the Debye-length was non-zero but nevertheless small compared to a characteristic scale of spatial variation. It was shown that in this limit the concentration can be described by a one dimensional nonlinear integro-differential equation which reduces to the linear diffusion equation describing ambipolar diffusion if all but the leading order terms in the ratio of Debye-length to a characteristic spatial scale are neglected. Numerical solution shows that the nonlinear integro-differential equation provides a better approximation to the true solution. The approach described here could be useful for other problems in electro-diffusion where the Debye length is small but may not be considered negligible in relation to other length scales.
Acknowledgement: Support from the NIH under grant R01EB007596 is gratefully acknowledged.

References

1. I. Rubinstein, Electro-Diffusion of Ions (SIAM, Philadelphia, PA., 1990).
2. P. Markowich, C. Ringhofer, and C. Schmeiser, Semiconductor Equations (Springer-Wien, New York, 2002).
3. N. Krall and A. Trivelpiece, Principles of Plasma Physics (McGraw-Hill, New York, 1973).
4. S. Ghosal, Annu. Rev. Fluid Mech. 38, 309 (2006).
5. M. Planck, Ann.phys.Chem. 39, 161 (1890).
6. D. Henderson, Z.phys.Chem. 59, 118 (1907).
7. D. Henderson, Z.phys.Chem. 63, 325 (1908).
8. H. Hickman, Chem. Eng. Sci. 25, 381 (1970).
9. J. Jackson, J. Phys. Chem. 78, 2060 (1974).
10. A. Nayfeh, Perturbation Methods (Pure and Applied Mathematics) (John Wiley and sons, New York, 2000).
113488