A Neumann Boundary Value Problem in Two-Ion Electro-Diffusion with Unequal Valencies

. In prior work, a series of two-point boundary value problems have been investigated for a steady state two-ion electro-diﬀusion model system in which the sum of the valencies ν + and ν − is zero. In that case, reduction is obtained to the canonical Painlev´e II equation for the scaled electric ﬁeld. Here, a physically important Neumann boundary value problem in the generic case when ν + + ν − (cid:54) = 0 is investigated. The problem is novel in that the model equation for the electric ﬁeld involves yet to be determined boundary values of the solution. A reduction of the Neumann boundary value problem in terms of elliptic functions is obtained for privileged valency ratios. A topological index argument is used to establish the existence of a solution in the general case, under the assumption ν + + ν − ≤ 0.


1.
Introduction. The theory of electro-diffusion originated in the liquid-junction theory of Nernst [1] and Planck [2]. It provides a macroscopic description of the transmission of charged particles through material barriers and has applications notably, in the modeling of biological membranes [3,4] and in electrochemistry [5].
Here, it proves convenient to partition the ions into m classes characterized by the same electric charge q j = q 0 ν j where q 0 is the unit of charge and ν j is a nonzero integral signed valency. The m-ion electro-diffusion model in steady régimes then reduces to the Nernst-Planck equations [6].
coupled to Gauss' equation Here, and where N i are the ion densities, N 0 is an arbitrary unit of ionic density, E is the electric field, T the temperature, κ the Boltzmann constant, andλ = [ κT /(4πq 2 0 N 0 ) ] 1/2 is the Debye length where is the dielectric constant. The c i are arbitrary constants of integration. The Painlevé analysis of the system (1)−(2) has recently been undertaken in [7].
Attention is restricted to the two-ion case and, in the notation of [8], we set n 1 = n + , n 2 = n − , ν 1 = ν + , ν 2 = ν − , c 1 = c + , c 2 = c − whence (1)−(2) yield The two-ion system (5)−(7) in the special case when ν + + ν − = 0 was investigated by Grafov and Chernenko [9] and independently by Bass [10]. An analogous system was subsequently derived independently in the context of semi-conductor theory by Kudryashov [11]. In both cases, reduction to the Painlevé II equation was obtained. This integrable connection has been exploited in [12] and [13] to apply a Bäcklund transformation sequentially to generate solutions of the Bass system. Addition of (5) and (6) together with use of (7) yields, on integration, where c = c + + c − and k is an arbitrary constant of integration. Elimination of n + between (7) and (8) then gives The condition that there is no net current in the junction yields [8] ν where D ± = u ± kT so that Accordingly, and (9) becomes If the junction has boundaries at x = 0 and x = δ then, on introduction of the scalings yields where now denotes the derivative with respect to x * and It is observed that (14) and (16) incorporate via c and k the boundary values n(0) = n + (0) + n − (0), n(1) = n + (1) + n − (1), together with y(0) and y(1). It is anticipated that the interface concentrations n + (0), n − (0), n + (1), n − (1) be known (Bass [10]). However, the boundary terms y(0) and y(1) dependent on the yet to be determined solution y remain. Insertion of (17), (18) into (16) yields The involvement of the boundary terms y(0) and y(1) in (20) poses a formidable impediment to its analysis and will be addressed in the sequal. It is anticipated that analogous procedures to those presented will be applicable mutatis mutandis to the corresponding (one-point) boundary value problem on the semi-infinite domain.

Reduction to elliptic integral formulation for privileged valency ratios.
Here, we consider the Neumann boundary value problem p (0) = p (1) = 0 (22) on the region [0,1] It is noted that the boundary conditions (22) imply, by virtue of Gauss' equation (7) that The latter conditions correspond to charge neutrality at the boundaries and can be imposed at the outset. These imply the necessary requirements on the interface concentration data. In particular, if, as in [8,10,11], ν + + ν − = 0 , it is seen that Here, we proceed with the generic case On introduction of the ansatz into the model equation (21) for the scaled electric field p, it is seen that, and (34) † It is noted that the corresponding arbitrary constants in the 3-ion case are required to be zero in Bass [14].
The Specializations with ν + ν − < 0, ν + + ν − = 0 These are or In the two canonical cases with valencies the electric field equation becomes, in turn or (27) where w is given in terms of elliptic integrals via or w = ± Bw + 2kw 2 /a 2 + Ew 4 (42) respectively. The positive sign is taken in the sequel.

Case I
Here, we proceed with the case (37) with ν + : ν − = 1 : −2 as obtained via Painlevé analysis in [7]. Insertion of (27) into (5) and (6) with c + = c − = 0 on integration yields (44) where the Gauss equation (7) shows that whence, the concentrations are given by The Neumann boundary conditions (22) require that Moreover, so that If we denote so that It is noted that the requirement (25) shows that the constants n + and n − are related by In the above, w is given via where it is required that the constants A, D and k be specified.
Case II Here, we consider the case (38) with ν + : ν − = 2 : −1 so that, in relations (43), (44) for the concentrations, whence The Neumann boundary conditions (43) require that be applied to Illustration Here, we set w(0) = w(1) = 1 (62) together with A = 1 so that, from (47), D = 2. The corresponding class of solutions of the two-ion system is then given by In the sequel, we apply an exact shooting method to a Neumann boundary value problem for the nonlinear equation (20) when reduction to an elliptic integral formalism is not available.
The proof will follow from a series of lemmas. First, we observe that, on searching for positive solutions, if we set z = y/γ, where γ = y(0), then (66)-(67) is equivalent to with conditions where In the original problem, the parameter l together with the other parameters, λ, D, and δ, are given, and we seek a solution y of (66) which satisfies the boundary conditions (67). As noted before, the fact that the differential equation (66) contains the unknown boundary values y(0) and y(1) makes the problem unconventional and cannot be adequately solved by traditional methods. To circumvent this difficulty, we consider the new problem (68) and (69). Here the parameter l is no longer given beforehand. Instead, we prescribe the parameters α and γ and proceed to solve the initial value problem (68) with the first two conditions in (69). By adjusting α appropriately, we can coerce the third condition in (69) to hold also, and then we recover l from the solution by solving (70). If the calculated l coincides with the given value of the parameter, we have our desired solution.
Thus, it suffices to find a pair (α, γ) such that the corresponding solution of the initial value problem satisfies z (1) = 0 and However, it is observed that an impediment arises on the fact that z may or may not be 'properly' defined throughout the entire interval [0, 1]. One situation is that z(t) can blow up to infinity before t reaches the endpoint 1 of the interval. Another possibility is that z(t) can hit the t-axis somewhere inside [0, 1] and then becomes negative subsequently. Let us define an 'endpoint' σ ∈ (0, 1] of the solution z in the following way: Thus, we are able to define a two-dimensional shooting operator T given by where L is given by Although physical considerations mean that we are interested only in those (α, γ) in the (open) first quadrant, we note parenthetically that T is defined also for γ = 0. We seek a pair (α, γ) such that T (α, γ) = (0, l): in that case it shall be seen that σ = 1 and γ > 0, and hence the corresponding z is a positive solution of (68)-(69) with α as in (70).
We shall make use of two comparison lemmas. The first is well known (see, for example, [16]) and the second is specific to our equation. and W (x) = F (x, W (x)) for x ∈ [0, a], where F is continuous and non-decreasing in the second variable for each fixed x ∈ [0, a]. If, in addition, it is assumed that Lemma 3.3. Let z be a solution of (68) with α > 0, and either γ > 0, or Then This implies that the righthand side term of (68) has positive z-derivative for z ≥ z(x 0 ). As it is also an increasing function of x, the result obviously follows. When γ = 0 and λ = αD, (68) becomes z = λz − αD. A direct computation gives The conclusion of the Lemma can then be verified directly.
Lemma 2 allows us to establish two fundamental facts about the shooting operator T . These are set down below: Proof. We start by proving the following claim: if z (σ) = 0, then 0 < z < 2 on [0, 1].
Continuity of T now follows from the standard continuous dependence result for ordinary differential equations.
Finally, if T (α, γ) = (0, l) then from the above claim, σ = 1 and the equality L = l implies that α satisfies (70), and y is therefore a solution of the original problem.
In order to prove the existence of a pair (α, γ) such that T (α, γ) = (0, l), we shall find a bounded domain C ⊂ (0, +∞) × [0, +∞) such that the topological index I of the curve T • ∂C, which is the image of the boundary of C under T , satisfies From the standard topological index theory, this implies that the equation T = (0, l) has at least one solution in C. More specifically, C shall be defined as the rectangle P QRS given by the vertices where α * and γ * are suitable constants to be chosen later.
Proof. If z (x 1 ) = 0 for some x 1 > x 0 , then either z(x 1 ) < z(x 0 ) or else from Lemma 3.3 we obtain z (x 1 ) > 0 and x 1 is a local minimum. In both cases, z attains a local maximum at some x 2 ∈ (x 0 , x 1 ) with z(x 2 ) ≥ z(x 0 ), and again we obtain z (x 2 ) > 0, a contradiction. Thus z does not vanish after x 0 , and then z (x) > 0 for x > x 0 . It is observed that, in particular, Case 1 cannot hold. Note that If Case 2 holds, or if Case 3 holds and z(1) ≥ 1, then (from the definition of L) L > 0, regardless of whether αD ≥ λ. or not. For the remaining case, when σ = 1 and z(1) < 1, we have The first inequality follows from Lemma 3.3, the equality in the second line from the differential equation (68) and the definition (71), and the inequality in the third line from the assumptions δ ≤ 1, z(1) < 1, and αD ≥ λ. Hence, we have L > 0.
The following lemmas provide a picture of the image of ∂C.
Lemma 3.6. The segment P Q is mapped one-to-one onto the segment P Q , where P = (0, 0) and Q = (−r, 0) for some r > 0.
Proof. Along the segment P Q, γ = 0 and z(x) is given by (73). Hence, T (P ) = (0, 0). The injectiveness property of T on P Q is actually not needed in the proof of Theorem 1. To prove it, we have to consider two cases. As we increase α from λ/D, initially we have Case 3, in which σ = 1. We can then use (73) to get to see that z (σ) is decreasing in α. However, after α reaches a critical value α 0 , i.e. for α > α 0 , Case 1 prevails. In this case, we multiply the differential equation z = λz − αD by z and integrate from x = 0 to x = σ to obtain Hence, z 2 (σ) is an increasing function of α. However, z (σ) is negative. Thus z (σ) is a decreasing function of α.
The next lemma shows that, except for the point P , the image of the segment P S lies in the first quadrant: Lemma 3.7. If α = λ/D and γ > 0 then z (σ) > 0 and L > 0.
Proof. If z is initially increasing, then 0 is a local minimum and Lemma 3.5 applies.
If z is initially decreasing, then as C ≤ 0 and z ≤ 0 in [0, ε] for some ε > 0, we obtain: Finally, if neither of the previous situations occurs, z would have more than one local minimum, and this contradicts Lemma 3.5.
The remaining two lemmas concern a convenient choice of α * and γ * . Proof. As z (σ) = 0, we are in Case 3 and it is readily seen that γ > 0. From Lemma 3.5, z cannot have a local minimum in [0, 1). Hence the global minimum of z is attained at the endpoint σ = 1, and we deduce that z (1) ≥ 0 and z is nonincreasing. From (68), we also deduce that z ≡ 1, then z(1) < 1 and which implies L > l. Again, the equality in the second line follows from (68) and (71).
In view of the preceding result, we fix a constant α * such that and proceed with the last lemma. Proof. We already know that T (S) lies on the first quadrant, so we may assume that α > λ/D. Then z < 0 in a neighborhood of 0 and z is initially decreasing. If z attains a local minimum at some point x 0 < σ, then Lemma 3.5 applies. Thus, it suffices to prove that z cannot be strictly decreasing all the time. Suppose, on the contrary, that z decreases strictly on [0, σ], then the term −Cγz(x)z (x) is nonnegative. Also, as z ≤ 1, we have that (1 − z 2 )z ≤ 2(1 − z), and hence from (68) we obtain: Next, fix a constant m such that α * D/(α * + λ) < m < 1. This can be done since D < 1 + 1 l , and the value of α * can be modified if necessary, as far as it satisfies (74).
Finally, define W as the solution of the linear problem Direct computation shows that where As α ≤ α * , it follows from the choice of m that C 1 > 0, and for γ * large enough we also have that ]. Since C 1 e γ * → +∞ as γ * → +∞, for γ * large enough, we get W (1) > 1, and this contradicts the fact that z ≤ 1.
Proof of Theorem 3.1. From the previous lemmas we conclude that the index of the curve T • ∂C relative to the point (0, l) is −1, and hence T (α, γ) = (0, l) for some (α, γ) ∈ C.  To help us to visualize the proof, we plot the image P Q R S of the rectangle P QRS under the mapping T in Figure 1, for a special case in which the physical parameters have been chosen to be λ = 1, D = 1, and δ = 1, and the constants α * and γ * in the definition of the rectangle P QRS have been chosen to be 6 and 4 respectively. The numerical experiment is done with the help of MATLAB.
The two distinctive kinks, one on the curve Q R and the other on R S represent the locations where there is a switch of case of the nature of the endpoint σ (between the three cases listed before Lemma 1). Let O denote the origin and A where the curve Q R intersects the L-axis. Then the topological index of P Q R S is 1 for every point that lies on the line segment OA. Hence, the original Neumann problem has a solution for these values of L. By increasing α * and γ * , more values of L will be covered (Lemmas 7 and 8).
If we reduce γ * to 3, the image of the side RS becomes R S . As shown in the figure, it does not lie entire in the first quadrant. This attests the fact that Lemma 8 only holds if γ * is sufficiently large.