## Abstract

The problem of the one-dimensional electro-diffusion of ions in a strong binary electrolyte is considered. The mathematical description, known as the Poisson–Nernst–Planck (PNP) system, consists of a diffusion equation for each species augmented by transport owing to a self-consistent electrostatic field determined by the Poisson equation. This description is also relevant to other important problems in physics, such as electron and hole diffusion across semiconductor 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 the 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 partial differential equation that provides a better approximation than the classical linear equation for ambipolar diffusion, but reduces to it in the appropriate limit.

## 1. Introduction

The one-dimensional electro-diffusion equations describing the evolution of ionic concentrations *n*(*x*,*t*), *p*(*x*,*t*) in a fully dissociated binary electrolyte with a self-consistent electrostatic potential *ϕ*(*x*,*t*) may be put in the form (Rubinstein 1990)
1.1
1.2and
1.3
where *u*=*u*_{p}/*u*_{n} is the ratio of mobilities of species *p* and *n* and *z*=*z*_{p}/*z*_{n} is the ratio of their valences. We will take *n* as the species with lesser mobility, so that *u*≥1. Since the two species are oppositely charged, *z* is negative. The above equations have been rendered dimensionless by scaling the concentrations by the concentration of *n* at infinitely far points (*n*_{*}) so that , the coordinate *x* by a length scale *λ*_{*} related to the Debye length, the time by a diffusion time and the potential *ϕ* by the thermal energy scale (*k*_{B}*T*)/(*ez*_{n}). In the above, *D*_{n} and *D*_{p} are the diffusivities of the two species that are related to the mobilities *u*_{n},*u*_{p} through the Einstein relation *D*_{n}/*u*_{n}=*D*_{p}/*u*_{p}=*k*_{B}*T*, where *k*_{B} is the Boltzmann constant and *T* the absolute temperature. The length scale *λ*_{*} referred to above is defined by the expression in terms of the electronic charge (*e*), permittivity of vacuum (*ϵ*_{0}) 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 *n* as .

Equations (1.1)–(1.3) and their modifications also describe other transport problems of interest in various areas of physics. For example, if , where *E*_{0} is an external field and is the perturbation in the potential, and if a term −*N*(*x*) representing the density of free charges is added to the right-hand side of equation (1.3), then these equations become the Van Roosbroeck equations describing the migration of charge carriers in solid state physics, where *n*, *p* and *N* represent the concentrations of electrons, holes and dopants (Markowich *et al.* 2002). If *N*(*x*) 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 (Rubinstein 1990). In plasma physics (Krall & Trivelpiece 1973), 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 equations (1.1)–(1.3). The generalization of these equations to three or more species in the presence of an applied electric field *E*_{0} describes the separation of charged macromolecules in capillary electrophoresis (Ghosal 2006) and other separation processes based on electric charge.

One of the simplest experimental realizations of equations (1.1)–(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 *t*=0, 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 (1890) 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 were 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 ∂_{xx}*ϕ* on the left-hand side of equation (1.3), so that it reduces to *p*=−*n*/*z*. This relation can now be used to eliminate the terms involving *ϕ* in equations (1.1) and (1.2). It then follows that *n* (as well as *p*=−*n*/*z*) obeys the diffusion equation
1.4
where *D*=*u*(1−*z*)/(1−*uz*), a result owing to Henderson (1907, 1908). Since *z*<0, clearly *u*≥*D*≥1. 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 equation (1.3), it must be retained in equations (1.1) and (1.2), a situation that appeared self-contradictory and led to some controversy until Hickman (1970) (see also Jackson 1974) provided a proper interpretation within the framework of an asymptotic theory based on expansion in the small parameter *kλ*_{D}, where *k* is a characteristic wave number and *λ*_{D} is the Debye length. Since the Debye length *λ*_{D} is normally a very small quantity (about 3 nm in a 0.01 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 nanometres. 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 equations (1.1)–(1.3).

In this paper, we use the method of multiple scales (Nayfeh 2000) to reduce equations (1.1)–(1.3) to a nonlinear one-dimensional system in the limit of long (compared with the Debye length) length scales and slow (compared with 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 equation for the concentration *n* emerges. Numerical solutions of this reduced equation are seen to agree better with those of the full system of equations (1.1)–(1.3), particularly at early times.

## 2. Asymptotic theory

We would like to study the behaviour of equations (1.1)–(1.3) at long length and time scales, under the boundary conditions that *n*,*p*,*ϕ* 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 (Nayfeh 2000), we introduce slow variables and *τ*=*ϵt* and suppose that *n*,*p*,*ϕ* depend solely on *ξ* and *τ*. Then equations (1.1)–(1.3) reduce to
2.1
2.2and
2.3
where *ϵ* is a small parameter. Expanding all dependent variables in asymptotic series in *ϵ*, such as *n*=*n*_{0}+*ϵn*_{1}+*ϵ*^{2}*n*_{2}+⋯ , substituting in equations (2.1)–(2.3) and equating the coefficients of *ϵ* on both sides, we get a series of equations starting with the lowest order ones
2.4
2.5and
2.6
If the last equation is used to eliminate the second terms from the first two, we get
2.7
and
2.8
with *D*=*u*(1−*z*)/(1−*uz*) representing ambipolar diffusion; if the time derivative terms are eliminated instead and the result integrated, we get
2.9
If in equation (2.9) we put , we get the formula of Planck (1890) for the potential drop across a liquid junction.

At the next order in *ϵ*, we have
2.10
2.11and
2.12
If we add equations (2.10) and (2.11) and use equation (2.6), the term in *ϕ*_{1} is eliminated. We then obtain after substituting *ϕ*_{0} from equation (2.9) and after some algebra
2.13
and
2.14
where *α* and *β* are positive constants defined by
2.15
and
2.16
Equations (2.7) and (2.13) may be combined into a single equation written in terms of our original variables *x* and *t*,
2.17
with an error that is higher than order *ϵ*^{2}. In arriving at equation (2.17), we have replaced the terms by and by , which is justified because doing so introduces an error in equation (2.17) that is higher than order *ϵ*^{2}. Similarly, combining equation (2.14) with equation (2.8), we get
2.18
with an error of order *ϵ*^{2}. The last term in equation (2.18) is due to the departure from charge neutrality.

## 3. Discussion

Let us now consider some of the consequences of equations (2.17) and (2.18). First, we consider weak perturbations from the constant values at infinity: and , where . This gives a hyper-diffusion equation for ,
3.1
Therefore, solutions of the form have growth rates *a*=−*k*^{2}*D*+*k*^{4}*α*. Since *α*>0, high-wavenumber modes are unstable. But equation (2.17) is not valid for modes as such modes violate the premise *k*≪1 on which the asymptotic theory is based. The instability could be either physical or an artefact of the approximations employed. Actually the latter is true, as can be seen by considering the linearized version of equations (1.1)–(1.3) with *u*^{−1}=0,
3.2
3.3and
3.4
If equation (3.3) is integrated and the result substituted in equation (3.4), an exact expression for *ϕ* may be found in terms of
3.5
3.6
3.7
where is Green’s function of the Helmholtz operator ∂_{xx}−*m*^{2}. If the exact inversion, that is, equation (3.6), is substituted in equation (3.2), one gets an integro-differential equation that in Fourier space yields a growth rate
3.8
The exact formula for *a* indicated by the = sign above shows that there is no high-wavenumber instability if the exact inversion, equation (3.6), is employed. The second part, indicated by the ∼ sign, shows the result of the approximate inversion, equation (3.7), 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 equation (3.8).

In the next section, we follow the time evolution of ion concentrations by numerically computing solutions to equation (2.17). In order to do this, the spurious high-wavenumber instability must be suppressed. There are a number of ways of doing this; we adopt the following approach. In the term in equation (2.17), we replace *α*∂_{xx} by (1−*α*∂_{xx})^{−1}−1, since in the long length scale limit we can write ∂_{xx}[1−*α*∂_{xx}]^{−1}∼∂_{xx}(1+*α*∂_{xx}+⋯ )=∂_{xx}+*α*∂_{xxxx}+⋯ . Thus, in place of equation (2.17) we have the nonlinear integro-differential equation
3.9
Equation (3.9) is asymptotically equivalent to equation (2.17) since they differ only by terms of order higher than *ϵ*^{2}, but it is free from the spurious high-wavenumber instability. By virtue of its construction, it also has the property that, when *u*^{−1}=0, the linearized version of equation (3.9) is the exact solution of equations (3.2)–(3.4). Thus, when *u* is large and amplitudes are small, the solutions of equation (3.9) 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. Equation (3.9) is also easy to simulate numerically; the integral is replaced by a numerical quadrature formula, and this is done in the next section.

## 4. Numerical solutions

In order to test the accuracy of the asymptotic equations, we solved a test problem at three levels of approximation:

— (A) using the ambipolar diffusion equation: equation (1.4);

— (B) using the nonlinear evolution equation: equation (3.9);

— (C) using the full Poisson–Nernst–Planck (PNP) equations: equations (1.1)–(1.3).

As a test problem, an initial profile was specified and *n*(*x*,*t*) was set to its initial values and held fixed at the left and right boundaries. The boundaries were set sufficiently far so as to have no practical effect on the solution. The parameters were fixed at *A*=1.0, *u*=6 and *z*=−1.

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 the relevant equations in each of cases A–C. The trapezoidal scheme was used to discretize the integral on the right-hand side of equation (3.9).

It is clear that, in case A, the variance *σ*^{2} would increase linearly with time, but not so in case B. Further, a distribution that is initially Gaussian remains so at future times under case A, but not under case B. Therefore, the excess kurtosis *γ* is expected to remain zero at all times in case A, but not in case B. Thus, the variance *σ*^{2} and excess kurtosis *γ* provide a convenient way to measure departure from the ambipolar diffusion limit characterized by equation (3.9).

In figure 1, the quantities *Δ*≡(2*D*)^{−1} *dσ*^{2}/*dt* and *γ* are plotted as a function of time (*t*) for cases B and C. In case A, as expected, *Δ* and *γ* remain constant in time. It is seen that, at sufficiently long times, solutions for cases A, B, and C all approach a common asymptotic limit. However, at shorter times, the solution for case B is in better accord with the full PNP equations represented by case C.

The qualitative nature of the time variation of *γ* and *Δ* may also be understood on the basis of equation (2.17). To see this, we restrict ourselves to small amplitudes, , with . Then, , so that equation (2.17) becomes
4.1
where the ⋯ indicates nonlinear terms of higher order. If we multiply equation (4.1) by *x*^{k} and integrate, we can generate ordinary differential equations for the moments , starting with *dm*_{0}/*dt*=*dm*_{1}/*dt*=0. The neglected sixth and higher derivative terms do not contribute for *k*≤4. Without loss of generality, we can assume *m*_{1}=0 so that here *σ*^{2}=*m*_{2}/*m*_{0} and . By combining the equations for the second and fourth moments, we derive
4.2
and
4.3
where ⋯ indicates higher order terms. Now, if *α*=0, then either *γ*=0 (if it is zero initially) or else *γ* monotonically decays to zero. When deviations from the ambipolar diffusion limit is considered (*α*>0), 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, equation (4.2) shows that, since *β*>0, the rate of increase of the variance is slightly less than 2*D*, but the deficit in the growth rate becomes progressively smaller at large times. This is indeed what is observed in figure 1.

## 5. Conclusion

In summary, the one-dimensional electro-diffusion problem under a self-consistent electrostatic field was considered in the limit where the Debye length was non-zero, but nevertheless small compared with a characteristic scale of spatial variation. It was shown that in this limit the concentration can be described by a nonlinear equation that reduces to the linear diffusion equation describing ambipolar diffusion if all but the leading order terms in the ratio of the Debye length to a characteristic spatial scale are neglected. The modified equation exhibits a linear instability at short length scales, but this is an artefact of the asymptotic approximation and may be removed without affecting the asymptotic validity of the equation. Numerical solution shows that the nonlinear equation provides a better approximation to the true solution than the linear diffusion equation with an effective diffusivity. The mathematical properties of nonlinear equations describing diffusive transport have been studied in a general setting by various authors (e.g. Andreu *et al.* 2005; Bellomo & Bellouquid 2006).

The formal condition for the validity of equation (3.9) is , where the bars on top of the variables indicate that they are in a dimensional form. In dimensionless notation, the condition of validity reduces to . It should be noted that equation (2.17) as well as the lower order effective diffusivity approximation fail 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 with the local value of the Debye length.

As figure 1 shows, for dimensionless times longer than a characteristic value of order unity, the solutions of equations (2.17) and (1.4) are essentially identical. This is observed for other initial concentration profiles as well, including the ‘step profile’ corresponding to the liquid junction problem though only the Gaussian case is shown here. In physical terms, a dimensionless time of order unity corresponds to a physical time μs with typical parameter values. Thus, the practical usefulness of equation (2.17) is generally quite limited, except perhaps as a way of estimating the size of the neglected terms in equation (2.7).

The range of applicability, however, becomes much broader if one considers electro-diffusion in more than one dimension, since this introduces the possibility of geometric boundaries characterized by some extrinsic length scale, *a*. If *a*∼*λ*_{D}, the ambipolar diffusion limit would be inadequate. One example of such a situation is diffusion of ions through nanometre-sized pores in biological or synthetic membranes (Schoch *et al.* 2008). Biological phenomena that are modelled using the twin assumptions of diffusive transport and local charge neutrality, such as the celebrated Goldman–Hodgkin–Katz voltage equation (Hille 2001), in membrane physiology could potentially be refined following the approach presented here. Another example involves the phenomenon of diffusiophoresis (Abécassis *et al.* 2008). In diffusiophoresis, a non-uniform salt concentration results in the establishment of intrinsic electric fields within the electrolyte owing to the ‘liquid junction’ effect, and this alters the diffusivity of charged colloidal particles. Thus, the diffusivity of colloidal particles depends not only on the gradient of their own concentration, but also on the gradient of concentration of the background electrolyte. The effect has been shown by Ibele *et al.* (2008) to produce nanoscale patterns by self-assembly reminiscent of self-organization in biological systems. The present analysis could be of value as a first step in the study of such multi-dimensional electro-diffusion problems with nanoscale features of interest in the biological sciences.

## Acknowledgements

Support from the National Institutes of Health (NIH) under grant R01EB007596 is gratefully acknowledged.

## Footnotes

- Received September 24, 2009.
- Accepted January 26, 2010.

- © 2010 The Royal Society