Boundary layer solutions of Charge Conserving Poisson-Boltzmann equations: one-dimensional caseThe work is done when C.-C. Lee and T.-C. Lin were visiting the Department of Mathematics of Penn State University (PSU) during the summer of 2014. They express sincere thanks to PSU for all the hospitality and the great working environment. The authors also thank Professors Bob Eisenberg, Ping Sheng, Chiun-Chuan Chen and Rolf Ryham for numerous interesting discussions and valuable comments. The authors also thank the anonymous referees for their helpful comments and suggestions.The research of C.-C. Lee was partially supported by the Ministry of Science and Technology of Taiwan under the grants NSC-101-2115-M-134-007-MY2 and MOST-103-2115-M-134-001. Y. Hyon is partially supported by the National Institute for Mathematical Sciences grant funded by the Korean government (No. B21501). T.-C. Lin is partially supported by National Center of Theoretical Sciences (NCTS) and the Ministry of Science and Technology of Taiwan grant NSC-102-2115-M-002-015-MY4 and MOST-103-2115-M-002-005-MY3. C. Liu is partially supported by the NSF grants DMS-1109107, DMS-1216938, DMS-1159937 and DMS-1412005.

Boundary layer solutions of Charge Conserving Poisson-Boltzmann equations: one-dimensional case††thanks: The work is done when C.-C. Lee and T.-C. Lin were visiting the Department of Mathematics of Penn State University (PSU) during the summer of 2014. They express sincere thanks to PSU for all the hospitality and the great working environment. The authors also thank Professors Bob Eisenberg, Ping Sheng, Chiun-Chuan Chen and Rolf Ryham for numerous interesting discussions and valuable comments. The authors also thank the anonymous referees for their helpful comments and suggestions. The research of C.-C. Lee was partially supported by the Ministry of Science and Technology of Taiwan under the grants NSC-101-2115-M-134-007-MY2 and MOST-103-2115-M-134-001. Y. Hyon is partially supported by the National Institute for Mathematical Sciences grant funded by the Korean government (No. B21501). T.-C. Lin is partially supported by National Center of Theoretical Sciences (NCTS) and the Ministry of Science and Technology of Taiwan grant NSC-102-2115-M-002-015-MY4 and MOST-103-2115-M-002-005-MY3. C. Liu is partially supported by the NSF grants DMS-1109107, DMS-1216938, DMS-1159937 and DMS-1412005.

Chiun-Chang Lee Department of Applied Mathematics, National Hsinchu University of Education, No. 521, Nanda Road, Hsinchu 300, Taiwan, email(chlee@mail.nhcue.edu.tw)    Hijin Lee Department of Mathematical Sciences, Korea Advanced Institute of Science and Technology, Daejeon, Republic of Korea 305-701, email(hijin@kaist.ac.kr)    YunKyong Hyon Division of Mathematical Models, National Institute for Mathematical Sciences, Daejeon, Republic of Korea 305-811, email(hyon@nims.re.kr)    Tai-Chia Lin Department of Mathematics, Center for Advanced Study in Theoretical Sciences (CASTS), National Center for Theoretical Sciences (NCTS), National Taiwan University, Taipei, Taiwan 10617, email(tclin@math.ntu.edu.tw)    Chun Liu Department of Mathematics, Pennsylvania State University, University Park, PA 16802, USA, email(liu@math.psu.edu)
Abstract

For multispecies ions, we study boundary layer solutions of charge conserving Poisson-Boltzmann (CCPB) equations [50] (with a small parameter ) over a finite one-dimensional (1D) spatial domain, subjected to Robin type boundary conditions with variable coefficients. Hereafter, 1D boundary layer solutions mean that as approaches zero, the profiles of solutions form boundary layers near boundary points and become flat in the interior domain. These solutions are related to electric double layers with many applications in biology and physics. We rigorously prove the asymptotic behaviors of 1D boundary layer solutions at interior and boundary points. The asymptotic limits of the solution values (electric potentials) at interior and boundary points with a potential gap (related to zeta potential) are uniquely determined by explicit nonlinear formulas (cannot be found in classical Poisson-Boltzmann equations) which are solvable by numerical computations.

Key words. charge conserving Poisson-Boltzmann equations, boundary layer, multispecies ions

AMS subject classifications. ???

1 Introduction

Almost all biological activities involve transport in ionic solutions, which involves various couplings and interactions of multiple species of ions. Many complicated types of electrolytes involved in biological processes, such as those in ion channel proteins, certain amino acids (movable side chain) are crucial to the functions of these ion channels. The electrostatic properties involving multispecies (at least three species) ions can be fundamentally different to those with only one or two species [4, 33]. To see such difference, we study charge conserving Poisson-Boltzmann (CCPB) equation for multispecies ions which is derived from steady state Poisson-Nernst-Planck systems with charge conservation law, and is the surface potential model for the generation of a surface charge density layer related to electric double layers [30, 50]. For simplicity of analysis, we consider a physical domain with the simplest geometry, and represent CCPB equation as follows:

 −ϵ2ϕ′′=N∑i=1zie0mi∫1−1e−zie0kBTϕe−zie0kBTϕ forx∈(−1,1), \hb@xt@.01(1.1)

where is the total concentration of species with the valence , is the (electrical) potential, is the elementary charge, is the Boltzmann constant, and is the absolute temperature. The parameter  , where is the dielectric constant of the electrolyte, is the thermal voltage, is the length of the domain , and is the appropriate concentration scale (cf. [42]). Furthermore, is known as the Debye length and is of order for the physiological cases of interest (cf. [7]). Thus we may assume as a small parameter tending to zero. Similar equations to (LABEL:p4) can also be obtained by the other variational method [53].

Under suitable scales on and , we let ’s be the valences of anions, i.e., , and ’s be the valences of cations, i.e., , . Then the total concentrations of anions and cations are approximately given as () and (), respectively. Hence equation (LABEL:p4) can be transformed into

 ϵ2ϕ′′ϵ(x) = N1∑k=1akαk∫1−1eakϕϵ(y)dyeakϕϵ(x)−N2∑l=1blβl∫1−1e−blϕϵ(y)dye−blϕϵ(x) forx∈(−1,1),

where ’s and ’s satisfy and .

Most of the physical and biological systems possess the charge neutrality (zero net charge). One may assume the pointwise charge neutrality i.e. at all points the anion and cation charges exactly cancel in order to make calculations easier in a free diffusion system (cf. [19] p. 319). Here we replace the pointwise charge neutrality by a weaker hypothesis called the global electroneutrality being represented as

 Global Electro-neutrality: N1∑k=1akαk=N2∑l=1blβl, \hb@xt@.01(1.3)

which means that the total charges of anions and cations are equal, where ’s and ’s are the valences, and ’s and ’s are the concentrations of anions and cations, respectively. Consequently, the CCPB equation (LABEL:2-eqn1) may satisfy (LABEL:g-neutral).

When one deals with more general (realistic) situations, such as when there are more than two species involved in the solution, situations become more subtle and complicated. Note that the equation (LABEL:2-eqn1) has nonlocal dependence on which is essentially different from the classical Poisson-Boltzmann (PB) equation as follows:

 ϵ2ϕ′′ϵ(x)=N1∑k=1akαk2eakϕϵ(x)−N2∑l=1blβl2e−blϕϵ(x)forx∈(−1,1). \hb@xt@.01(1.4)

Here ’s and ’s are bulk concentration of anions and cations, respectively. In equation (LABEL:2-eqn1), ’s and ’s are for total concentration of anions and cations, respectively. For notation convenience, we use the same notations ’s and ’s in equations (LABEL:2-eqn1) and (LABEL:2-eqnw), but with different physical meaning. In this paper, we shall show different asymptotic behaviors of the CCPB equation (LABEL:2-eqn1) and the PB equation (LABEL:2-eqnw) for various constants satisfying (LABEL:g-neutral). The main goal of this paper is to compare the CCPB equation (LABEL:2-eqn1) and the PB equation (LABEL:2-eqnw) under the hypothesis (LABEL:g-neutral). Such a difference can be clarified in Theorems LABEL:2-thm:w and LABEL:2-thm5, see also, Remark LABEL:pbn-pb.

Boundary effects are important in a wide range of applications and provide formidable challenges [23, 25]. For CCPB equations, the main issue is how boundary conditions effect the solution values (electric potentials) at interior and boundary points. One may use the Neumann boundary condition for a given surface charge distribution and the Dirichlet boundary condition for a given surface potential (cf. [1]). Here we consider a Robin boundary condition [6, 24, 30, 29, 46, 47, 48, 51] for the electrostatic potential at is given by

 ϕϵ(1)+ηϵϕ′ϵ(1)=ϕ+0,ϕϵ(−1)−ηϵϕ′ϵ(−1)=ϕ−0. \hb@xt@.01(1.5)

where , are extrachannel electrostatic potentials and is the coefficient depending on the dielectric constant [36, 37], and related to the surface capacitance. The parameter ratio can be viewed as a measure of the Stern layer thickness, where and are the effective permittivity and the capacitance of the Stern layer, respectively (cf. [6]). Thus we may regard as the ratio of the Stern-layer width to the Debye screening length. Similar discussion can also be found in [13] and [41]. To see the influence of on the asymptotic behavior of ’s, we consider the limit to be either a non-negative constant or infinity.

Suppose i.e. , where is a non-negative constant. Then we show that the solution of (LABEL:2-eqn1) with (LABEL:2-eqn2) satisfies and for , where and can be uniquely determined by (LABEL:1.19.1)-(LABEL:1.19.3) which imply that the value is changed with respect to (see Figure LABEL:fig:layer). Moreover, the potential difference is decreasing to (cf. Theorem LABEL:2-thm5). Note that as the parameter goes to zero, the solution  has a boundary layer producing the potential gap affected by Stern and Debye (diffuse) layers and related to zeta potential (cf. [22]) which plays an important role in ionic fluids. However, for the PB equation (LABEL:2-eqnw), the value must be zero which is independent of and (cf. Theorem LABEL:2-thm:w). This shows the difference of the CCPB equation (LABEL:2-eqn1) and the PB equation (LABEL:2-eqnw) which can also be observed by numerical experiments (See Figure LABEL:fig:F_1 and Table LABEL:table:D_c in Section LABEL:2-mpnpsec:5). Furthermore, numerical computations give several conditions to let the profile of function to become monotone decreasing and increasing (Figure LABEL:fig:ct_3_12_allm and LABEL:fig:ct_4_all in Section LABEL:2-mpnpsec:5) and non-monotone (Figure LABEL:fig:ct_5_14_all_z in Section LABEL:2-mpnpsec:5).

In [30], we studied the CCPB equation (LABEL:2-eqn1) for case of , and , i.e., the case of one anion and one cation species with monovalence. In this case, equation (LABEL:2-eqn1) can be rewritten as

 ϵ2ϕ′′ϵ(x)=nϵ(x)−pϵ(x) forx∈(−1,1), \hb@xt@.01(1.6) nϵ(x)≡αeϕϵ(x)∫1−1eϕϵ(y)dyand pϵ(x)≡βe−ϕϵ(x)∫1−1e−ϕϵ(y)dy, \hb@xt@.01(1.7)

where and represent (pointwise) concentrations of anion and cation species, respectively. When holds (the electroneutral case), we had shown previously that for . Moreover, the CCPB equation (LABEL:2-eqn2014)-(LABEL:2013-1013-3) and the conventional PB equation have same asymptotic behavior (cf. Theorem 1.4 of [30]). In order for the readers to compare those with the results in the current paper, most results of [30] are summarized in Appendix. To certain degrees, it also justifies why in many situations, PB equation provides more or less expected solutions. On the other hand, we consider the non-electroneutral case, i.e. . Without loss of generality, we assume i.e. which means that the total concentration of anion species is less than that of cation species. Then we prove that for , but (cf. (LABEL:2014-0910-2)). This shows that electroneutrality holds true in the interior of , but non-electroneutrality occurs at the boundary points . Furthermore, the extra charges are accumulated near the boundary points (see Theorem  LABEL:NE-thm).

The mixture of monovalent and divalent ions such as , , and plays the most important roles for vital biological processes. For instance, opening and closing of ionic channels is accomplished by escape or entry of into the channels (cf. [18]). The voltage may depend on the concentration of (cf. [19]). Differences in ionic concentrations create a potential gap across the cell membrane that drives ionic currents (cf. [26] P. 34). To see how the voltage i.e. (electrical) potential depends on , we may use the equation (LABEL:2-eqn1) with , , and to describe the mixture of (or ), and ions, where , and . In Theorem LABEL:2-thm5 (ii), we prove that when the electro-neutrality holds, that is, , the solution of (LABEL:2-eqn1) satisfies for and is uniquely determined by (LABEL:1.19.1) and

 1−e3ccoshtecsinhc=β1β2=[Na+][Ca2+]>0 \hb@xt@.01(1.8)

where , and is a negative constant (see Remark LABEL:1.4). The formula (LABEL:ca1) shows that the interior potential (voltage) is increased if the boundary potential is fixed and the ratio is increased e.g.  is decreased and is fixed. Furthermore, Theorem LABEL:2-thm5 is also applicable to the other cases with multi-species ions including multivalent and polyvalent ions so the formula (LABEL:ca1) can be generalized to

 z−e(1+z)c(sinh(zt)sinht)2ecsinhc=β1β2, \hb@xt@.01(1.9)

for and (see Remark LABEL:1.4). Note that (LABEL:ca2) shows how the value depends on the value . Such a result cannot be found in the PB equation (LABEL:2-eqnw).

1.1 Asymptotic behavior of the PB equation (LABEL:2-eqnw)-(LABEL:2-eqn2)

The PB equation (LABEL:2-eqnw) with the boundary condition (LABEL:2-eqn2) can be regarded as the Euler-Lagrange equation of the energy functional

 EPBϵ[u]=12∫1−1(ϵ2|u′|2+f(u))dx+ϵ22ηϵ[(ϕ−0−u(−1))2+(ϕ+0−u(1))2], \hb@xt@.01(1.10)

for , where

 f(s)=N1∑k=1αkeaks+N2∑l=1βle−bls fors∈R. \hb@xt@.01(1.11)

For the PB equation (LABEL:2-eqnw) with the boundary condition (LABEL:2-eqn2), we study the asymptotic behavior of the solution of (LABEL:2-eqnw) as approaches zero. The boundary condition (LABEL:2-eqn2) plays a crucial role on the monotonicity of . Here we consider three cases for the signs of and : (a) , (b)  and (c) . Then the corresponding results are stated as follows:

Theorem 1.1

Assume . Let be the solution of equation (LABEL:2-eqnw) with the boundary condition (LABEL:2-eqn2). Then

1. For , exponentially converges to zero as goes to zero;

2. If , then is convex on and for . Moreover, there exists such that for , attains the minimum at an interior point of .

3. If , then is concave on and for . Moreover, there exists such that for , attains the maximum at an interior point of .

4. If , then is monotone on and .

5. If and , then uniquely determined by

 |ϕ+0−ˆt|=γ(f(ˆt)−f(0))1/2andmin{0,ϕ+0}≤ˆt≤max{0,ϕ+0}, \hb@xt@.01(1.12)

where is defined by (LABEL:id:2-015). Moreover, is decreasing in if and increasing in if .

1.2 The main results

In this section we present the main results, which are about the asymptotic behavior of the solution of (LABEL:2-eqn1) and (LABEL:2-eqn2) as goes to zero, in our research of CCPB equation. The CCPB equation (LABEL:2-eqn1) with the boundary condition (LABEL:2-eqn2) can be regarded as the Euler-Lagrange equation of the energy functional

 Eϵ[u] = ∫1−1ϵ22|u′|2dx+N1∑k=1αklog∫1−1eakudx+N2∑l=1βllog∫1−1e−bludx +ϵ22ηϵ[(ϕ−0−u(−1))2+(ϕ+0−u(1))2],

for . The existence and uniqueness for the solution of (LABEL:2-eqn1) and (LABEL:2-eqn2) is the following proposition:

Proposition 1.2

There exists a unique solution of the equation (LABEL:2-eqn1) with the boundary condition (LABEL:2-eqn2).

The proof of the above Proposition LABEL:2-thm1 can be easily obtained from the arguments of [30] (see Appendix therein) and [31].

Suppose and . Then Proposition LABEL:2-thm1 implies the solution of (LABEL:2-eqn1) and (LABEL:2-eqn2) must be trivial and . To study the nontrivial solution of (LABEL:2-eqn1) and (LABEL:2-eqn2), it is sufficient to assume . Replacing by for any constant , one may remark that the equation (LABEL:2-eqn1) is invariant. Consequently, without loss of generality, we may assume hereafter.

When , i.e., the global electroneutral case, Theorem LABEL:2-thm2 shows that is uniformly bounded to and that exponentially approaches zero in as tends to zero. Thus, it is expected that there exists a constant such that all interior values of tends to as goes to zero. Along with Lebesgue’s dominated convergence theorem, we have

 limϵ↓0∫1−1eakϕϵdx=2eakc,limϵ↓0∫1−1e−blϕϵdx=2−blc, \hb@xt@.01(1.14)

and then the energy functional (LABEL:2014-0821-energy) with approaches to the energy functional as follows (up to a constant independent of ):

 ˆEϵ[ϕϵ]=12∫1−1(ϵ2|ϕ′ϵ|2+f(ϕϵ−c))dx+ϵ22ηϵ[(ϕ−0−ϕϵ(−1))2+(ϕ+0−ϕϵ(1))2], \hb@xt@.01(1.15)

where is defined by (LABEL:id:2-015). Here we have used (by (LABEL:2014-0902)) and the approximation with to get

 log∫1−1eakϕϵdx∼12∫1−1eak(ϕϵ−c)dx+log(2eakc)−1as0<ϵ≪1.

Similarly, we have

 log∫1−1e−blϕϵdx∼12∫1−1e−bl(ϕϵ−c)dx+log(2e−blc)−1as0<ϵ≪1.

Therefore, we show that in the case of global electroneutrality (LABEL:g-neutral), the energy functional (LABEL:2014-0821-energy) approaches (LABEL:2014-0821-PB), which has the same form as the PB energy functional (LABEL:20140818-1).

The asymptotic behavior of ’s at boundary may depend on the scale of . Here we study two cases for the scale of : (i) and (ii) , where is a nonnegative constant. Then the relation between the boundary value limits and the interior value limit are demonstrated as follows:

Theorem 1.3

Assume and . Let be the solution of equation (LABEL:2-eqn1) with the boundary condition (LABEL:2-eqn2).Then

 limϵ↓0ϕϵ(−1)=−t,limϵ↓0ϕϵ(1)=t and limϵ↓0ϕϵ(x)=c for x∈(−1,1),

where and are determined as follows:

1. If   , then .

2. If   and , then uniquely solves the following equations:

 ϕ+0−t = γ(f(t−c)−f(0))1/2, \hb@xt@.01(1.16) f(t−c) = f(−t−c), \hb@xt@.01(1.17) |c|< t ≤ϕ+0. \hb@xt@.01(1.18)

Moreover, writing and in (ii), we have

• , and , where is uniquely determined by .

• and both are decreasing on .

Formally, using in as tends to zero, equation (LABEL:2-eqn1) may approach to the following PB equation:

 ϵ2ϕ′′ϵ(x)=N1∑k=1akαk2eak(ϕϵ(x)−c)−N2∑l=1blβl2e−bl(ϕϵ(x)−c)% forx∈(−1,1), \hb@xt@.01(1.19)

which may give results of Theorem LABEL:2-thm5 by formal asymptotic analysis. However, in this paper, we focus on rigorous mathematical analysis and provide the proof of Theorem LABEL:2-thm5 in Section LABEL:2-mpnpsec:2.

Theorem LABEL:2-thm5(i) shows that there is no boundary layer and uniformly in as if . Theorem LABEL:2-thm5(ii) assures the existence of boundary layers. Furthermore, Theorem LABEL:2-thm5 (ii-A) and (ii-B) represent the ratio of Stern screening length to the Debye screening length affects the boundary and interior potentials: (a) The decrease of results in the increase of (the potential difference between the boundary and interior); (b) If , the potential difference may approach zero. Notice that the formula (LABEL:1.25) is quite different from (LABEL:1.19.1)-(LABEL:1.19.3). This may show the difference between solutions of the CCPB equation (LABEL:2-eqn1) and the PB equation (LABEL:2-eqnw).

Remark 1.1
• Theorem LABEL:2-thm:w (ii) and (iii) show that as , the solution of the PB equation (LABEL:2-eqnw) may lose the monotonicity. However, the solution of the CCPB equation (LABEL:2-eqn1) always keeps the monotonicity (see Remark LABEL:2-rk1.2 (i)). This provides the difference between solutions of the CCPB equation (LABEL:2-eqn1) and the PB equation (LABEL:2-eqnw).

• For equation (LABEL:2-eqn1), the values (interior potential) and (boundary potential) depend on each other and satisfy precise formulas (LABEL:1.19.1)-(LABEL:1.19.3). However, for equation (LABEL:2-eqnw), interior potential and boundary potentials (determined by (LABEL:1.25)) are independent to each other.

Remark 1.2

When , , , and , we may get (LABEL:ca1) from (LABEL:1.19.1) and (LABEL:1.19.2). Moreover, (LABEL:ca2) can also be derived from (LABEL:1.19.1) and (LABEL:1.19.2) for the case that , , , and . By (LABEL:ca1) and (LABEL:ca2), it is easy to check that for . Then can be regarded as an decreasing function to . Consequently, by Theorem LABEL:2-thm5 (iv), is increasing to .

When , , and , further asymptotic behavior of near the boundary describing the boundary layers is stated as follows:

Theorem 1.4

Assume , , and . Under the same hypotheses of Theorem LABEL:2-thm2 and Theorem LABEL:2-thm5(ii), the asymptotic behavior of near the boundary can be represented by

 ϕ+1,ϵ(x)≤ϕϵ(x)≤ϕ+2,ϵ(x) forx∈(y+ϵ,1), \hb@xt@.01(1.20) ϕ−1,ϵ(x)≤ϕϵ(x)≤ϕ−2,ϵ(x) forx∈(−1,y−ϵ), \hb@xt@.01(1.21)

where satisfy , and

 ϕ+i,ϵ(x) = c+log{A+i,ϵ+B+i,ϵcsch2[C+i,ϵϵ(1−x)+logD+i,ϵ]}, \hb@xt@.01(1.22) ϕ−i,ϵ(x) = c+log{A−i,ϵ−B−i,ϵsech2[C−i,ϵϵ(1+x)+logD−i,ϵ]},i=1,2. \hb@xt@.01(1.23)

Here , , and , , are constants depending on such that , , and as goes to zero.

In the case of , , and , we may solve equation (LABEL:2-eqn1-asp) precisely and get the form of (LABEL:2-id:w13) and (LABEL:2-id:w13add) near and , respectively. One may remark how the values and affect the asymptotic behavior of near the boundary .

When (the non-electroneutral case), the asymptotic behavior for the solution , and of the equation (LABEL:2-eqn2014)-(LABEL:2013-1013-3) with the boundary condition (LABEL:2-eqn2) is stated as follows:

Theorem 1.5

Assume and . Let be the solution of the equation (LABEL:2-eqn2014)-(LABEL:2013-1013-3) with the boundary condition (LABEL:2-eqn2) and . Then

• When and , there exists a positive constant depending on and such that and

 α2−λϵ(κ)≤nϵ(x)≤pϵ(x)≤β2+λϵ(κ),forx∈[−1+ϵκ,1−ϵκ]. \hb@xt@.01(1.24)

Moreover, we have

 limϵ↓0nϵ(±1)=0and limϵ↓0ϵ2pϵ(±1)=(α−β)28, \hb@xt@.01(1.25) limϵ↓0supx∈[−1+ϵκ,1−ϵκ]∣∣∣nϵ(x)−α2∣∣∣ = limϵ↓0supx∈[−1+ϵκ,1−ϵκ]∣∣∣pϵ(x)−α2∣∣∣=0, \hb@xt@.01(1.26) limϵ↓0∫−1+ϵκ−1nϵ(x)dx = limϵ↓0∫11−ϵκnϵ(x)dx=0, \hb@xt@.01(1.27) limϵ↓0∫−1+ϵκ−1pϵ(x)dx = limϵ↓0∫11−ϵκpϵ(x)dx=β−α2. \hb@xt@.01(1.28)
• Let be any compact subset of . When is sufficiently small, the asymptotic expansion of in with the exact leading-order term and second-order term is given as follows:

 ϕϵ(x)−ϕϵ(±1)=log1ϵ2+log(α−β)24α+oϵ(1),forx∈K, \hb@xt@.01(1.29)

where denotes as a small quantity tending to zero as goes to zero.

Similar results also hold for .

Remark 1.3
• To exclude the boundary layer of with thickness (cf. Theorem 1.6 of [30]), we consider integrals of and over the interval , where is independent of . Note that and can be represented by (see (LABEL:2013-1013-3)), and that Theorem LABEL:NE-thm(ii) implies , and . This shows that as approaches zero, both the total concentrations of anion and cation species in the bulk tend to the same positive constant , while the total concentrations of anion and cation species in the region ( which is next to the boundary with thickness ) tend to zero and positive constant , respectively.

• We want to emphasize that Theorem LABEL:NE-thm(ii) improves the asymptotic behavior of shown in Theorem 1.5 of our previous paper [30].

Following results play important roles throughout this paper.

• Multiplying the equation (LABEL:2-eqn1) by , (LABEL:2-eqn1) may be transformed into

 ϵ22ϕ′2ϵ(x) = N1∑k=1αk∫1−1eakϕϵ(y)dyeakϕϵ(x)+N2∑l=1βl∫1−1e−blϕϵ(y)dye−blϕϵ(x) +Cϵ,

where is a constant depending on .

• Differentiating (LABEL:2-eqn1) to and multiplying it by ,

 ϵ2ϕ′′′ϵ(x)ϕ′ϵ(x) \hb@xt@.01(1.31) =⎛⎝N1∑k=1a2kαk∫1−1eakϕϵ(y)dyeakϕ