## Abstract

A Painlevé II model derived out of the classical Nernst–Planck system is applied in the context of boundary value problems that describe the electric field distribution in a region *x*>0 occupied by an electrolyte. For privileged flux ratios of the ion concentrations, the auto-Bäcklund transformation admitted by the Painlevé II equation may be applied iteratively to construct exact solutions to classes of physically relevant boundary value problems. These representations involve, in turn, either Yablonskii–Vorob’ev polynomials or classical Airy functions. The requirement that the electric field distribution and ion concentrations in these representations be non-singular imposes constraints on the physical parameters. These are investigated in detail along with asymptotic properties.

## 1. Introduction

The physical processes involved in the migration of material particles across boundaries such as permeable membranes have fundamental importance in diverse areas such as electrochemistry, semiconductor theory and biology. The classical Nernst–Planck system (Nernst 1888; Planck 1890) represents a basic model descriptive of these processes involving the transport of charged particles through material barriers. An analysis of electrochemical phenomena at interfaces as modelled by the Nernst–Planck system coupled with Poisson’s equation for the electric field was undertaken by Bass (1964*a*,*b*). Detailed discussion of the appropriate boundary conditions imposed by physical considerations was contained therein.

The role of the Nernst–Planck system in the modelling of ionic transport through biological membranes has been described by Bass (1971) and Bass & Moore (1967), and subsequently in monographs by Cole (1968) and Hille (2001). The problem of the prediction of ionic currents in biological membranes as based on the Poisson–Nernst–Planck (PNP) model remains a subject of much current research interest (q.v. Singer & Norbury 2009 and references therein). In particular, in the absence of exact results, both singular perturbation methods (Liu 2005; Singer & Norbury 2009) and numerical analysis treatments (Kato 1995; Horno *et al.* 1990) have received much attention.

The PNP system has also been established as a basic model in semiconductor theory (Markovich 1984; Markovich & Ringhofer 1984; Markovich *et al.* 1990). Thus, importantly, methods developed in other ionic channel contexts may be portable *mutatis mutandis* to the analysis of semiconductor devices.

The major impediment to our analytic understanding of the Nernst–Planck system as a predictive model in the earlier-mentioned electrochemical, biological or semiconductor contexts remains its underlying nonlinearity. This renders it, in general, analytically intractable. Nonetheless, it turns out that, remarkably, analytic progress in terms of integrability can indeed be made for certain valency classes and ‘quantized’ sequences of ionic flux ratios. Such ‘selective’ situations will be the concern of the present paper with particular attention to electrolytic boundary value problems on *x*≥0. It is noted that selectivity, namely the process that allows the migration of certain types of ions across channels but precludes that of other types of ionic species, is an important observed characteristic of electrolytic permeation processes (Hille 2001).

It was Bass (1964*a*) who, in a study of electrolytic phenomena in the presence of a planar boundary, showed that the Nernst–Planck system, augmented by Poisson’s equation, reduces in a particular two-ion case to consideration of the classical Painlevé (1902) II equation. In his subsequent analysis, Bass relied upon an approximate solution to the Painlevé II equation in which the nonlinear term was neglected, thereby precluding a complete physical description of the electrolytic process. However, at that time, the mechanism to generate sequences of exact solutions of the Painlevé II equation was not known. The procedure, in fact, involves the application of a Bäcklund transformation originally obtained by Lukashevich (1971) and subsequently rediscovered in the context of modern soliton theory by Boiti & Pempinelli (1979) and Airault (1979). It is recalled that Bäcklund transformations have extensive applications in both nonlinear continuum mechanics and soliton theory (q.v. Rogers & Shadwick 1982; Rogers & Schief 2002). In Conte *et al.* (2007), Painlevé analysis was systematically applied to the general steady-state *m*-ion coupled electrodiffusion system as set down by Leuchtag (1981). The physical parameters for which the system admits a general solution that is single-valued were thereby delimited. It was established, remarkably, that the Painlevé II equation recurs in a canonical manner in *m*-ion reductions with *m*>2 for selective valency ratio classes.

Rogers *et al.* (1999) returned to the boundary value problem originally posed and treated approximately by Bass (1964*a*), which describes the electric field distribution in a region *x*>0 occupied by an electrolyte. The auto-Bäcklund transformation admitted by the Painlevé II equation was applied iteratively to construct exact representations for the electric field distribution for boundary value problems, wherein the ratio of fluxes of the positive and negative ions adopts one of an infinite sequence of values. Here, the singularity structure of these classes of Bäcklund-generated exact representations is analysed in order to ascertain the range of their physical application. These hierarchies of exact solutions of the Painlevé II equation involve the logarithmic derivative of either Yablonskii–Vorob’ev polynomials (Yablonskii 1959; Vorob’ev 1965) or classical Airy functions so that physical considerations require that parameter values be restricted in such a way as to avoid their zeros.

It turns out that the distribution of complex zeros of the Yablonskii–Vorob’ev polynomials has a remarkably regular structure and their study in connection with integrable Painlevé hierarchies is of current research interest. Thus, Clarkson & Mansfield (2003) have exhibited highly symmetric ‘triangular’ confinement regions in the complex plane for the standard Yablonskii–Vorob’ev polynomials. This elegant work has most recently been followed by consideration of generalized Yablonskii–Vorob’ev polynomials and the Painlevé II hierarchy by Kudryashov & Demina (2007, 2008).

Here, our concern is with the distribution of real zeros of Yablonskii–Vorob’ev polynomials and the application of classical results for the Airy function to isolate conditions for the regular solution of the two-ion boundary value problem as originally posed in Bass (1964*a*).

## 2. The Poisson–Planck–Nernst model

A steady-state electrolytic boundary value problem is considered wherein a strong, dilute, symmetric -valent electrolyte occupies the region *x*>0 bounded by the electrode interface *x*=0. The macroscopic electric field *E*(*x*) exerts a force on each of the positive or negative ions, where *e* denotes the elementary charge.

The positive and negative ions are supposed to have concentrations *c*_{+}(*x*) and *c*_{−}(*x*), respectively, giving rise to the electric field ** E**=

*E*

**via Poisson’s equation 2.1 where 2.2 and 2.3 with**

*i**ε*the dielectric constant, so that 2.4 In the absence of volume sources of ions in strong electrolytes, the numbers of ions of each sort are conserved in each volume element, so that, in the steady state to be examined here, 2.5 for

*x*>0, where

*D*

_{±}are the diffusivity coefficients, here taken as constants, while

*u*

_{±}denote the ionic mobilities. In addition, Einstein’s relation, valid for each type of ion, yields 2.6 where

*k*

_{B}is Boltzmann’s constant and

*T*is the absolute temperature.

In the plane geometry under consideration, the steady-state conservation laws (2.5) on use of the Einstein relations give, on integration,
2.7
and
2.8
It is noted that the constants of integration *A*_{±} are such that
2.9
where *Φ*_{+} and *Φ*_{−} are the fluxes of the ion concentrations *c*_{+} and *c*_{−}.

The general solution of the system under consideration, namely equations (2.7) and (2.8) augmented by Poisson’s equation (2.4), contains, in addition to the constant weighted fluxes *A*_{±}, three constants of integration. For the present boundary value problem as posed on *x*>0, the pertinent boundary conditions at the electrode surface *x*=0 involve
2.10
It is noted that the latter pair of boundary conditions are equivalent to specification of the boundary concentrations
2.11

In the following, it will be shown that exact solutions of the nonlinear coupled electrolytic system consisting of equations (2.7) and (2.8) together with equation (2.4) may be constructed for any one of an infinite sequence of privileged ‘weighted’ flux ratios
2.12
The corresponding electric field is non-singular while the concentrations are positive and non-singular in accordance with physical requirements. The electric field
2.13
on the boundary may be prescribed arbitrarily while the total edge concentration
2.14
may be chosen within a set *S*, which is determined by *E*_{0} and *R*.

## 3. The Painlevé II reduction

Here, we review the reduction of the electrolytic system to the second Painlevé equation *P*_{II}. This was originally obtained in Bass (1964*a*).

On addition of Nernst–Planck equations (2.7) and (2.8) and the use of Poisson’s equation (2.4), it is seen that
3.1
integration of which yields
3.2
where
3.3
denote the boundary data on the electrode. Subtraction of equations (2.7) and (2.8) gives
3.4
where equation (2.4) shows that
3.5
whence
3.6
Elimination of the total concentration *n*=*c*_{+}+*c*_{−} in equation (3.6) via relation (3.2) now yields, as in Bass (1964*a*),
3.7
With boundary data *E*_{0} and *c*_{±}(0), the edge concentration *n*_{0}=*c*_{+}(0)+*c*_{−}(0) is known, while the boundary electric field gradient *E*′_{0} is given by Poisson’s equation. Alternatively, if *E*_{0},*E*′_{0} and *n*_{0} are prescribed, then the separate ion concentrations on the boundary are given via
3.8

Introduction of the change of variables
3.9
where
3.10
into equation (3.7) now leads to the standard Painlevé (1902) II equation
3.11
with
3.12
Here, we are concerned with the generic case *A*_{+}+*A*_{−}≠0. The ion concentrations are then given by
3.13
with
3.14
and
3.15
Thus, if, in accordance with equation (3.9), we constrain the two constants of integration in *Y* (*z*) by
3.16
then the general solution of the original electrolytic system is obtained via the Painlevé reduction with five arbitrary parameters *E*_{0},*n*_{0},*A*_{±} and *Y* ′(*γ*) as required. Moreover, since the sign of *β* coincides with that of *A*_{+}+*A*_{−}, one is concerned with solutions on the half-space
or
In this connection, it is observed that the relation (3.14) shows that
3.17
whence, on the one hand, the total concentration cannot be bounded on the half-space *z*>*γ* and, on the other hand, the electric field and the total concentration cannot be bounded simultaneously on the half-space *z*<*γ*. This reflects the well-documented physical problems inherent in making the transition to as discussed by Bass (1964*a*). Thus, a physically realistic model necessarily involves consideration of a bounded region (Bass 1964*c*). This aspect will be discussed subsequently. It is emphasized that the results to be presented here on the application of the Bäcklund-generated solutions of the Painlevé II equation in the electrolytic context are also, in principle, pertinent to boundary value problems posed on bounded regions such as those encountered in bio-membrane theory.

In what follows, it proves natural to prescribe the data
3.18
so that the ‘total weighted flux’ *A*_{+}+*A*_{−} is given by
3.19
while the total concentration *n*_{0} on the boundary is determined parametrically in terms of *γ* by
3.20

## 4. Iterative application of a Bäcklund transformation

It turns out, remarkably, that all known exact solutions of the Painlevé II equation may be generated by an elegant auto-Bäcklund transformation. This result was originally obtained by Lukashevich (1971) and is encapsulated in the following theorem.

## Theorem 4.1

*If Y (z;C) is a solution of P*_{II} *with parameter C, then*
4.1
*is a solution of P*_{II} *with* *, C≠−1/2.*

In the present electrolytic context, this produces a nonlinear superposition principle involving the physical variables *E*,*c*_{+} and *c*_{−} embodied in (q.v. Rogers *et al.* 1999) the following corollary.

## Corollary 4.2

*The electric field* *associated with the solution* *is given by*
4.2
*where*
4.3

Iterative application of the earlier-mentioned Bäcklund transformation to appropriate seed solutions of *P*_{II} leads to an infinite sequence of exact solutions expressed as logarithmic derivatives of either Yablonskii–Vorob’ev polynomials or Airy functions. These sequences, which may be encoded in quadratic recurrence relations, are now analysed in the present electrolysis context. To this end, it is convenient to exploit the symmetry
4.4
of the governing system to restrict attention to the case
4.5
without loss of generality.

### (a) The rational solutions

In Rogers *et al.* (1999), the Bäcklund transformation was applied iteratively to the base solution corresponding to the seed electrically neutral state
4.6
associated with *C*=0 to obtain a sequence of rational solutions. In general, if
4.7
then iterative application of the Bäcklund transformation may be used to establish that the scaled electric field
4.8
constitutes a particular solution of *P*_{II}, wherein the functions *P*_{k} are defined in terms of the quadratic recurrence relation
4.9
Remarkably, the *P*_{k} are monic polynomials of degree (*k*+1)(*k*+2)/2 with integer coefficients and each term in *P*_{k} has the same degree modulo 3 (Vorob’ev 1965; Clarkson 2003; Noumi 2004). In particular, the Yablonskii–Vorob’ev polynomials *P*_{0},…,*P*_{3} read
4.10
It may be shown that the associated scaled concentrations *p*_{+}[*k*] and *p*_{−}[*k*] are given by
4.11
where the rational functions *p*[*k*] are defined by
4.12
The above generate rational expressions for the scaled electric field and ion concentrations as follows:
4.13
It is noted that rational solutions of the ordinary differential equation (*P*_{34}) satisfied by *p*_{±} have been discussed in Clarkson (2003).

In order to address the physical significance of the above rational solutions, it is necessary to investigate their singularity structure. To this end, it is observed that the degree of the monic polynomials *P*_{k} implies that
4.14
Since the concentrations *c*_{±} are required to be positive, the above solutions are only applicable on the half-space
4.15
so that
4.16
Thus, if *z*_{k} denotes the largest real zero of the polynomials *P*_{k},*P*_{k−1},*P*_{k−2} and *P*_{k−3} for each , that is
4.17
(with *P*_{−3}:=1), then any
4.18
guarantees that the scaled electric field *Y* [*k*] is non-singular and that the scaled concentrations *p*_{±}[*k*] are non-singular and positive on the half-space *z*>*γ*. This is indicated in figure 1 for *k*=0,1,2.

The preceding analysis gives rise to the following construction of non-singular rational solutions of the original electrolytic system. For any ‘selective’ discrete value
4.19
of the ratio of weighted fluxes, the scaled solution (*Y* [*k*],*p*_{±}[*k*]) is known. The modulus of the electric field *E*_{0} on the boundary *x*=0 may be chosen arbitrarily and its sign is that of *Y* [*k*](*γ*) for any fixed *γ*>*z*_{k} by virtue of relation (3.20), that is
4.20
The corresponding physical solution
4.21
is then unique with
4.22
so that, indeed, *E*_{0}=*αY* [*k*](*γ*) (cf. equation (3.16)). The total concentration *n*_{0} at the boundary is given by
4.23
and may be adjusted by varying the parameter *γ*. Since *Y* [*k*] and *p*_{+}[*k*]−*p*_{−}[*k*] are bounded, the electric field *E* and its gradient *E*′ are bounded and, in fact, vanish at infinity, while the total concentration *c*_{+}+*c*_{−} grows linearly as by virtue of equation (4.14). The class of exact solutions so constructed contains one discrete and two continuous parameters, namely *A*_{−}/*A*_{+},*E*_{0} and *γ*. The precise connection with the natural physical boundary value problem that requires the prescription of the boundary values *E*_{0},*E*_{0}′ and *n*_{0}, that is, the involved nature of the mapping
4.24
remains a matter of investigation.

### (b) Airy-type solutions

If *C*=1/2, then *P*_{II} admits the solution
4.25
where *ϕ* is any solution of the (scaled) Airy equation
4.26
A chain of exact solutions of *P*_{II} for half-integers
4.27
may then be generated via the Bäcklund transformation (4.1) with this seed. Indeed, on iteration of the Bäcklund transformation, it may be shown that *P*_{II} admits an infinite sequence of exact solutions
4.28
where the sequence of functions (*ϕ*_{l})_{l≥−1} is defined by the recurrence relation
4.29
with initial values
4.30
The *ϕ*_{k} turn out to be homogeneous polynomials of degree *k*+1 in *ϕ* and *ϕ*′ with
4.31
where Ai and Bi are the Airy functions of the first and second kind, respectively, and *a*,*b* are arbitrary constants of integration.

It proves convenient to express the above solutions in terms of the related quantities
4.32
and
4.33
The solutions of *P*_{II} are then given by
4.34
wherein the sequence (*ψ*_{l})_{l≥−1} is determined by the recurrence relation
4.35
together with the initial values *ψ*_{−1}=1 and *ψ*_{0}=1. The first three non-trivial terms read
4.36
The corresponding scaled concentrations may be conveniently expressed in terms of the same sequence (*ψ*_{l}) according to
4.37
where
4.38
On use of the recurrence relation (4.35), we can write the latter expression as
4.39
and thereby obtain the Airy-type solutions
4.40

The general solution of Airy’s equation in standard form has infinitely many negative zeros, whence we deduce that *Φ*(*z*) is non-singular on the half-space
4.41
where *z*_{0} is the first positive zero of *ϕ*(*z*). This implies that the Airy-type solutions of *P*_{II} as generated by the Bäcklund transformation can only have physical application if
4.42
Furthermore, for the general solution *ϕ* of Airy’s equation, we have the asymptotic results
4.43
as . From this, we conclude that
4.44
and then, using equation (4.39), that
4.45
Thus, in the case *b*≠0, the scaled concentrations are asymptotically negative and hence unphysical. Accordingly, we set *b*=0 and proceed only with the case
4.46
without loss of generality.

The asymptotic formulae (4.43) and (4.45) together with the asymptotic property
4.47
as , show that, for any *k*, there exists a *z*_{k} such that the scaled electric field is non-singular and the scaled concentrations are non-singular and positive on *z*<*z*_{k}. The Airy-type solution behaviour is illustrated in figure 2 for *k*=1,2,3. The corresponding physical variables
4.48
may now be reconstructed in the manner described in the previous section. However, here,
4.49
and *E*_{0} and *Y* (*γ*) are required to have opposite signs by virtue of constraint (4.42) and relation (3.19).

### (c) Boundedness

In the above case, it is evident that the concentrations *c*_{±} vanish as but that the electric field *E* is unbounded. To ensure that the Airy-type solutions are amenable to a physical interpretation, one could either consider two electrodes at a finite distance from each other or alternatively investigate how the electric field is ultimately joined up with a constant (ohmic) field in an electrically neutral convective region (cf. Bass 1964*a*). It appears that the latter case is natural in the current context. Thus, we assume that the diffusion layer governed by the earlier-mentioned solutions is confined to the region 0<*x*<*ξ* and Ohm’s law in the form
4.50
holds in the convective region , where , and *i* is the current density. Since the current density must be continuous at the (notional) interface *x*=*ξ* between the diffusion layer and the convective region, the expression
4.51
valid throughout the diffusion layer implies that the electric field *E*^{+}=*E*(*ξ*+) just outside the diffusion layer is given by
4.52
This is to be compared with the electric field *E*^{−}=*E*(*ξ*−) in the diffusion layer as the field jump at *x*=*ξ* corresponds to a surface charge density *ω* distributed over the interface and encapsulated in the relation
4.53

In the present situation, relations (3.19) and (4.49) lead to the expression
4.54
in terms of the data (*k*,*E*_{0},*γ*), while
4.55
If, as in Bass (1964*a*), we introduce the Debye length
4.56
and recall that
4.57
then we may introduce a generalized Debye length according to
4.58
Since may be made arbitrarily small for an appropriate choice of *γ*, it is natural to replace the usual assumption of *ξ*≫*κ*^{−1} by
4.59
that is, the distance between the electrode and the interface is assumed to be ‘large’ compared with the generalized Debye length. Relation (4.55)_{2} then implies that since *γ*<*z*_{0}≈2.946 and hence
4.60
in the expression (4.55) for *E*^{−}. In particular, the asymptotic formulae (4.43), (4.45) and (4.47) or direct inspection of the integrated conservation laws (2.7) and (2.8) give rise to
4.61
for so that equation (4.52) becomes
4.62
4.63
by virtue of equation (4.49).

## 5. General considerations

The steady-state Nernst–Planck two-ion system treated in Bass (1964*a*) and in the present work is privileged by virtue of its integrable reduction to the Painlevé II equation. It is natural to investigate whether integrable reductions arise in other multi-ion systems. In this connection, it proves convenient to partition the ions into *m* classes characterized by the same electric charge *q*_{j}=*q*_{0}*v*_{j}, where *q*_{0} is the unit of charge and *v*_{j} is a non-zero signed valency. The *m*-ion electrodiffusion system in steady régimes then reduces to the following system of *m*+1 coupled first order equations (Leuchtag 1981)
5.1
where *x* is the coordinate normal to the planar boundary, *p* is the electric field and *n*_{j} is the number of ions with the same charge *q*_{j}=*q*_{0}*v*_{j}. It is noted that the system (5.1) admits the first integral (cf. equation (3.2))
5.2

In Conte *et al.* (2007), classical tests with origins in work by Kowaleski (1889*a*,*b*), Painlevé (1900) and Gambier (1910) were applied to the nonlinear system (5.1) and integrable reductions were explicitly shown to exist for certain valency classes. The corresponding solutions turn out to be given in terms of either Painlevé II transcendents or Jacobi elliptic functions. In the two-ion case, reductions were obtained for the three ‘selective’ valency ratios
5.3
The latter case leads to the Painlevé II reduction obtained by Bass (1964*a*).

### (a) The two-ion case

It proves instructive to examine the general two-ion case *m*=2 in more detail. Thus, if we set *n*_{1}=*n*_{+},*n*_{2}=*n*_{−},*v*_{1}=*v*_{+},*v*_{2}=*v*_{−},*A*_{1}=*A*_{+} and *A*_{2}=*A*_{−}, then the system (5.1) adopts the form
5.4
It is noted that an analogous system in the case *v*_{+}+*v*_{−}=0 has been set down in the context of semiconductor theory by Kudryashov (1997).

The first integral (5.2) now yields
5.5
and elimination of *n*_{+} between equations (5.4)_{3} and (5.5) then gives the electric field gradient
5.6
Accordingly,
5.7
In the case *v*_{+}+*v*_{−}=0 as considered in Bass (1964*a*) and Kudryashov (1997), reduction to *P*_{II} is evident. However, if *v*_{+}+*v*_{−}≠0, then the electric field equation (5.7) is analytically intractable except for the selective valency classes *v*_{+}:*v*_{−}=1:−2 and *v*_{+}:*v*_{−}=1:2 as isolated in Conte *et al.* (2007).

### (b) Boundary conditions

Even in the case of the integrable Painlevé II reduction when *v*_{+}+*v*_{−}=0, the side conditions imposed on bounded regions produce an additional impediment to analytic progress. It is remarked that such two-point boundary value problems arise naturally in the study of solubility processes in neuronal membranes (Bass & Moore 1967; Thompson 1994). In that context, if the membrane occupies the region 0≤*x*≤*ξ* with cell interior *x*<0 and exterior *x*>*ξ* occupied by appropriate physiological solutions, then the boundary data are commonly taken as the boundary concentrations together with the net current density. If, as in Thompson (1994), there is no current in the junction, then
5.8
where *D*_{±} denotes the respective diffusivity coefficient, and hence
5.9
Accordingly, the electric field equation (5.7) becomes
5.10
where equation (5.5) shows that
5.11
and
5.12
Thus, it is seen that equation (5.11) incorporates via *A*_{+}+*A*_{−} and *K* the total edge concentrations *n*(0)=*n*_{+}(0)+*n*_{−}(0) and *n*(*ξ*)=*n*_{+}(*ξ*)+*n*_{−}(*ξ*), together with the boundary values of the electric field *p*(0) and *p*(*ξ*). Whereas it is anticipated that the interface concentrations *n*_{+}(0),*n*_{−}(0),*n*_{+}(*ξ*) and *n*_{−}(*ξ*) be known (Bass & Moore 1967), the solution-dependent boundary terms *p*(0) and *p*(*ξ*) remain. The intrusion of the latter boundary data poses a formidable obstacle to analytic progress. In recent work, Amster *et al.* (submitted) addressed this problem in the case when electrical neutrality in the reservoirs imposes the Neumann boundary conditions
5.13
A novel two-dimensional exact shooting approach was adopted and a practical algorithm for the numerical solution of the nonlinear two-point Neumann boundary value problem was presented.

## Footnotes

- Received November 23, 2009.
- Accepted January 14, 2010.

- © 2010 The Royal Society