IOP PUBLISHING

JOURNAL OF PHYSICS A: MATHEMATICAL AND THEORETICAL

doi:10.1088/1751-8113/45/45/455204

J. Phys. A: Math. Theor. 45 (2012) 455204 (13pp)

Unsteady analytical solutions to the
Poisson–Nernst–Planck equations
Johannes Schönke
Interdisciplinary Center for Scientific Computing, University of Heidelberg, INF 368,
D-69120 Heidelberg, Germany
E-mail: johannes.schoenke@iwr.uni-heidelberg.deReceived 22 March 2012, in final form 26 September 2012
Published 29 October 2012
Online at stacks.iop.org/JPhysA/45/455204
Abstract
It is shown that the Poisson–Nernst–Planck equations for a single ion species
can be formulated as one equation in terms of the electric field. This previously
not analyzed equation shows similarities to the vector Burgers equation and is
identical with it in the one dimensional case. Several unsteady exact solutions
for one and multidimensional cases are presented. Besides new mathematicalinsights which these first known unsteady solutions give, they can serve as test
cases in computer simulations to analyze numerical algorithms and to verify
code.
PACS numbers: 02.30.Jr, 82.45.Gj(Some figures may appear in colour only in the online journal)

1. IntroductionMore than 120 years ago, Nernst [1] and Planck [2] found the proper mathematical description
for a process which we call electrodiffusion, but the mathematical analysis of the Poisson–
Nernst–Planck (PNP) equations is still far from complete. In full 3D space only the steadyequations and the long time behaviour of unsteady solutions have been examined in terms ofexistence, stability and uniqueness, see e.g. [3]. For the 1D case some deeper insights into
the structure of solutions were gained. Especially the steady two species system was studied
qualitatively in [4] and static solutions for such a system with insulating (zero flux) boundary
conditions have been derived in [5]. Further, unsteady solutions for the linearized equations
were found [6]. To our knowledge, there is no publication known to give explicit unsteadysolutions to the full (nonlinear) PNP system.
In this paper the analysis of the single species PNP system and the resulting unsteadysolutions serve mainly two purposes: first, to contribute to the mathematical analysis byshowing the connections between the PNP system and an equation similar (in 1D identical)
to Burgers equation; second, to present a set of explicit unsteady solutions which can be
1751-8113/12/455204+13$33.00 © 2012 IOP Publishing Ltd

Printed in the UK & the USA

1

J. Phys. A: Math. Theor. 45 (2012) 455204

J Schönke

used in a test suite for the analysis and verification of numerical codes. The second point
seems especially important if one looks at the (currently inevitable) procedure for unsteadynumerical tests, where the original PNP equations need to be modified to result in a previously
constructed, artificial solution, see e.g. [7].
The rest of the paper is organized as follows. In section 2 we give a short introduction to
the PNP equations. In section 3 we perform the analysis of the single species system, whereas
the boundary conditions are discussed in section 4. The 1D case is investigated in section 5,
where also several exemplary solutions are given. The construction of 2D/3D solutions out of
1D solutions is discussed in section 6. In section 7 some comments on the spherical symmetriccase can be found. Section 8 gives the conclusions.
2. The PNP model
Ions drifting through a static solvent can be described by the PNP system, consisting of the
Nernst–Planck equation (stating the conservation of mass)∂ni+ ∇ · Fi = 0(1)∂twith the ion flux
Fi = −Di (∇ni + zi ni ∇φ) ,(2)

and the Poisson equation∇ · (ϵ∇φ) = −e2 n∗ !
zi ni .
ϵ0 kT i(3)

ni is the relative concentration (with respect to a scaling concentration n∗ ), zi is the valence andDi the (possibly position dependent tensorial) diffusion coefficient of the ion species i. φ is
the relative electric potential energy with respect to the thermal energy (φ = eU/kT with U in
Volts, for room temperature φ = 1 corresponds to about 25 mV). ϵ is the relative permittivity(which again may be position dependent and tensorial) and T is the temperature of the solvent.
Equations (1) and (3) constitute a nonlinear system of i + 1 partial differential equations.
Equation (3) is of elliptic type, while (1) by itself represents a system of i uncoupled advection–diffusion equations with the variable advective velocity field vi = −Di zi ∇φ.
3. The single species system
If we only consider a single ion species (with concentration n and valence z) and further
assume the diffusion coefficient as well as the relative permittivity to be spatially independentand isotropic, we can write the system (1) and (3) as∂n= D ∇ · (∇n + zn∇φ)
(4)
∂t
zn
$φ = − 2(5)λwith the characteristic length scale"
ϵϵ0 kTλ=.
e2 n∗For a monovalent species(|z| = 1) the right-hand side of (5) is directly connected to the Debye#
length λ2Deb = λ2 / z2i ni , therefore we have $φ = −1/(zλ2Deb ).

2

J. Phys. A: Math. Theor. 45 (2012) 455204

J SchönkeWith the scaling of length and time coordinates through the introduction of dimensionlessvariables x = r/λ (where r is the original space variable) and τ = tD/λ2 in (4) and (5) we get∂n= ∇ · (∇n + zn∇φ)(6)∂τ$φ = −zn .

(7)∇ · E = zn .

(8)Introducing the electric field E = −∇φ, equation (7) reads

If we use (8) to eliminate n from (6) and use the equality of mixed partials for E (this equality
relies on the existence and continuity of the respective mixed derivatives of E, which translateswith (8) to the existence and continuity of ∂n/∂t) we can integrate (6) once in space and arriveat∂E= ∇ ∇ · E − zE (∇ · E) + ∇ × G,
(9)∂τ
where the last term ∇ × G is a gauge field with an arbitrary and possibly time dependent vector
field G(x, τ ) (resulting from the spatial integration). With the transformations
E′ = −zEandG′ = −zG(10)plus the fact that E is a conservative field by definition (so we can use the fact that∇ ∇ · E = ∇ · ∇E = $E), equation (9) reduces to∂E′= $E′ + E′ (∇ · E′ ) + ∇ × G′ .
(11)
∂τ
′
Even for a vanishing source term ∇ × G , this is unfortunately not the vector Burgers equation,
since E′ (∇ · E′ ) ̸= (E′ · ∇ )E′ . Therefore the vectorial extension of the classical Hopf–
Cole transformation is of no use here. Although no comparable general solution procedure
is outlined in this work, we present several exact solutions to (11) in the following sections.
There were found no further references to equations like (11) in the literature.
4. Boundary conditions
4.1. Nernst–Planck equation
On the boundary ∂' of the considered domain ', a condition for the ion flux (2) may be given
for the Nernst–Planck equation (1) through specification of the (possibly time dependent)
normal component of the ion flux F over the boundaries, i.e. n · F = f (τ ). Translated to our
new variable E′ the condition reads$%n · ∂E′ /∂τ − ∇ × G′ = f (τ ) on ∂'.
(12)
We distinguish two relevant cases. The first one is a zero flux (or closed) system, f = 0. The
second one is an in or outflow condition f ̸= 0 (corresponding to an electric current withregard to a given boundary surface). With respect to the concentration n the flux condition isindeed a Robin boundary condition (i.e. a weighed combination of Dirichlet and Neumannboundary conditions) as can be seen from (2).
4.2. Poisson equation
For the Poisson equation (3), a Dirichlet boundary value has to be prescribed at at least one
point on ∂', since the potential is only defined up to a spatial constant. Apart from that either
Dirichlet or Neumann conditions can be prescribed.
3

J. Phys. A: Math. Theor. 45 (2012) 455204

J Schönke

4.3. Conditions for analytic solutions
Normally, analytic solutions are given on an unbounded domain. If we therefore choose aboundary ∂', the solution itself defines the boundary condition type and its values on ∂'.
Depending on the solution, the construction of different boundary conditions is possible andwill be discussed in the given examples below.
5. Solutions in the 1D caseIn one space dimension (11) reads′∂E ′∂ 2E ′′ ∂E+E=+ g(τ ),
(13)
∂τ
∂x2
∂x
where g(τ ) is the (spatially independent) 1D remnant of the gauge field ∇ × G′ . With thetransformations (cf [8])& τ′g(τ̄ ) dτ̄ =: E ′ − )(τ )
(14)
η=E −
τ0

y=x+

& τ
τ0

(τ − τ̄ ) g(τ̄ ) dτ̄ =: x + γ (τ )(15)

we arrive at the classical (scalar) Burgers equation
∂η
∂ 2η
∂η
= 2 +η
.
(16)
∂τ
∂y
∂y
From the physical point of view, the transformations (14) and (15) describe the action of an
external, spatially constant (but possibly time dependent) electric field ) which moves theions through space (done by the shift γ in the spatial coordinate).
One family of solutions to (16) can be found through the well-known Hopf–Coletransformation
2 ∂ψ
∂
,
(17)
η = 2 (ln ψ ) =∂yψ ∂yresulting in the heat equation
∂ψ
∂ 2ψ
=
.
(18)
∂τ
∂y2
Note that from the definition of E ′ in 1D (E ′ = z ∂φ/∂x) we haveφ(x, τ ) = 1z [2 ln ψ (x + γ , τ ) + x)] + φ0 (τ ),(19)∂η/∂τ = f (τ ) on ∂'.

(20)

with an arbitrary function φ0 (τ ). Finally, the concentration is given by2 ∂2n(x, τ ) = − 2 2 [ln ψ (x + γ , τ )].
z ∂x
The flux boundary condition (12) for the Nernst–Planck equation translates with (14) toGiven one of the several solutions to (18), we can now construct solutions to the Poisson–
Nernst–Planck system.
5.1. Examples in 1D
The given examples below refer to the system (6) and (7), the solutions are therefore given asn(x, τ ) and φ(x, τ ).
4

J. Phys. A: Math. Theor. 45 (2012) 455204

J Schönke

5.1.1. Spatially homogeneous outflow. The simplest possible example is that of a spatiallyhomogeneous evolution of the concentration, i.e. n = n(τ ). Although we could derive eventhis result with the help of the formalism developed in the previous sections, it is easier todeduce from (8) that in this case the electric field must have the form E = znx. From thedefinition of E and (6) then follows the ODE dn/dτ = −z2 n2 which has the solution
1,n(τ ) =
θ +τwith an arbitrary constant θ . The solution describes a symmetric flow away from x = 0. The
respective potential iszz x2
φ(x, τ ) = − n(τ )x2 + φ0 (τ ) = −+ φ0 (τ ) .
(21)
2
2 θ +τ
This example is well suited as a simple numerical test. We start at τ = 0 with a constantconcentration n(0) = 1/θ as the initial condition. Setting one boundary at x = 0, we define
a closed (no flux) boundary there, since the flux is given as F (x, τ ) = z2 n2 (τ )x. At thesecond boundary, say x = x0 , we then just have to set the (time dependent) outflow conditionf (τ ) = z2 n2 (τ )x0 (cf section 4)1. For the potential φ there are two possible boundary condition
settings which produce the same result (up to a free constant): we could either use the valuesof (21) at both boundaries for a Dirichlet condition, or we set φ to an arbitrary value at x0(say, φ(x0 , τ ) = 0) and use the Neumann zero condition at x = 0 (∂φ/∂x|0 = 0) which is an
inherent property of the solution (21).
5.1.2. Closed system. In this more advanced examples we describe the full formalism needed
to construct an analytic solution. At first, we have to choose a solution for (18), in this example
it reads (cf [9] for an extensive list)2ψ (y, τ ) = e−C τ sin(Cy + A) + B,

(22)

where A, B, C are free parameters. Then we can use (17) to find the respective solution η(y, τ )to the Burgers equation (16). Choosing g = 0 in (13) reduces (14) and (15) to E ′ = η andx = y. Using (10) we get the electric fieldcos(Cx + A)
2C
.
(23)
E(x, τ ) = −z BeC2 τ + sin(Cx + A)
With (8), the solution for the concentration then is2n(x, τ ) =2C2 1 + BeC τ sin(Cx + A)'
( ,z2 BeC2 τ + sin(Cx + A) 2

(24)

and the potential is given with (22) and (19) as
2
2
φ(x, τ ) = ln[e−C τ sin(Cx + A) + B] + φ0 (τ ) .
(25)zThis solution is periodic in space with the period P := 2π /C. Further, it has the zero flux
property (cf section 4) at x0 ( j) = (π j − A + π /2)/C , j ∈ Z. Therefore we get a closed(insulated) system, if the boundaries are positioned at two of these zero flux points. The
distance between two consecutive x0 is P/2. For B ! 1 the solution has singularities. For
B > 1 there always exist regions where n < 0 (if a full spatial period is considered), meaningthat the ions behave like oppositely charged particles in those regions. Even if this property
appears strange from a physical perspective, it causes no problems if we want to use this
1 From the point of view of a numerical test it is more convenient here to prescribe Neumann zero conditions for n
at all boundaries.

5

J. Phys. A: Math. Theor. 45 (2012) 455204

J Schönke

0

n and

-1

-2

-3

-4

-5

-5

-4

-3

-2

-1

0

1

2

3

4

5

x
Figure 1. The concentration (24) and the potential (25) for the closed system example. n(x, τ )(solid lines) and φ(x, τ ) (dotted lines) are shown for τ = 0 (red), 1 (blue), 3 (green) and 10(yellow).

example as a numerical test (except numerical schemes which rely on the positivity of n, ofcourse); there is no mathematical reason for n to be strictly positive. Since n can be negative,0= x0 ( j1,2 ) is
it is possible that the total charge Q between two arbitrary zero flux points x1,2
zero. It is the case here, because from (8) and (23) we find Q to be proportional to
& x20
$
%
$
%
z
n(x, τ ) dx = E x20 , τ − E x10 , τ = 0.
x10

As a numerical example, the choice A = −0.5π , B = 1.2, C = 0.2π covers a full period
between the zero flux points x0 (±1) = ±5. If these positions are chosen as boundaries, theexample presents a closed system with vanishing total charge. The concentration and thepotential are depicted in figure 1, where φ0 (τ ) is chosen such that φ(±5, τ ) = 0.
5.1.3. Pure advection. We give another example, where an otherwise stationary solution
is advected through space with constant velocity. If ψ is given as a simple product ofspace and time dependence, η(y) does not depend explicitly on time, cf (17), (but possibly
implicitly through y). From the linear combination of two basic solutions to (18) given by6

J. Phys. A: Math. Theor. 45 (2012) 455204

J Schönkeψ (y, τ ) = exp(A2 τ ± Ay) we can therefore construct a stationary solution η(y) to (16) withthe use of (17)η(y) = 2A tanh(Ay) .

(26)

To give an example for the usage of the gauge field g(τ ), i.e. the inclusion of an external
electric field, we choose the simple case of a delta distribution g(τ ) = −vδ(τ ) with a constantv. From (14) and (15) (with τ0 = 0) we get the external field2 ) = −v, the shift γ = −vτ andtherefore y = x − vτ . The electric field then reads
1
E(x, τ ) = − [2A tanh(A(x − vτ )) − v].
zWith (8), the concentration is
2A2[1 − tanh2 (A(x − vτ ))].
(27)
z2
The concentration profiles are advected through space with the velocity v. The potential isgiven with (19) as1(28)
φ(x, τ ) = [2 ln cosh(A(x − vτ )) − xv] + φ0 (τ ).
z
For x → ±∞ the potential becomes asymptotically linear with the slope (±2A − v)/z.
The concentration profiles are kept ‘in shape’ due to electric forces (proportional to ∇φ)at the domain boundaries, in contrast to the vanishing boundary forces in the previous
example (section 5.1.2) where the profiles flattened out with time. Regarding the negativityof the concentrations (n(x, τ ) < 0 everywhere) we refer to the comments from the previous
section 5.1.2.
In figure 2 the concentration and the potential are shown with A = v = 1 in the interval[−3, 5], where φ0 (τ ) is chosen such that φ(5, τ ) = 0 and the other Dirichlet boundary valueφ(−3, τ ) is then computed according to (28). Boundary conditions for n may be specified
directly as Dirichlet with values from (27), or as flux conditions with f (τ ) calculated from(20) and (26).
n(x, τ ) = −6. Solutions in the 2D/3D case
6.1. Construction from 1D solutions
Looking back at (11), it is clear that for multidimensional problems in Cartesian coordinates
the only coupling between the different dimensions (coordinate directions) comes from thedivergence term ∇ · E′ . This gives the opportunity to construct one family of solutions in
2D/3D from the 1D solutions. Assuming that the solution E′ is identical in each dimension
and these identical components are denoted by Ẽ ′ , it follows immediately from the form of thedivergence in N-dimensional Cartesian space, that (11) can be written for each component as′∂ 2 Ẽ ′∂ Ẽ ′′ ∂ Ẽ
=+ g̃(τ ).
+
N
Ẽ
∂τ
∂x2
∂x
(The gauge field also needs identical components g̃, of course.) Under these special
assumptions equation (11) has indeed the same structure as the vector Burgers equation.
With a transformation similar to (10)E ′ = N Ẽ ′andg = N g̃2This simplest example of a gauge transformation represents indeed one of the symmetry groups for the solutions
to (16).
7

J. Phys. A: Math. Theor. 45 (2012) 455204

J Schönke

1

0

n and

-1

-2

-3

-4

-3

-2

-1

0

1

2

3

4

5

x
Figure 2. The concentration (27) and the potential (28) for the pure advection example. n(x, τ )(solid lines) and φ(x, τ ) (dotted lines) are shown for τ = 0 (red) and 2 (blue).

we then produce equation (13) again. Therefore, every solution in 1D (see section 5.1) can beused to construct new solutions in 2D/3D!
Assuming the N-dimensional space vector is denoted by x with Cartesian components xk ,(k = 1, . . . , N) and the 1D solution is given (e.g. like in section 5.1.2 with (22) as ψ (x, τ ), (23)
as E(x, τ ) and (24) as n(x, τ )). Then the electric field E has the components Ek = E(xk , τ )/N,and the concentration is given asNn(x, τ ) = n(x1 , . . . , xN , τ ) =1 !
n(xk , τ ),
N k=1

(29)

and analogous the potential readsN

φ(x, τ ) = φ(x1 , . . . , xN , τ ) =1!
φ(xk , τ ) .
N k=1

(30)

The concentration and the potential for the closed system example (cf section 5.1.2, with thesame numerical values A, B, C) are shown in 2D in figures 3 and 4.
8

J. Phys. A: Math. Theor. 45 (2012) 455204

J Schönke
n
0.2

n
0.2
0.0
-0.2
-0.4
-0.6
-0.8
-1.0

0.0

-0.2

5

-0.4
4
3

-0.6

2
1
x2

-0.8

0
-1
-2
-3
-4
-5

-5

-4

-3

-2

-1

0

1

2

3

4

5

-1.0

x1

Figure 3. The concentration (29) for the closed system example in 2D. n(x, τ ) = n(x1 , x2 , τ ) isshown for τ = 1.

6.2. Other cases equivalent to Burgers equation
As in the case we have just examined, certain assumptions about the solutions to (11) will leadto a formal equivalence of (11) and the vector Burgers equation. Another such example occurs
if the electric field E has only one direction but its absolute value E may vary arbitrarily inspace, i.e. (again in Cartesian coordinates) E(x) = (E(x), 0, 0). But this case has no physicalrelevance for our problem since this field is only conservative (i.e. can be written as ∇φ) if it
depends exclusively on x1 (the coordinate of the field direction). But this is nothing else but
the 1D case again.
7. Spherical symmetric case
Equation (11) written in spherical coordinates, assuming spherical symmetry (i.e. the angularcomponents are zero), is given as
)
*
∂Er′
1 ∂
2E ′
∂E ′
E′ ∂
g(τ )
= 2
r2 r − 2 r + 2r (r2 Er′ ) + 2 ,
(31)
∂τ
r ∂r
∂r
r
r ∂r
r
where Er′ is the radial component of the electric field and r is the radius. The gauge fieldg(τ )/r2 is restricted to the given radial dependency (since it has to be divergence free).

9J. Phys. A: Math. Theor. 45 (2012) 455204

J Schönke

φ

0.0

φ
0.0
-0.5
-1.0
-1.5
-2.0
-2.5

-0.5

-1.0

5
-1.5

4
3
2
1
x2

-2.0
0
-1
-2

-2.5

-3
-4
-5

-5

-4

-3

-2

-1

0

1

2

3

4

5

x1

Figure 4. The potential (30) for the closed system example in 2D. φ(x, τ ) = φ(x1 , x2 , τ ) is shownfor τ = 1.

A dimensional analysis shows that there exists a family of self-similar solutions to (31),i.e. the equation is invariant under the following scaling transformation:Er′ (r, τ ) = τ −1/2 H(s)ands = τ −1/2 r.

(32)

This is only possible if we set g(τ ) = τ −1/2 g0 with an arbitrary constant g0 . Introducing(32) into (31) we arrive at an ODE which can be further reduced by the transformation
H(s) = [K(s) + g0 ]/s2 resulting
*
)
s dK
d2 K
K + g02− K = 0.
(33)
+
− +
ds2
s2
s
2 ds

Several analytic methods to solve (33) fail and a Lie-group analysis (to find possible hidden
symmetries) yielded no useful results. But the reduction to an ODE offers at least simpler andindependent ways to produce numerical solutions to the problem.
7.1. The limit s ≫ 1
Additionally,
we consider the limit s ≫ 1 for (33). Recalling (32), this condition translates to√r ≫ τ , meaning ‘for large (and with time increasingly larger) radii’. In this limit the bracket

10

J. Phys. A: Math. Theor. 45 (2012) 455204

J Schönkein (33) reduces to s/2 (provided K is bounded by s2 for s → ∞, which is valid if Er′ does notdiverge for r → ∞) and we obtain the linear ODE

s dK
d2 K
+
− K = 0.
(34)
2
ds
2 ds
The general solution of (34) is given by the linear combination
,
*
) 2
*
+
)
√ s2
s
2
+ 1 [1 − erf(s/2)] − s e−s /4 + B
+1 ,
(35)
π
K=A
2
2
with free constants A and B and the error function erf(x). The first term A(...) in (35) decreases
to zero exponentially for s → ∞. Therefore, we need to consider only the second term B(...)
in our limit. Back substitution of K to get H and Er′ together with (10) finally gives the electric
field
√ ,
+
B
1
(B + g0 ) τEr (r, τ ) = −
√ +
z 2 τ
r2
and with (8) the concentration
1 B
n(r, τ ) = − 2 √ .
z r τThis solution describes a rarefaction process, the concentration (n > 0 if B < 0) decreases
steadily with time and the ion flux is pointing outward.
8. Conclusions

We have uncovered a connection between the single species Poisson–Nernst–Planck (PNP)
system and a nonlinear equation for the respective electric field. This equation (11) may beconsidered a member of the family of nonlinear, vector-valued drift-diffusion equations. The
further investigation of this equation regarding the differences compared to the vector Burgers
equation (the most famous member of this family) may provide additional insight into thestructure of (11) and its solutions.
A general procedure has been derived to find solutions to the problem in one dimension,
because in this case the equation can be brought to the form of the classical (scalar) Burgers
equation. Thus, we were able to give several time dependent solutions and realized different
types of boundary conditions, including zero flux (insulating) conditions.
The especially simple form of the divergence operator in Cartesian coordinates offered
the possibility to construct multidimensional solutions out of 1D solutions, where the spatial
dependency is identical in each coordinate direction. It is noted that this procedure is notpossible for the vector Burgers equation.
Finally, an exemplary approach to find solutions in the spherical symmetric case has beenoutlined, showing that there are various promising ways to search for further solutions to the
problem.
Since the PNP system has been investigated by the scientific community for more thana century, we want to emphasize the novelties in this paper through a comparison with therecent literature.
The analytic solutions given here are the first discovered time dependent (non-stationary)exact solutions at all. It is explicitly stated in the introduction of [6] that there are no known
exact analytical solutions to the full time dependent PNP system yet. Further, there is anextensive analysis and an in-depth historical review of the electrodiffusion problem in [10],where no indications of any known exact solutions are given. This is the case for [10] itself
and all thematically relevant references therein. Other recent publications about theoretical
11

J. Phys. A: Math. Theor. 45 (2012) 455204

J Schönke

aspects of electrodiffusion also do not contain any information on exact solutions, see, e.g.
[11–13].
A second novelty is the connection between the single species PNP system and the equation(11) for the electric field. If this connection would already be known, it might find its way intoa review like [10] but it is not given there or in any of the above mentioned references. Besides
journal articles, such a connection (if known) is likely to be found in respective textbooks onelectrodiffusion. One of the up-to-date books on this topic is [14] which includes a detailed
physical and mathematical analysis of the PNP system, but a connection to equation (11) could
not be found. Other classical textbooks like [15, 16] or [17] contain no information regarding
this point either.
Another novelty is the equation (11) itself. A detailed search in respective journals (e.g.
‘Journal of Differential Equations’ or ‘Journal of Mathematical Physics’) shows that theliterature on vector-valued partial differential equations (PDE) is generally rare except for
famous cases like e.g., the Navier–Stokes equation. A discussion of an equation like (11) was
not found. A further extremely rich source of information on nonlinear PDEs is the handbook[8] which contains thousands of equations and references to relevant publications from thelast centuries (fortunately including lots of the vast amount of Russian contributions to the
topic). Nothing comparable to equation (11) was found in [8]. Another brilliant textbook is[18], covering up-to-date developments in the area of nonlinear PDE’s with a huge number ofexamples from all fields of theoretical physics. Again, no references to an equation like (11)could be found.
Therefore, since (11) has not been discussed before (although it apparently describes aphysical process, namely the dynamics of the electric field for a single species electrolyte),
the other properties of (11) derived in this work are new, too. This includes the construction ofsolutions in the multidimensional case from 1D solutions (cf section 6) as well as the findings
for the spherical symmetric case in section 7.
Acknowledgments
This work was funded through a grant from the German Ministry of Education and Research(BMBF, 01GQ1003A).
References[1] Nernst W 1889 Die elektromotorische Wirksamkeit der Ionen Z. Phys. Chem. 4 129[2] Planck M 1890 Ueber die Erregung von Electricität und Wärme in electrolyten Ann. Phys. Chem. 39 161–86[3] Biler P, Hebisch W and Nadzieja T 1994 The Debye system: existence and large time behavior of solutions
Nonlinear Anal. 23 1189–209[4] Park J-H and Jerome J W 1997 Qualitative properties of steady-state Poisson–Nernst–Planck systems:
mathematical study SIAM J. Appl. Math. 57 609–30[5] Golovnev A and Trimper S 2010 Steady state solution of the Poisson–Nernst–Planck equations Phys. Lett.
A 374 2886–9[6] Golovnev A and Trimper S 2009 Exact solution of the Poisson–Nernst–Planck equations in the linear regime
J. Chem. Phys. 131 114903[7] Benzhuo L, Holst M J, McCammon J A and Zhou Y C 2010 Poisson–Nernst–Planck equations for simulating
biomolecular diffusion-reaction processes: I. Finite element solutions J. Comput. Phys. 229 6979–94[8] Polyanin A D and Zaitsev V F 2004 Handbook of Nonlinear Partial Differential Equations (Boca Raton, FL:
CRC Press)[9] Polyanin A D 2002 Handbook of Linear Partial Differential Equations for Engineers and Scientists (Boca
Raton, FL: CRC Press)
12

J. Phys. A: Math. Theor. 45 (2012) 455204

J Schönke[10] Bazant M Z, Thornton K and Ajdari A 2004 Diffuse-charge dynamics in electrochemical systems Phys. Rev.
E 70 021506[11] Kilic M, Bazant M and Ajdari A 2007 Steric effects in the dynamics of electrolytes at large applied voltages: I.
Double-layer charging Phys. Rev. E 75 021502[12] Rubinstein I, Zaltzman B, Futerman A, Gitis V and Nikonenko V 2009 Reexamination of electrodiffusion timescales Phys. Rev. E 79 021506[13] Jäckle J 2007 The causal theory of the resting potential of cells J. Theor. Biol. 249 445–63[14] Kontturi K, Murtomäki L and Manzanares J A 2008 Ionic Transport Processes: In Electrochemistry and
Membrane Science (Oxford: Oxford University Press)[15] Sten-Knudsen O 2002 Biological Membranes: Theory of Transport, Potentials and Electric Impulses(Cambridge: Cambridge University Press)[16] Weiss T F 1996 Cellular Biophysics: Transport vol 1 (Cambridge, MA: MIT Press)[17] Rubinstein I 1990 Electro-Diffusion of Ions (Philadelphia, MA: SIAM)[18] Debnath L 2011 Nonlinear Partial Differential Equations for Scientists and Engineers 3rd edn (Boston, MA:
Birkhäuser)

13