國立臺灣大學 數學系 預印本 Department of Mathematics, National Taiwan University

### www.math.ntu.edu.tw/ ~ mathlib/preprint/2013- 07.pdf

## A new approach to the Lennard-Jones potential and a new model: PNP-steric equations

### Tai-Chia Lin and Bob Eisenberg

### April 19, 2013

### model: PNP-steric equations

### Tai-Chia Lin

^{∗}### and Bob Eisenberg

^{†}**Abstract**

A class of approximate Lennard-Jones (LJ) potentials with a small parameter is found whose Fourier transforms have a simple asymptotic behavior as the parameter goes to zero.

When the LJ potential is replaced by the approximate LJ potential, the total energy functional becomes simple and exactly the same as replacing the LJ potential by a delta function. Such a simple energy functional can be used to derive the Poisson-Nernst-Planck equations with steric eﬀects (PNP-steric equations), a new mathematical model for the LJ interaction in ionic solutions. Using formal asymptotic analysis, stability and instability conditions for the 1D PNP-steric equations with the Dirichlet boundary conditions for one anionic and cationic species are expressed by the valences, diﬀusion constants, ionic radii and coupling constants.

This is the ﬁrst step to study the dynamics of solutions of the PNP-steric equations.

**1** **Introduction**

The Lennard-Jones (LJ) potential, a well-known mathematical model for the interaction be- tween a pair of ions, has important applications in many ﬁelds of biology, chemistry and physics (cf. [46]). Such a capable model can be represented by

*Ψ (x) =* *C*1

*r*^{12} *−C*2

*r*^{6} for *r =|x| > 0, x ∈ R*^{d}*,*

*where C*_{1}*, C*_{2} *are positive constants related to ﬁnite ion size, and d≤ 3 is the spatial dimension*
(cf. [36]). The inverse twelfth-power term is the repulsive term of Ψ and the inverse sixth-power
*term is the attractive term of Ψ. The LJ potential Ψ can be extended to an l− m LJ potential*
given as follows:

Ψ*l,m**(x) =* *C*1

*r*^{m}*−C*2

*r** ^{l}* for

*r =|x| > 0, x ∈ R*

^{d}*,*

*where l and m are any positive constants with m > l > d. All the mathematical arguments of the*
*LJ potential Ψ here can easily be generalized to the l− m LJ potential for m > l > d, C*1*> 0, and*
*C*_{2}*≥ 0.*

To compute the energy of ions interacting by the LJ potential, the following energy functional is considered:

*E**LJ**[c**i**, c**j*] =

∫∫

R^{d}*×R*^{d}

*Ψ (x− y) c**i**(x) c**j**(y) dxdy ,* (1.1)
*for nonnegative functions c**i* *and c**j* which denote the distribution (concentration) functions of the
*ith and jth ion species (cf. [29]). Similar energy functionals for Coulomb interactions can be found*
*in [9]. Note that the energy functional E**LJ* *is for two ion species if i* *̸= j, but for only one ion*
*species if i = j. Because the LJ potential Ψ is singular at the origin, the functional E**LJ* is a
singular integral deﬁned by

*E**LJ**[c**i**, c**j*] = lim

*σ**→0*

∫∫

R^{d}*×R*^{d}

*Ψχ**σ**(x− y) c**i**(x) c**j**(y) dxdy ,* (1.2)

*∗*Department of Mathematics, National Taiwan University, Taida Institute for Mathematical Sciences (TIMS),
No.1, Sec.4, Roosevelt Road, Taipei 106, Taiwan, email: tclin@math.ntu.edu.tw

*†*Department of Molecular Biophysics Physiology Rush Medical Center, 1653 West Congress, Parkway, Chicago,
IL 60612, USA, email: beisenbe@rush.edu

1

*for c**i**, c**j* *∈ L*^{2}+

(R* ^{d}*)

={

*f* *∈ L*^{2}(
R* ^{d}*)

*: f≥ 0*}

*, where χ**σ**= χ**σ**(z) is the characteristic function of*
the exterior ball{

*z∈ R** ^{d}*:

*|z| > σ*}

*. In [16, 17, 18, 28], the energy functional E*_{LJ}*with C*_{2}= 0 is
*used, but here both C*2*= 0 and C*2*> 0 are considered.*

Conventionally, the system of Poisson-Nernst-Planck (PNP) equations, a model of ion transport, plays a crucial role in the study of many physical and biological phenomena (cf. [2, 3, 6, 7, 10, 12, 14, 15, 30, 32, 40, 41, 42, 45, 48, 52, 53, 59, 60]). However, when ions are crowded in a narrow channel, the PNP equations become unreliable because the ion-size eﬀect becomes important, but the PNP equations represent ions as point particles without size (cf. [1, 4, 8, 11, 18, 21, 22, 26, 33, 37, 38, 50, 51, 56, 58, 61]). Hence the PNP equations need to be modiﬁed in order to describe solutions where the ion-size eﬀect is important. Biological solutions are mixtures containing divalents in which ion size eﬀects are always important (cf. [17, 19, 20]).

To modify the PNP equations, many eﬀorts have been made to combine the energy functional of the PNP equations with the other exclusion terms which may come from liquid state theory and density functional theory [5, 23, 34, 38, 40, 47]. Being related to liquid state theory, the Lennard-Jones (LJ) potential is often used as an approximate model of the van der Waals force (cf. [29, 46]). However, it seems that no one had ever used the LJ potential to modify the PNP equations before the pioneering works of Eisenberg and Liu who derived the PNP equations with size eﬀects i.e. equations (2.6)-(2.8) by combining the repulsive term of the LJ potential and the energy functional of the PNP equations (cf. [16, 28]).

Equations (2.6)-(2.8) importantly generalize the PNP equations and numerically simulate the selectivity of ion channels which can not be obtained by solving the PNP equations alone (cf. [62, 63, 64, 65, 66]). Nevertheless, due to imposing the LJ potential, equations (2.6)-(2.8) become a complicated system of diﬀerential-integral equations having no numerical eﬃciency (cf. [17, 28]) and allowing no theoretical result either. The goal of this paper is to approximate the LJ potential and simplify equations (2.6)-(2.8) into the PNP-steric equations, a new mathematical model for the LJ interaction in ionic solutions. Instead of singular integrals of equations (2.6)-(2.8), the PNP- steric equations are composed of the PNP equations and nonlinear diﬀerential terms with coupling constants. Numerical simulations of the PNP-steric equations are presented in [25], which shows the numerical eﬃciency to simulate the selectivity of ion channels previously studied in Monte Carlo simulations with results comparable to a wide range of experiments.

The main diﬃculty of numerical simulation of the PNP-steric equations is how to choose cou-
pling constants suitably. This motivates us to study a simple case (one anion and one cation
species) of the PNP-steric equations with the Dirichlet boundary conditions using suitable asymp-
*totic expansions to see the eﬀect of coupling constants g**ij*’s. The formal asymptotic analysis gives
*the stability and instability conditions represented by g**ij*’s which can be regarded as the ﬁrst step
in a series of analyses. (see Section 4 and 5). More theoretical results will be done soon which may
be useful for the choice of coupling constants in order to do further numerical simulations on ion
channels.

This paper has two major parts: one is the approximation of the LJ potential using band-limited functions and the other is the stability and instability conditions for the PNP-steric equations with the Dirichlet boundary conditions using asymptotic expansions. The approximate LJ potentials are introduced in Section 2.1 and the detailed mathematical arguments are stated in Section 3.

The PNP-steric equations are derived in Section 2.2 and the stability and instability conditions for the PNP-steric equations are proved in Section 4 and 5, respectively.

**2** **Preliminaries**

**2.1** **The approach to Lennard-Jones potential**

Because the LJ potential Ψ does not have a Fourier transform, it is diﬃcult to study the
*energy functional E**LJ* *directly. When the spatial frequency variable ξ is bounded, the Fourier*
transform d*Ψχ**σ**(ξ) tends to inﬁnity as σ goes to zero, but by the Riemann-Lebesgue Lemma,*

*|ξ|→∞*lim *Ψχ*d*σ**(ξ) = 0 for all σ > 0. Hence the asymptotic behavior of the Fourier transform dΨχ**σ**(ξ)*
is dominated at bounded spatial frequencies and negligible at high spatial frequencies. This gives
a reason to cut oﬀ high spatial frequencies in order to see the asymptotic behavior of the Fourier
transform d*Ψχ*_{σ}*(ξ). From a physical and biological point of view, it seems obvious that particularly*

small spatial features cannot be used by evolution to produce biological function and so should
*be absent in analysis. The high spatial frequency cut-oﬀ function φ**σ* is a band-limited function
deﬁned by

*φ**σ**(x) = (1− χ**σ**−γ**(ξ))*^{v}*,* (2.1)

*for x, ξ∈ R*^{d}*, where v denotes the inverse Fourier transform. Obviously,*
ˆ

*φ**σ**(ξ) = 1− χ**σ*^{−γ}*(ξ) =*

{ 1 if *|ξ| ≤ σ*^{−γ}*,*

0 if *|ξ| > σ*^{−γ}*,* (2.2)

*for ξ* *∈ R*^{d}*, where hat denotes the Fourier transform, 0 < γ < 1 is a constant independent of σ*
*and ξ.*

Generically, band-limited functions can handle spatial information locally in frequency domains
and play important roles in Fourier analysis of mathematics and have many applications to en-
gineering, physics and statistics (cf. [31, 35, 55]). Physically, it is obvious that the band-limited
approximation is reasonable if not inevitable. After all, derivations of PNP replace Maxwell’s
equations with Poisson’s equation, and neglect ’capacitance to ground’ (capacitive coupling be-
tween ionic solutions and nearby ground planes) always present in experiments. These and many
other eﬀects at high frequencies are better attenuated to zero than approximated irrationally (i.e.,
without known error bounds) as they are in treatments with unlimited bandwidth. Here a new ap-
*proach to the LJ potential Ψ uses the spatially band-limited function φ** _{σ}*to deﬁne the approximate
potential Ψ as follows:

Ψ_{σ}*(z) = (Ψχ*_{σ}*) ⋆ φ*_{σ}*(z)* for *z∈ R*^{d}*,* (2.3)
*where the asterisk is the standard convolution, φ**σ*is the spatially band-limited function deﬁned in
*(2.1) and χ**σ* is the characteristic function of the exterior ball*{z ∈ R** ^{d}*:

*|z| > σ}. Note that χ*

*σ*

*(z)*and 1

*− χ*

*σ*

^{−γ}*(ξ) are the characteristic functions of*{

*z∈ R** ^{d}*:

*|z| > σ*} and {

*ξ∈ R** ^{d}*:

*|ξ| ≤ σ*

*} , respectively, and both of them extend to the entire spaceR*

^{−γ}

^{d}*as the small parameter σ goes to zero.*

*As σ goes to zero, χ**σ**∼ 1, 1 − χ**σ*^{−γ}*∼ 1, and*

Ψc_{σ}*(ξ) = dΨχ*_{σ}*(ξ)φ*c_{σ}*(ξ) = dΨχ*_{σ}*(ξ) [1− χ**σ*^{−γ}*(ξ)]∼ bΨ (ξ) .*
which implies Ψ*σ**∼ Ψ. Hence formally,*

∫

R^{d}

*c*_{i}*(x) (Ψ*_{σ}*⋆ c*_{j}*) (x) dx∼*

∫

R^{d}

*c*_{i}*(x) (Ψ ⋆ c*_{j}*) (x) dx = E*_{LJ}*[c*_{i}*, c** _{j}*]

*This shows how we approximate the energy functional E*

*.*

_{LJ}*The approximate energy functional E** _{LJ,σ}* is deﬁned by

*E*

*LJ,σ*

*[c*

*i*

*, c*

*j*] =

∫∫

R^{d}*×R*^{d}

Ψ*σ**(x− y) c**i**(x) c**j**(y) dxdy =*

∫

R^{d}

*c**i**(x) (Ψ**σ* *⋆ c**j**) (x) dx.*

*As σ goes to zero, the functional E*_{LJ,σ}*tends to the functional E** _{LJ}* if the following hypothesis
holds:

(H) lim

*σ**→0+*

[∫

R^{d}*Ψχ*d_{σ}*(ξ) [1− χ**σ**−γ**(ξ)]c*b_{i}*(ξ)c*b_{j}*(ξ)dξ−*∫

R^{d}*Ψχ*d_{σ}*(ξ)c*b_{i}*(ξ)c*b_{j}*(ξ)dξ*
]

*= 0.*

(see Proposition 3.1). Here the meaning of ‘approximate’ is diﬀerent from that of conventional
approximation theory. Note that the characteristic function 1*−χ**σ** ^{−γ}* truncates the high frequencies

*|ξ| > σ** ^{−γ}*, but still preserves the spatial frequencies of order

*|ξ| ∼ σ*

^{−θ}*(for all 0 < θ < γ) tending*

*to inﬁnity as σ goes to zero. By standard theorems of Fourier analysis (cf. [54]), the functional*

*E*

*LJ,σ*

*satisﬁes E*

*LJ,σ*

*[c*

*i*

*, c*

*j*] =∫

R^{d}*c*b*i**(ξ) \*Ψ*σ**⋆ c**j**(ξ) dξ =*∫

R^{d}*c*b*i**(ξ) c*Ψ*σ**(ξ)c*b*j**(ξ) dξ On the other hand,*
(2.2) and (2.3) imply cΨ*σ**(ξ)∼ C*1 *ω**d*

12*−d**σ*^{d}^{−12}*as σ goes to zero, where ω**d* *is the surface area of d*
dimensional unit ball (see the proof of Proposition 3.2). Hence

*E**LJ,σ**[c**i**, c**j*]*∼ C*1

*ω**d*

12*− dσ*^{d}^{−12}

∫

R^{d}*c*b*i**(ξ)c*b*j**(ξ) dξ = C*1

*ω**d*

12*− dσ*^{d}^{−12}

∫

R^{d}

*c**i**(x)c**j**(x) dx .*
*This shows that the functional E** _{LJ,σ}*is asymptotically close to another functional ˜

*E*

*as follows:*

_{LJ,σ}*E**LJ,σ**[c**i**, c**j*]*∼ ˜E**LJ,σ**[c**i**, c**j**] = C*1*S**σ*

∫

R^{d}

*c**i**(x)c**j**(x) dx ,* (2.4)

*where S**σ* = _{12}^{ω}_{−d}^{d}*σ*^{d}* ^{−12}*. Thus the functional ˜

*E*

*LJ,σ*can be regarded as an approximate energy

*functional to E*

*LJ*

*. Note that the constant C*1 comes from the repulsive term of the LJ potential,

*but the constant C*2 (for the attractive term of the LJ potential) does not aﬀect the leading term of the asymptotic behavior of cΨ

*σ*. This may support the work of [16, 17, 18, 28] which only use the repulsive term of the LJ potential to describe ionic interactions. However, the eﬀects of the attractive terms remain to be investigated. They may have been selected by evolution to produce qualitative behavior of importance to biology. Note that the functional ˜

*E*

*LJ,σ*is much simpler than

*the functional E*

*LJ*and it can be expressed as

*E*˜_{LJ,σ}*[c*_{i}*, c** _{j}*] =

∫∫

R^{d}*×R*^{d}

*C*_{1}*S*_{σ}*δ*_{0}*(x− y)c**i**(x) c*_{j}*(y) dxdy*

*being the same as replacing the LJ potential Ψ by a delta function C*_{1}*S*_{σ}*δ*_{0}in the energy functional
*E*_{LJ}*, where δ*_{0}(*·) is the standard delta function concentrating at the origin. Such a simple energy*
functional ˜*E** _{LJ,σ}* can be used to derive the PNP-steric equations as a new model of ionic solutions.

**2.2** **PNP-steric equations**

The Poisson-Nernst-Planck (PNP) equations consist of the Nernst-Planck equations coupled with the Poisson equation being expressed as follows:

*∂c*_{i}

*∂t* =_{k}^{D}^{i}

*B**T**∇ ·*(

*c*_{i}*∇*^{δE}_{δc}^{pnp}* _{i}* )

*,* *i = 1,· · · , N,*

*−∇ · (ε∇ϕ) = ρ*0+

∑*N*
*i=1*

*z**i**ec**i**,*

*where E**pnp*is the energy functional of the PNP equations given by
*E*_{pnp}*[c*_{1}*,· · · , c**N**, ϕ] :=*∫

R^{d}

[
*k*_{B}*T*

∑*N*
*i=1*

*c*_{i}*log c** _{i}*+

^{1}

_{2}(

*ρ*_{0}+

∑*N*
*i=1*

*z*_{i}*ec** _{i}*
)

*ϕ*
]

*dx*

*Here N is the number of ion species, c**i* *is the distribution function, D**i* is the diﬀusion constant,
*and z**i* *is the valence of the ith ion species, respectively. Besides, ϕ is the electrostatic potential,*
*ε is the dielectric constant, ρ*0 *is the permanent (ﬁxed) charge density of the system, k**B* is the
*Boltzmann constant, T is the absolute temperature and e is the elementary charge. More precisely,*
the PNP equations are denoted as

*∂c**i*

*∂t* *= D**i**∇ ·*(

*∇c**i*+_{k}^{z}^{i}^{e}

*B**T**c**i**∇ϕ*)

*,* *i = 1,· · · , N,*

*−∇ · (ε∇ϕ) = ρ*0+

∑*N*
*i=1*

*z**i**ec**i*

(2.5)

*To include the hard sphere repulsion of ions, the energy functional E**pnp*is modiﬁed by adding
the energy functional

∑*N*
*i,j=1*

*E*_{LJ}*[c*_{i}*, c*_{j}*] with the constants C*_{1}= ^{1}_{2}*ϵ*_{ij}*(a*_{i}*+ a** _{j}*)

^{12}

*and C*

_{2}= 0, where

*c*

*i*

*and c*

*j*

*are the distribution (concentration) functions of the ith and jth ion species with the radii*

*a*

*i*

*and a*

*j*

*, respectively. Then the modiﬁed energy functional E*

*mpnp*becomes

*E**mpnp**[c*1*,· · · , c**N**, ϕ] :=*

∫

R^{d}

(
*k**B**T*

∑*N*
*i=1*

*c**i**log c**i*+1
2

(
*ρ*0+

∑*N*
*i=1*

*z**i**ec**i*

)
*ϕ*

)
*dx*

+1 2

∑*N*
*i,j=1*

∫∫

R^{d}*×R*^{d}

*ϵ*_{ij}*(a*_{i}*+ a** _{j}*)

^{12}

*|x − y|*^{12} *c*_{i}*(x) c*_{j}*(y) dxdy ,*

Using energy variational analysis (cf. [16, 28]), the modiﬁed PNP equations are the following equations

*∂c*_{i}

*∂t* = _{k}^{D}^{i}

*B**T**∇ ·*(

*c**i**∇*^{δE}_{δc}^{mpnp}* _{i}* )

*,* *i = 1,· · · , N,*

i.e. *∂c**i*

*∂t* +*∇ · J**i* *= 0 ,* *i = 1,· · · , N,* (2.6)
coupled with the Poisson equation

*−∇ · (ε∇ϕ) = ρ*0+

∑*N*
*i=1*

*z**i**ec**i**,* (2.7)

*where ﬂux J**i* is

*J**i*=*−D**i**∇c**i**−D**i**c**i*

*k**B**Tz**i**e∇ϕ −* *D**i**c**i*

*k**B**T*

∑*N*
*j=1*

*∇*

∫

R^{d}

*ϵ**ij**(a**i**+ a**j*)^{12}

*|x − y|*^{12} *c**j**(y) dy .* (2.8)
However, equations (2.6)-(2.8) are diﬃcult to investigate theoretically because they are partial
diﬀerential-integral equations with singular integrals. Moreover, due to the eﬀect of high (Fourier)
frequencies, the numerical computations of equations (2.6)-(2.8) may lose accuracy and become
ineﬃcient (cf. [28]).

Instead of equations (2.6)-(2.8), a simple model can be derived by replacing the energy func- tional

∑*N*
*i,j=1*

*E**LJ**[c**i**, c**j**] (with C*1=^{1}_{2}*ϵ**ij**(a**i**+a**j*)^{12}*and C*2= 0) by the approximate energy functional

∑*N*
*i,j=1*

*E*˜*LJ,σ**[c**i**, c**j*] = 1

2*ϵ**ij**(a**i**+ a**j*)^{12}*S**σ*

∫

R^{d}

*c**i**(x)c**j**(x) dx*

*deﬁned in (2.4). Then the energy functional E**mpnp* can be approximated by

*E**σ**[c*1*,· · · , c**N**, ϕ] :=*

∫

R^{d}

(
*k**B**T*

∑*N*
*i=1*

*c**i**log c**i*+1
2

(
*ρ*0+

∑*N*
*i=1*

*z**i**ec**i*

)
*ϕ*

)
*dx*

+1 2

∑*N*
*i,j=1*

*ϵ**ij**(a**i**+ a**j*)^{12}*S**σ*

∫

R^{d}

*c**i**(x) c**j**(x)dx ,*

*which is much simpler than the energy functional E** _{mpnp}*. Using energy variational analysis (cf. [16,
28]), a new model called the PNP-steric equations is expressed by

*∂c*_{i}

*∂t* =_{k}^{D}^{i}

*B**T**∇ ·*(

*c**i**∇*^{δE}_{δc}^{σ}* _{i}*)

*,* *i = 1,· · · , N,*

i.e. *∂c**i*

*∂t* +*∇ · J**i**= 0 ,* (2.9)

coupled with the Poisson equation

*−∇ · (ε∇ϕ) = ρ*0+

∑*N*
*i=1*

*z**i**ec**i**,* (2.10)

*where ﬂux J**i* is

*J**i*=*−D**i**∇c**i**−D*_{i}*c*_{i}

*k**B**Tz**i**e∇ϕ −* *D*_{i}*c*_{i}*k**B**TS**σ*

∑*N*
*j=1*

*ϵ**ij**(a**i**+ a**j*)^{12}*∇c**j**.* (2.11)

*Here the symmetry ϵ**ij* *= ϵ**ji*has been assumed for notation convenience. The PNP-steric equations
are of convection-diﬀusion type with the following energy dissipation law:

*d*

*dtE**σ**[c*1*,· · · , c**N**, ϕ] =−*

∫

R^{d}

∑*N*
*i=1*

*D**i**c**i*

*k*_{B}*T|∇ (k**B**T log c**i**+ z**i**eϕ + µ**i*)*|*^{2}*dx*

*where µ**i*= _{δc}^{δ}

*i*

∑*N*
*j,k=1*

*E*˜*LJ,σ**[c**j**, c**k**] = S**σ*

∑*N*
*j=1*

*ϵ**ij**(a**i**+ a**j*)^{12}*c**j* is the chemical potential. On the other
hand, the PNP-steric equations have more nonlinear diﬀerential terms than the PNP equations
so they can simulate the selectivity of ion channels eﬃciently (cf. [25]). Note that the selectivity
of ion channels can not be found by using the PNP equations. This shows that the PNP-steric
equations are more capable than the PNP equations and much simpler than equations (2.6)-(2.8).

*For the case of two species ions (i.e. N = 2) with one anionic and cationic species, the index*
*j = 1, 2 is replaced by j = n, p for notation convenience and the PNP-steric equations (2.9)-(2.11)*
are presented as

*∂c**n*

*∂t* *=D**n*

[

*∇ ·*
(

*∇c**n*+ *z**n**e*
*k*_{B}*Tc**n**∇ϕ*

)

+ ˜*S**σ**∇ · (g**nn**c**n**∇c**n**+ g**np**c**n**∇c**p*)
]

*,* (2.12)

*∂c**p*

*∂t* *=D**p*

[

*∇ ·*
(

*∇c**p*+ *z**p**e*
*k*_{B}*Tc**p**∇ϕ*

)

+ ˜*S**σ**∇ · (g**pp**c**p**∇c**p**+ g**np**c**p**∇c**n*)
]

*,* (2.13)

*−∇ · (ε∇ϕ) =ρ*0*+ z**n**ec**n**+ z**p**ec**p**,* (2.14)

*for x∈ R*^{d}*, t > 0, where the function c*_{n}*is for anion concentration, c** _{p}* is for cation concentration,

*S*˜

*σ*=

_{k}^{1}

*B**T* *S**σ**, g**np* *= ϵ*12*(a*1*+ a*2)^{12}*, g**nn* *= ϵ*11*(2a*1)^{12} *and g**pp* *= ϵ*22*(2a*2)^{12}. Due to the spatial
*dimension d≤ 3, the constant ˜S**σ**∼ σ*^{d}^{−12}*becomes a large quantity tending to inﬁnity as σ goes*
to zero. Let ˜*ε =* ^{k}_{e}* ^{B}*2

^{T}*ε, ˜δ = 1/ ˜S*

*σ*

*, τ = t/˜δ, ˜ϕ(x, τ ) =*

_{k}

^{e}*B**T**ϕ(x, t) and ˜c**j**(x, τ ) = c**j**(x, t) for j = n, p.*

Then the equations (2.12)-(2.14) are transformed into 1

*D*_{n}

*∂˜c**n*

*∂τ* =˜*δ∇ ·*(

*∇˜c**n**+ z** _{n}*˜

*c*

_{n}*∇ ˜ϕ*)

+*∇ · (g**nn*˜*c*_{n}*∇˜c**n**+ g** _{np}*˜

*c*

_{n}*∇˜c*

*p*

*) ,*(2.15) 1

*D*_{p}

*∂˜c**p*

*∂τ* =˜*δ∇ ·*(

*∇˜c**p**+ z**p**c*˜*p**∇ ˜ϕ*)

+*∇ · (g**pp*˜*c**p**∇˜c**p**+ g**np*˜*c**p**∇˜c**n**) ,* (2.16)

*−∇ · (˜ε∇ ˜ϕ) =˜ρ*0*+ z**n*˜*c**n**+ z**p*˜*c**p**,* (2.17)
*for x∈ R*^{d}*, τ > 0, where ˜ρ*_{0}=^{ρ}_{e}^{0}. Note that ^{∂c}_{∂t}* ^{j}* = ˜

*δ*

^{−1 ∂˜c}

_{∂τ}

^{j}*for j = n, p, and τ*

*∼ 1 is equivalent to*

*t∼ ˜δ = 1/ ˜S*

*σ*

*being a small time scale. For notation convenience, removing tilde and replacing τ*

*by t, the equations (2.15)-(2.17) become*

1
*D**n*

*∂c*_{n}

*∂t* *=δ∇ · (∇c**n**+ z*_{n}*c*_{n}*∇ϕ) + ∇ · (g**nn**c*_{n}*∇c**n**+ g*_{np}*c*_{n}*∇c**p**) ,* (2.18)
1

*D**p*

*∂c*_{p}

*∂t* *=δ∇ · (∇c**p**+ z*_{p}*c*_{p}*∇ϕ) + ∇ · (g**pp**c*_{p}*∇c**p**+ g*_{np}*c*_{p}*∇c**n**) ,* (2.19)

*−∇ · (ε∇ϕ) =ρ*0*+ z*_{n}*c*_{n}*+ z*_{p}*c*_{p}*,* (2.20)

*for x∈ R*^{d}*, t > 0.*

Instead of the entire space R* ^{d}*, here the spatial domain is considered as an one-dimensional

*(d = 1) interval (−1, 1), and then the equations (2.18)-(2.20) are changed as the following equations:*

1
*D**n*

*∂c**n*

*∂t* *=δ* *∂*

*∂x·*
(*∂c**n*

*∂x* *+ z**n**c**n*

*∂ϕ*

*∂x*
)

+ *∂*

*∂x*
(

*g**nn**c**n*

*∂c**n*

*∂x* *+ g**np**c**n*

*∂c**p*

*∂x*
)

*,* *for x∈ (−1, 1) , t > 0 ,*
1

*D**p*

*∂c*_{p}

*∂t* *=δ* *∂*

*∂x*
(*∂c*_{p}

*∂x* *+ z*_{p}*c*_{p}*∂ϕ*

*∂x*
)

+ *∂*

*∂x*
(

*g*_{pp}*c*_{p}*∂c*_{p}

*∂x* *+ g*_{np}*c*_{p}*∂c*_{n}

*∂x*
)

*,* *for x∈ (−1, 1) , t > 0 ,*

*−ε∂*^{2}*ϕ*

*∂x*^{2} *=ρ*_{0}*+ z*_{n}*c*_{n}*+ z*_{p}*c*_{p}*,* *for x∈ (−1, 1) , t > 0 .*

*Setting ε = 1 and ρ*0= 0, these equations become
1

*D**n*

*∂c*_{n}

*∂t* *=δ* *∂*

*∂x*
(*∂c*_{n}

*∂x* *+ z*_{n}*c*_{n}*∂ϕ*

*∂x*
)

+ *∂*

*∂x*
(

*g*_{nn}*c*_{n}*∂c*_{n}

*∂x* *+ g*_{np}*c*_{n}*∂c*_{p}

*∂x*
)

*,* *for x∈ (−1, 1) , t > 0 ,*
(2.21)
1

*D*_{p}

*∂c**p*

*∂t* *=δ* *∂*

*∂x*
(*∂c**p*

*∂x* *+ z**p**c**p*

*∂ϕ*

*∂x*
)

+ *∂*

*∂x*
(

*g**pp**c**p*

*∂c**p*

*∂x* *+ g**np**c**p*

*∂c**n*

*∂x*
)

*,* *for x∈ (−1, 1) , t > 0 ,*
(2.22)

*−∂*^{2}*ϕ*

*∂x*^{2} *=z*_{n}*c*_{n}*+ z*_{p}*c*_{p}*,* *for x∈ (−1, 1) , t > 0 ,*

(2.23)
*where δ is a small parameter tending to zero. For simplicity, we consider the following Dirichlet*
boundary conditions:

*c**n*=*−*^{B}_{z}_{n}^{1}*, c**p*= ^{B}_{z}^{1}

*p**, ϕ = B*3 *as x = 1 ,*
*c**n*=*−*^{B}_{z}_{n}^{2}*, c**p*= ^{B}_{z}^{2}

*p**, ϕ = B*4 *as x =−1 ,*

(2.24)

*where B*_{k}*, k = 1,· · · , 4 are constants. Here we only focus on the equations (2.21)-(2.23) with the*
boundary condition (2.24) but not (2.18)-(2.20) on the entire space. Using the formal asymptotic
analysis, we may get the stability and instability conditions of the equations (2.21)-(2.23) with the
*boundary condition (2.24) which show the eﬀect of coupling constants g** _{ij}*’s.

To see the solution of the equations (2.21)-(2.23), the following asymptotic expansions are used:

*c**n* *= c**n,0**+ δc**n,1**+ δ*^{2}*c**n,2*+*· · · ,*
*c**p**= c**p,0**+ δc**p,1**+ δ*^{2}*c**p,2*+*· · · ,*
*ϕ = ϕ*_{0}*+ δϕ*_{1}*+ δ*^{2}*ϕ*_{2}+*· · · .*
*Then the zeroth order solution (c**n,0**, c**p,0**, ϕ*0) satisﬁes

1
*D**n*

*∂c**n,0*

*∂t* = *∂*

*∂x*
(

*g**nn**c**n,0*

*∂c**n,0*

*∂x* *+ g**np**c**n,0*

*∂c**p,0*

*∂x*
)

*,* (2.25)

1
*D**p*

*∂c**p,0*

*∂t* = *∂*

*∂x*
(

*g**pp**c**p,0*

*∂c**p,0*

*∂x* *+ g**np**c**p,0*

*∂c**n,0*

*∂x*
)

*,* (2.26)

*−∂*^{2}*ϕ*_{0}

*∂x*^{2} *= z*_{n}*c*_{n,0}*+ z*_{p}*c*_{p,0}*,* (2.27)

*for x∈ (−1, 1) and t > 0. As g**np*= 0, both (2.25) and (2.26) become the standard porous medium
*equations (cf. [57]) hence for g**np**̸= 0, the equations (2.25) and (2.26) can be regarded as a coupled*
system of porous medium equations.

The steady state equations of (2.25) and (2.26) are denoted as
(*g*_{nn}*c*_{n,0}*c*^{′}_{n,0}*+ g*_{np}*c*_{n,0}*c*^{′}* _{p,0}*)

_{′}*(x) = 0 ,* (2.28)

(*g**pp**c**p,0**c*^{′}_{p,0}*+ g**np**c**p,0**c*^{′}* _{n,0}*)

_{′}*(x) = 0 ,* (2.29)

*for x* *∈ (−1, 1), where a prime mark (′) is denoted as diﬀerentiation to the spatial variable x.*

Electroneutrality (which means the charge of anions is equal to that of cations) holds in most
biological systems (cf. [67]). A particular assumption called pointwise electroneutrality which
means the charge of anions is everywhere equal to that of cations appears in the zeroth order
equation [24] as a pleasingly natural physical approximation. To ﬁnd solutions of (2.28)-(2.29)
*with pointwise electroneutrality, it is natural to assume that c*_{n,0}*and c** _{p,0}*have the following form

*c*

*n,0*

*(x) = z*

*p*

*w(x) ,*

*c*

*p,0*

*(x) =−z*

*n*

*w(x)*for

*x∈ (−1, 1) ,*(2.30)

*so the total charge of c*

_{n,0}*and c*

_{p,0}*becomes zero i.e. z*

_{n}*c*

_{n,0}*(x) + z*

_{p}*c*

_{p,0}*(x) = 0 for x*

*∈ (−1, 1).*

*Note that z*_{n}*and z** _{p}* are the associated valences, the charge on an individual ion. To solve the

*equations (2.28) and (2.29), the function w = w(x) must satisfy* (
*w*^{2})_{′′}

*(x) = 0 for x∈ (−1, 1),*
*and hence w = w(x) =√*

*K*_{1}*x + K*_{0}*, where K*_{0} *and K*_{1} are constants determined by the Dirichlet
boundary conditions as follows:

*−z**n**c**n,0**(1) = z**p**c**p,0**(1) = B*1*,* *−z**n**c**n,0*(*−1) = z**p**c**p,0*(*−1) = B*2*,* (2.31)
*and B*1*and B*2are the positive constants for the Dirichlet boundary conditions in (2.24). Because
*w(x) =* *√*

*K*1*x + K*0 *for x* *∈ [−1, 1], then by (2.30) and (2.31), it is obvious that* *√*

*K*1*+ K*0 =
*w(1) =* _{−z}^{B}^{1}

*n**z** _{p}* and

*√*

*K*0*− K*1*= w(−1) =* _{−z}^{B}_{n}^{2}_{z}* _{p}* which imply

*K*_{1}= *B*_{1}^{2}*− B*2^{2}

*2z*_{n}^{2}*z*^{2}* _{p}* and

*K*

_{0}=

*B*

^{2}

_{1}

*+ B*

^{2}

_{2}

*2z*

_{n}^{2}

*z*

^{2}

_{p}*.*Hence

*w(x) =*

√*(B*^{2}_{1}*− B*2^{2}*)x + B*_{1}^{2}*+ B*_{2}^{2}

*√*2*|z**n**z**p**|* *> 0* for *x∈ [−1, 1] ,* (2.32)
*being non-constant if B*1*̸= B*2*. Combining (2.27) and (2.30), the electric potential ϕ*0satisﬁes the
following equation:

*−ϕ** ^{′′}*0

*(x) = 0*for

*x∈ (−1, 1) .*(2.33)

The equation (2.33) is a standard diﬀerential equation solved as follows:

*ϕ*0*(x) = a*0*x + b*0 for *x∈ (−1, 1) ,* (2.34)
*where a*0 *and b*0 *are constants determined by the boundary conditions of ϕ*0 given in (2.24). The
*ﬁrst order solution (c**n,1**, c**p,1**, ϕ*1*) = (c**n,1**(x, t), c**p,1**(x, t), ϕ*1*(x, t)) satisﬁes*

1
*D*_{n}

*∂c**n,1*

*∂t* = *∂*

*∂x*
(*∂c**n,0*

*∂x* *+ z*_{n}*c*_{n,0}*∂ϕ*0

*∂x*
)

+ *∂*

*∂x*
[

*g*_{nn}*∂(c**n,0**c**n,1*)

*∂x* *+ g** _{np}*
(

*∂c*

*p,0*

*∂x* *c*_{n,1}*+ c*_{n,0}*∂c**p,1*

*∂x*
)]

*,*
(2.35)
1

*D**p*

*∂c**p,1*

*∂t* = *∂*

*∂x*
(*∂c**p,0*

*∂x* *+ z**p**c**p,0*

*∂ϕ*0

*∂x*
)

+ *∂*

*∂x*
[

*g**pp*

*∂(c**p,0**c**p,1*)

*∂x* *+ g**np*

(*∂c**n,0*

*∂x* *c**p,1**+ c**p,0*

*∂c**n,1*

*∂x*
)]

*,*
(2.36)

*−∂*^{2}*ϕ*1

*∂x*^{2} *= z**n**c**n,1**+ z**p**c**p,1**,* (2.37)

*for x∈ (−1, 1), t > 0. By (2.30), the equations (2.35) and (2.36) become*
1

*D**n*

*∂c**n,1*

*∂t* *= z**p**(w*^{′}*+ z**n**wϕ*^{′}_{0})* ^{′}*+

*∂*

*∂x*
[

*g**nn**z**p*

*∂(w c**n,1*)

*∂x* *+ g**np*

(

*−z**n**w*^{′}*c**n,1**+ z**p**w∂c**p,1*

*∂x*
)]

*,* (2.38)
1

*D**p*

*∂c*_{p,1}

*∂t* =*−z**n**(w*^{′}*+ z**p**wϕ*^{′}_{0})* ^{′}*+

*∂*

*∂x*
[

*−z**n**g**pp*

*∂(w c** _{p,1}*)

*∂x* *+ g**np*

(

*z**p**w*^{′}*c**p,1**− z**n**w∂c*_{n,1}

*∂x*
)]

*,* (2.39)
*for x∈ (−1, 1), t > 0.*

*Let ψ = z**n**c**n,1**+ z**p**c**p,1**and φ =−z**n**c**n,1**+ z**p**c**p,1**i.e. c**n,1*= _{2z}^{1}

*n**(ψ− φ) and c**p,1*=_{2z}^{1}

*p**(ψ + φ).*

Then after some algebraic calculations, the equations (2.38) and (2.39) can be transformed into
*the equations for ψ and φ as follows:*

*∂ψ*

*∂t* *= (D*_{n}*− D**p**) z*_{n}*z*_{p}*w*^{′′}*(x) + z*_{n}*z*_{p}*(D*_{n}*z*_{n}*− D**p**z*_{p}*) [w(x)ϕ*^{′}_{0}*(x)]** ^{′}* (2.40)
+ ( ˜

*d*

*np*

*− ˜d)∂*

^{2}

*[w(x)φ]*

*∂x*^{2} + (¯*g + ¯g**np*) *∂*

*∂x[w*^{′}*(x)ψ] + (¯g− ¯g**np*) *∂*

*∂x*
[

*w(x)∂ψ*

*∂x*
]

*,*
and

*∂φ*

*∂t* =*− (D**n**+ D**p**) z**n**z**p**w*^{′′}*(x)− z**n**z**p**(D**n**z**n**+ D**p**z**p**) [w(x)ϕ*^{′}_{0}*(x)]** ^{′}* (2.41)
+ (¯

*g*

*np*+ ¯

*g)∂*

^{2}

*[w(x)φ]*

*∂x*^{2} + ( ˜*d**np**− ˜d)* *∂*

*∂x[w*^{′}*(x)ψ]− ( ˜d**np*+ ˜*d)* *∂*

*∂x*
[

*w(x)∂ψ*

*∂x*
]

*,*

where ¯*g**np**, ˜d**np**, ¯g and ˜d are constants deﬁned as follows:*

¯

*g**np* =^{1}_{2}*g**np**(D**p**z**p**− D**n**z**n**) ,*
*d*˜* _{np}* =

^{1}

_{2}

*g*

_{np}*(D*

_{p}*z*

_{p}*+ D*

_{n}*z*

_{n}*) ,*

¯

*g* =^{1}_{2}*(D*_{n}*z*_{p}*g*_{nn}*− D**p**z*_{n}*g*_{pp}*) ,*
*d*˜ =^{1}_{2}*(D*_{n}*z*_{p}*g*_{nn}*+ D*_{p}*z*_{n}*g*_{pp}*) .*

(2.42)

*The equations (2.40) and (2.41) depend on valences z**i**, diﬀusion constants D**i* and coeﬃcients
*g*_{ij}*∼ (a**i**+ a** _{j}*)

^{12}

*(i, j = n, p) related to ionic radii. In particular, these equations can represent the*diﬀerence between NaCl, KCl and CaCl

_{2}.

For notational convenience, the equations (2.40) and (2.41) can be expressed as

*∂*

*∂t*

*ψ*

*φ*

=

*f*1

*f*2

+(

*L*11 *L*12

*L*21 *L*22

) *ψ*

*φ*

* ,* (2.43)

*where f*1 *and f*2 are external force functions given by

*f*1*(x) = (D**n**− D**p**) z**n**z**p**w*^{′′}*(x) + z**n**z**p**(D**n**z**n**− D**p**z**p**) [w(x)ϕ*^{′}_{0}*(x)]*^{′}*,* (2.44)
*f*_{2}*(x) =− (D**n**+ D*_{p}*) z*_{n}*z*_{p}*w*^{′′}*(x)− z**n**z*_{p}*(D*_{n}*z*_{n}*+ D*_{p}*z*_{p}*) [w(x)ϕ*^{′}_{0}*(x)]*^{′}*,* (2.45)
*and L*_{ij}*, i, j = 1, 2 are diﬀerential operators deﬁned by*

*L*11*ψ* = (¯*g + ¯g**np*) *∂*

*∂x[w*^{′}*(x)ψ] + (¯g− ¯g**np*) *∂*

*∂x*
[

*w(x)∂ψ*

*∂x*
]

*,* (2.46)

*L*_{12}*φ* = ( ˜*d*_{np}*− ˜d)∂*^{2}*[w(x)φ]*

*∂x*^{2} *,* (2.47)

*L*21*ψ* = ( ˜*d**np**− ˜d)* *∂*

*∂x[w*^{′}*(x)ψ]− ( ˜d**np*+ ˜*d)* *∂*

*∂x*
[

*w(x)∂ψ*

*∂x*
]

*,* (2.48)

*L*_{22}*φ* = (¯*g** _{np}*+ ¯

*g)∂*

^{2}

*[w(x)φ]*

*∂x*^{2} *.* (2.49)

When ¯*g > ¯g** _{np}*, the asymptotic stability of system (2.43) (i.e. the equations (2.40) and (2.41))
is proved in Corollary 4.2 (see Section 4). Thus the equations (2.38) and (2.39) have asymptotic
stability if the condition ¯

*g > ¯g*

_{np}*holds. However, if g*

_{nn}*= g*

_{pp}*= g*

*= 0, then system (2.43)*

_{np}*becomes ψ*

_{t}*= f*

_{1}

*and φ*

_{t}*= f*

_{2}

*which imply (ψ, φ) = (ψ*

_{0}

*+ tf*

_{1}

*, φ*

_{0}

*+ tf*

_{2}

*) for t > 0, where ψ*

_{0}

*= ψ|*

_{t=0}*and φ*0

*= φ|*

_{t=0}*are initial conditions of ψ and φ, respectively. Hence the asymptotic stability of*

*system (2.43) is gone. This shows that the eﬀect of coupling constants g*

*ij*’s changes the asymptotic stability of system (2.43).

Note that the operator

( *L*11 *L*12

*L*21 *L*22

)

=

( *L*11 0
0 *L*22

)

becomes diagonal if ˜*d**np**− ˜d =*
*d*˜* _{np}*+ ˜

*d = 0 i.e. ˜d*

*= ˜*

_{np}*d = 0. By (2.42), the condition ˜d*

*= ˜*

_{np}*d = 0 is equivalent to D*

_{p}*z*

_{p}*+ D*

_{n}*z*

*= 0*

_{n}*and D*

_{n}*z*

_{p}*g*

_{nn}*+ D*

_{p}*z*

_{n}*g*

_{pp}*= 0 being fulﬁlled if D*

_{p}*= D*

_{n}*, z*

*=*

_{p}*−z*

*n*

*and g*

_{nn}*= g*

*. In Section 5, the*

_{pp}*case of D*

_{p}*= D*

_{n}*= D > 0, z*

*=*

_{p}*−z*

*n*

*= z > 0 and g*

_{nn}*= g*

_{pp}*= g > 0 (i.e. symmetry electrolytes)*

*is considered in order to get the instability condition g < g*

*.*

_{np}**3** **The approximate LJ potential**

*In this section, we study the approximation of the functional E** _{LJ}* deﬁned in (1.2) and written as
follows:

*E*_{LJ}*[c*_{i}*, c** _{j}*] = lim

*σ**→0*

∫∫

R^{d}*×R*^{d}

[(Ψ_{12}*− Ψ*6*) χ*_{σ}*] (x− y) c**i**(x) c*_{j}*(y) dxdy,* (3.1)

*for c**i**, c**j* *∈ L*^{2}+

(R* ^{d}*)

={

*f* *∈ L*^{2}(
R* ^{d}*)

*: f≥ 0*}

. Here Ψ12*(z) = C*1*|z|** ^{−12}* is the repulsion term and
Ψ6

*(z) = C*2

*|z|*

^{−6}*is the attraction term of the LJ potential, and χ*

*σ*

*= χ*

*σ*

*(z) is the characteristic*function of the exterior ball{

*z∈ R** ^{d}*:

*|z| > σ*}

satisfying

*χ**σ**(z) =*

1 if *|z| > σ,*
0 if *|z| ≤ σ.*

The main idea here is to approximate \Ψ*m**χ**σ**(ξ) (m = 6, 12) the Fourier transform of the kernel*
Ψ*m**χ**σ* by \Ψ*m**χ**σ**φ*ˆ*σ**(ξ) in the frequency space, where hat denotes Fourier transform and φ**σ* is a
band-limited function satisfying

ˆ

*φ**σ**(ξ) = 1− χ**σ*^{−γ}*(ξ)* for *ξ∈ R*^{d}*.* (3.2)
Here 1*− χ**σ** ^{−γ}* is deﬁned by

1*− χ**σ*^{−γ}*(ξ) =*

1 if *|ξ| ≤ σ*^{−γ}*,*
0 if *|ξ| > σ*^{−γ}*,*

*where σ > 0 is a small parameter tending to zero, and 0 < γ < 1 is a constant independent of ξ*
*and σ. Note that 1− χ**σ** ^{−γ}* denotes the characteristic function of the ball

*{ξ ∈ R*

*:*

^{d}*|ξ| ≤ σ*

^{−γ}*}*expanding the entire spaceR

^{d}*as σ goes to zero.*

By the standard formulas of Fourier analysis, it is obvious that

Ψ\_{m}*χ*_{σ}*(ξ) ˆφ*_{σ}*(ξ) =*[(Ψ* _{m}*\

*χ*

_{σ}*) ⋆ φ*

_{σ}*](ξ) ,*(3.3)

*for ξ*

*∈ R*

*, where the asterisk denotes convolution. The approximate kernel of Ψ*

^{d}*is given as follows:*

_{m}*K**m,σ**(z) = (Ψ**m**χ**σ**) ⋆ φ**σ**(z)* for *z∈ R*^{d}*,* *m = 6, 12 ,* (3.4)
*and the approximate energy functional E**LJ,σ* given by

*E**LJ,σ**[c**i**, c**j*] :=

∫∫

R^{d}*×R*^{d}

*(K**12,σ**− K**6,σ**) (x− y)c**i**(x)c**j**(y)dxdy,* (3.5)

*for c**i**, c**j* *∈ L*^{2}+

(R* ^{d}*)

. By standard theorems of Fourier analysis (cf. [54]), it is easy to check that
the functions Ψ*m**χ**σ* *and φ**σ**are of L*^{2}(R* ^{d}*)

*∩ L*

*(R*

^{∞}

^{d}*) which implies K*

*m,σ*

*∈ L*

^{2}(R

*)*

^{d}*∩ L*

*(R*

^{∞}*). The*

^{d}*following proposition is for the approximation of E*

*LJ,σ*

*to E*

*LJ*.

**Proposition 3.1. Assume c***i**, c**j* *∈ L*^{2}+(R^{d}*) satisfy the following hypothesis:*

*(H)* lim

*σ**→0+*

[∫

R^{d}*Ψχ*d_{σ}*(ξ) [1− χ**σ**−γ**(ξ)]c*b_{i}*(ξ)c*b_{j}*(ξ)dξ−*∫

R^{d}*Ψχ*d_{σ}*(ξ)c*b_{i}*(ξ)c*b_{j}*(ξ)dξ*
]

*= 0 ,*
*where 0 < γ < 1 is a constant independent of σ and ξ. Then*

*σ*lim*→0+**{E**LJ,σ**[c**i**, c**j*]*− E**LJ**[c**i**, c**j*]*} = 0 .* (3.6)
*Note that the hypothesis (H) is achievable at least for functions c**i*’s satisfying *c*b*i**(ξ) = 0 for*

*|ξ| > σ*^{−γ}*and i = 1,· · · , N, which means all the high frequencies of c**i*’s have been cut oﬀ. Due
to the strong singularity of Ψ, d*Ψχ**σ**(ξ) has no asymptotic behavior like the right side of (3.12) as*
*σ goes to zero, especially for|ξ| ∼ σ** ^{−1}* (see Remark 3.3). This motivates us to replace d

*Ψχ*

*σ*

*(ξ) by*

*Ψχ*d

_{σ}*(ξ) [1− χ*

*σ*

^{−γ}*(ξ)] which is a kind of truncation on the ξ variable and has a simple asymptotic*formula (3.12). The truncation may lose the eﬀect of high frequencies

*|ξ| > σ*

^{−γ}*(0 < γ < 1), but*still involve a large part of the steric eﬀects because it keeps the eﬀect of frequencies like

*|ξ| ∼ σ*

^{−α}*for 0 < α < γ tending to inﬁnity as σ goes to zero. Please note that the main goal of our works is*to simplify the model of Liu and Eisenberg [16, 28]. Here we present a simpliﬁed model which is easy to study and useful to understand the selectivity of ion channels [25].