# On the asymptotic limit of steady state Poisson–Nernst–Planck equations with
steric effects

Jhih-Hong Lyu , Tai-Chia Lin Department of Mathematics, National Taiwan
University, Taipei 10617, Taiwan (d06221001@ntu.edu.tw).Department of
Mathematics, National Taiwan University, Taipei 10617, Taiwan; National Center
for Theoretical Sciences, Mathematics Division, Taipei 10617, Taiwan
(tclin@math.ntu.edu.tw).

###### Abstract

When ions are crowded, the effect of steric repulsion between ions becomes
significant and the conventional Poisson–Boltzmann (PB) equation (without
steric effect) should be modified. For this purpose, we study the asymptotic
limit of steady state Poisson–Nernst–Planck equations with steric effects
(PNP-steric equations). By the assumptions of steric effects, we transform
steady state PNP-steric equations into a PB equation with steric effects (PB-
steric equation) which has a parameter $\Lambda$ and positive constants
$\lambda_{i}$’s (depend on the radii of ions and solvent molecules). The
nonlinear term of PB-steric equation is mainly determined by a Lambert type
function which represents the concentration of solvent molecules. As
$\Lambda=0$, the PB-steric equation becomes the conventional PB equation but
as $\Lambda>0$, a large $\Lambda$ makes the steric repulsion (between ions and
solvent molecules) stronger. This motivates us to find the asymptotic limit of
PB-steric equation as $\Lambda$ goes to infinity. Under the Robin (or Neumann)
boundary condition, we prove theoretically and numerically that the PB-steric
equation has a unique solution $\phi_{\Lambda}$ which converges to solution
$\phi^{*}$ of a modified PB (mPB) equation as $\Lambda$ tends to infinity. Our
results show that the limiting equation of PB-steric equation (as $\Lambda$
goes to infinity) is a mPB equation which has the same form (up to scalar
multiples) as those mPB equations in [9, 10, 28, 32, 33, 34, 37]. Therefore,
the PB-steric equation can be regarded as a generalized model of mPB
equations.
Key words. PNP-steric equations, mPB equation, steric effects
AMS subject classifications. 35C20, 35J60, 35Q92

##  1 Introduction

As a well-known model of ion transport, the Poisson–Nernst–Planck (PNP)
equations are important in the study of many physical and biological phenomena
(cf. [4, 5, 12, 17, 18]). The PNP equations represent ions as point particles
without size. However, in crowded ions, steric repulsions may appear due to
ion sizes so the PNP equations should be modified (cf. [1, 3, 6, 7, 8, 11, 14,
15, 19, 21, 23, 24, 26, 27, 28, 30, 38, 39]). To include the steric effect of
ion sizes, we use the approximate Lennard–Jones potential to develop a new
model called the Poisson–Nernst–Planck equations with steric effects (PNP-
steric equations) which have parameters depending on ion radii (cf. [35]).
Recently, the PNP-steric equations have been used to study the selectivity and
gating of ion channels [16, 24, 25], the anodic dissolution at high current
densities [29], and the slit-shaped nanopore conductance [14]. Hence the PNP-
steric equations become a useful model which describes ion transport with
steric effects.

For the mixture of ions and solvent molecules, the PNP-steric equations can be
denoted as

$\displaystyle\frac{\partial c_{i}}{\partial
t}+\nabla\cdot\mathcal{J}_{i}=0,\quad\quad i=0,1,\dots,N,$ (1.1)
$\displaystyle\mathcal{J}_{i}=-D_{i}\nabla
c_{i}-\frac{D_{i}c_{i}}{k_{B}T}z_{i}\text{e}\nabla\phi-\frac{D_{i}c_{i}}{k_{B}T}\sum_{j=0}^{N}g_{ij}\nabla
c_{j},$ (1.2)
$\displaystyle-\nabla\cdot\left(\epsilon\nabla\phi\right)=\rho_{0}+\sum_{i=0}^{N}z_{i}\text{e}c_{i},$
(1.3)

where $N$ is the number of ion species, $\mathcal{J}_{i}$ is the flux density,
$D_{i}$ is the diffusion constant for $i=0,1,\dots,N$. In addition, $c_{0}$ is
the concentration of solvent molecules with the valence $z_{0}=0$, the radius
$a_{0}$, and $c_{i}$ is the concentration of the $i$th ion species with the
valence $z_{i}\neq 0$ and the radius $a_{i}$ for $i=1,\dots,N$. Moreover,
$\phi$ is the electrostatic potential, $\epsilon$ is the dielectric function,
$\rho_{0}$ is the permanent charge density function, $k_{B}$ is the Boltzmann
constant, $T$ is the absolute temperature, e is the elementary charge and

$g_{ij}=S_{\sigma}\epsilon_{ij}(a_{i}+a_{j})^{12}\geq 0,$ (1.4)

where $g_{ij}$ represents the strength of steric repulsion, $\epsilon_{ij}$ is
the energy coupling constant, and $S_{\sigma}\sim\sigma^{d-12}$ for the
dimension $d\leq 3$ (cf. [35]).

To get the steady state of equations (1.1)–(1.3), we set $\mathcal{J}_{i}=0$
for $i=0,1,\dots,N$. Then equation (1.2) implies

$D_{i}\ln
c_{i}+\frac{D_{i}}{k_{B}T}z_{i}\text{e}\phi+\frac{D_{i}}{k_{B}T}\sum_{j=0}^{N}g_{ij}c_{j}=\mu_{i},\quad
i=0,1,\dots,N,$ (1.5)

where $\mu_{i}$ is constant for $i=0,1,\dots,N$. Without loss of generality,
we set $D_{i}=k_{B}T=\text{e}=1$ and equations (1.3) and (1.5) become

$\displaystyle\ln
c_{i}+z_{i}\phi+\sum_{j=0}^{N}g_{ij}c_{j}=\mu_{i}\quad\text{for}~{}i=0,1,\dots,N,$
(1.6)
$\displaystyle-\nabla\cdot\left(\epsilon\nabla\phi\right)=\rho_{0}+\sum_{i=1}^{N}z_{i}c_{i}\quad\text{in}~{}\Omega,$
(1.7)

where $\Omega$ is a bounded smooth domain in $\mathbb{R}^{d}$. Under the
assumptions of $g_{ij}$ and $\mu_{i}$ (see the assumptions of steric effects
(A1) and (A2) below), equation (1.6) has unique smooth solutions
$c_{i}=c_{i}(\phi)$, for $\phi\in\mathbb{R}$ and $i=0,1,\dots,N$, so equation
(1.7) becomes a single nonlinear elliptic equation (see equation (1.12)).
Notice that equation (1.6) may have multiple solutions if $g_{ij}$’s satisfy
the other (different from (A1)) conditions in [36]. Hereafter, the boundary
condition of equation (1.7) is either the Robin boundary condition

$\phi+\eta\frac{\partial\phi}{\partial\nu}=\phi_{bd}\quad\text{on}~{}\partial\Omega,$
(1.8)

or the Neumann boundary condition

$\frac{\partial\phi}{\partial\nu}=0\quad\text{on}~{}\partial\Omega,$ (1.9)

where $\phi_{bd}\in\mathcal{C}^{2}(\overline{\Omega})$ is the extra
electrostatic potential and $\eta\geq 0$ is a constant related to the surface
dielectric constant.

Volume exclusion is considered to develop the modified Poisson–Boltzmann (mPB)
equation with steric effects (cf. [10]). When volume exclusion occurs, ions
and solvent molecules are well separated. This motivates us to assume that
ions and solvent molecules have strong repulsion in order to separate ions and
solvent molecules. To describe strong repulsion of ions and solvent molecules,
we set $g_{ij}=\Lambda\tilde{g}_{ij}$ for $i,j=0,1,\dots,N$, where $\Lambda$
is a large parameter tending to infinity and each $\tilde{g}_{ij}$ is a
positive constant independent of $\Lambda$. The main goal of this paper is to
study the asymptotic limit of equations (1.6) and (1.7) as $\Lambda$ goes to
infinity. The singularity of matrix $\tilde{G}=[\tilde{g}_{ij}]$ plays a
crucial role on this problem. Suppose that $\tilde{G}$ is nonsingular with
inverse matrix $\tilde{G}^{-1}$ and ${\bf
c}_{\Lambda}=[c_{0,\Lambda},\dots,c_{N,\Lambda}]^{\mathsf{T}}$ is a solution
of equation (1.6), where $c_{i,\Lambda}$ is smooth function of $\phi$ for
$i=0,1,\dots,N$. Then we differentiate the equation (1.6) with respect to
$\phi$ and obtain $(D_{\Lambda}+\Lambda\tilde{G})\frac{\mathrm{d}{\bf
c}_{\Lambda}}{\mathrm{d}\phi}=-{\bf z}$, where $D_{\Lambda}$ is a diagonal
matrix with diagonal entries $(c_{i,\Lambda})^{-1}$ and ${\bf
z}=[z_{0},\cdots,z_{N}]^{\mathsf{T}}$. Suppose that $D_{\Lambda}$ is uniformly
bounded for $\Lambda\geq 1$, and $c_{i,\Lambda}$ converges to $c_{i}^{*}$ in
$\mathcal{C}^{1}(\mathbb{R})$ as $\Lambda$ tends to infinity. Then ${\bf
c}^{*}=[c_{0}^{*},\cdots,c_{N}^{*}]^{\mathsf{T}}$ satisfies

$\frac{\mathrm{d}{\bf
c}^{*}}{\mathrm{d}\phi}=-\lim_{\Lambda\to\infty}(D_{\Lambda}+\Lambda\tilde{G})^{-1}{\bf
z}=\lim_{\Lambda\to\infty}\frac{\mathrm{d}{\bf
c}_{\Lambda}}{\mathrm{d}\phi}=-\lim_{\Lambda\to\infty}\Lambda^{-1}\tilde{G}^{-1}(I+\Lambda^{-1}D_{\Lambda}\tilde{G}^{-1})^{-1}{\bf
z}={\bf 0},$

and hence, each $c_{i}^{*}$ becomes constant which is trivial. Therefore, when
$D_{\Lambda}$ is uniformly bounded for $\Lambda\geq 1$, matrix $\tilde{G}$
must be singular in order to have nontrivial limiting functions
$c_{i}^{*}=c_{i}^{*}(\phi)$ for $i=0,1,\dots,N$.

In the rest of this paper, we use the assumptions of steric effects given by

  * (A1)

$g_{ij}=\Lambda\lambda_{i}\lambda_{j}$ for $i,j=0,1,\dots,N$,

  * (A2)

$\mu_{i}=\Lambda\lambda_{i}\tilde{\mu}_{0}+\hat{\mu}_{i}$ for $i=0,1,\dots,N$,

where $\lambda_{i}$, $\tilde{\mu}_{0}$ and $\hat{\mu}_{i}$ are positive
constants independent of $\Lambda$. Note that (A1) assures that
$\tilde{G}=[\lambda_{i}\lambda_{j}]$ is a singular Gramian matrix which is
rank one. By (A1) and (A2), equation (1.6) has unique solutions
$c_{i,\Lambda}=c_{i,\Lambda}(\phi)$ satisfying

$\displaystyle
c_{i,\Lambda}=(c_{0,\Lambda})^{\lambda_{i}/\lambda_{0}}\exp(\bar{\mu}_{i}-z_{i}\phi)\quad\text{for}~{}\phi\in\mathbb{R}~{}\text{and}~{}i=1,\dots,N,$
(1.10) $\displaystyle\ln
c_{0,\Lambda}+\Lambda\lambda_{0}\sum_{j=0}^{N}\lambda_{j}(c_{0,\Lambda})^{\lambda_{j}/\lambda_{0}}\exp(\bar{\mu}_{j}-z_{j}\phi)=\mu_{0}\quad\text{for}~{}\phi\in\mathbb{R},$
(1.11)

where
$\bar{\mu}_{i}=\hat{\mu}_{i}-\frac{\lambda_{i}}{\lambda_{0}}\hat{\mu}_{0}$ is
a constant independent of $\Lambda$ for $i=0,1,\dots,N$. Then we apply the
implicit function theorem on (1.11) and obtain that $c_{0,\Lambda}(\phi)$ is a
positive smooth function (see Proposition 2.1). Hence equation (1.7) becomes a
nonlinear elliptic equation which is called a Poisson–Boltzmann equation with
steric effects (PB-steric equation) and is expressed as

$-\nabla\cdot(\epsilon\nabla\phi_{\Lambda})=\rho_{0}+f_{\Lambda}(\phi_{\Lambda})\quad\text{in}~{}\Omega,$
(1.12)

where function $f_{\Lambda}=f_{\Lambda}(\phi)$ is denoted as

$f_{\Lambda}(\phi)=\sum_{i=1}^{N}z_{i}c_{i,\Lambda}(\phi)\quad\text{for}~{}\phi\in\mathbb{R}.$
(1.13)

Notice that by (1.10), (1.11) and (1.13), functions $f_{\Lambda}$ is the
nonlinear term of PB-steric equation (1.12), and is mainly determined by the
Lambert type function $c_{0,\Lambda}$ (given by (1.11)) which represents the
concentration of solvent molecules. The main difficulty of this paper is to
establish the convergence of $c_{i,\Lambda}$ and $\phi_{\Lambda}$ as $\Lambda$
tends to infinity.

By the implicit function theorem on Banach spaces (cf. [13, Theorem 15.1]), we
prove that $c_{i,\Lambda}$ converges to $c_{i}^{*}$ in space
$\mathcal{C}^{m}\left[a,b\right]$ as $\Lambda$ tends to infinity for
$m\in\mathbb{N}$ and $a<b$ (see Propositions 2.5 and 2.6). Here function
$c_{i}^{*}$ satisfies

$\displaystyle
c_{i}^{*}(\phi)=(c_{0}^{*}(\phi))^{\lambda_{i}/\lambda_{0}}\exp\left(\bar{\mu}_{i}-z_{i}\phi\right)>0\quad\text{for}~{}\phi\in\mathbb{R}~{}\text{and}~{}i=1,\dots,N,$
(1.14)
$\displaystyle\sum_{i=0}^{N}\lambda_{i}c_{i}^{*}\left(\phi\right)=\sum_{i=0}^{N}\lambda_{i}(c_{0}^{*}(\phi))^{\lambda_{i}/\lambda_{0}}\exp(\bar{\mu}_{i}-z_{i}\phi)=\tilde{\mu}_{0}\quad\text{for}~{}\phi\in\mathbb{R},$
(1.15)

where
$\bar{\mu}_{i}=\hat{\mu}_{i}-\frac{\lambda_{i}}{\lambda_{0}}\hat{\mu}_{0}$.
Thus function $f_{\Lambda}$ also converges to $f^{*}$ in space
$\mathcal{C}^{m}\left[a,b\right]$ as $\Lambda$ goes to infinity for
$m\in\mathbb{N}$ and $a<b$, where

$f^{*}(\phi)=\sum_{i=1}^{N}z_{i}c_{i}^{*}(\phi)$ (1.16)

(see Figure 1). Note that for $\Lambda>0$, function $f_{\Lambda}$ is strictly
decreasing and unbounded on the entire space $\mathbb{R}$ (see Propositions
2.2 and 2.4), and so does function
$f_{PB}(\phi)=\sum_{i=1}^{N}z_{i}\exp(-z_{i}\phi)$ of the conventional
Poisson–Boltzmann (PB) equation
$-\nabla\cdot(\epsilon\nabla\phi)=\rho_{0}+f_{PB}(\phi)$ in $\Omega$. Besides,
function $f^{*}$ is also strictly decreasing but bounded on $\mathbb{R}$ (see
Propositions 2.7 and 2.8), and so does function $f_{mPB}(\phi)$ of the mPB
equation (1.18) in [9, 10, 28]. Obviously, $f_{\Lambda}$ cannot uniformly
converge to $f^{*}$ on $\mathbb{R}$ as $\Lambda$ goes to infinity. Hence we
only have uniform convergence of $f_{\Lambda}$ on any bounded interval
$\left[a,b\right]$ but not the entire space $\mathbb{R}$. To get the
asymptotic limit of equation (1.12), we firstly have to prove the uniform
boundedness of $\phi_{\Lambda}$ (the solution of equation (1.12)) with respect
to $\Lambda$ (see Lemmas 3.2 and 4.1) in order to use the convergence of
$f_{\Lambda}$ in space $\mathcal{C}^{m}[a,b]$ for $m\in\mathbb{N}$ and $a<b$.
Here equation (1.10) is crucial for the proof of Lemmas 3.2 and 4.1. Notice
that $\rho_{0}=\rho_{0}(x)$ may be any non-zero function and the boundary
condition of equation (1.12) may be the Robin or Neumann boundary condition
but not the Dirichlet boundary condition so one cannot simply use the maximum
principle on equation (1.12) to prove Lemmas 3.2 and 4.1. By Lemmas 3.2 and
4.1, we apply the $W^{2,p}$-estimate (cf. [2, Theorem 15.2]), the Schauder
estimate (cf. [22, Theorem 6.30]) and the uniqueness of solution of equation
(1.17) to prove that $\phi_{\Lambda}$ converges to $\phi^{*}$ in space
$\mathcal{C}^{m}(\Omega)$ as $\Lambda$ tends to infinity for $m\in\mathbb{N}$,
where $\phi^{*}$ satisfies

$-\nabla\cdot(\epsilon\nabla\phi^{*})=\rho_{0}+f^{*}(\phi^{*})\quad\text{in}~{}\Omega$
(1.17)

with the Robin and Neumann boundary conditions (1.8) and (1.9), respectively.

Figure 1: Curves $1$–$3$ are profiles of function $f_{\Lambda}$ (defined in
(1.13)) with $\Lambda=10,20,40$, and curve $4$ is the profile of function
$f^{*}$ (defined in (1.16)).

Now we state the following theorems.

######  Theorem 1.1.

Let $\Omega\subsetneq\mathbb{R}^{d}$ be a bounded smooth domain,
$\epsilon\in\mathcal{C}^{\infty}(\overline{\Omega})$ be a positive function,
$\rho_{0}\in\mathcal{C}^{\infty}(\overline{\Omega})$,
$\phi_{bd}\in\mathcal{C}^{2}(\partial\Omega)$, $z_{0}=0$, and $z_{i}z_{j}<0$
for some $i,j\in\\{1,\dots,N\\}$. Assume that (A1) and (A2) hold true. Then
the PB-steric equation (1.12) with the Robin boundary condition (1.8)
satisfies

$\lim_{\Lambda\to\infty}\left\|\phi_{\Lambda}-\phi^{*}\right\|_{\mathcal{C}^{m}(\Omega)}=0\quad\text{for}~{}m\in\mathbb{N},$

where $\phi^{*}$ is the solution of equation (1.17) with the Robin boundary
condition (1.8).

######  Theorem 1.2.

Let $\Omega\subsetneq\mathbb{R}^{d}$ be a bounded smooth domain,
$\epsilon\in\mathcal{C}^{\infty}(\overline{\Omega})$ be a positive function,
$z_{0}=0$, and $z_{i}z_{j}<0$ for some $i,j\in\\{1,\dots,N\\}$. Assume that
(A1) and (A2) hold true. Then

  1. (i)

Suppose that $\rho_{0}$ is constant. Then the PB-steric equation (1.12) with
the Neumann boundary condition (1.9) has a unique constant solution
$\phi_{\Lambda}$. Moreover if $\rho_{0}\in(-M^{*},-m^{*})$, then
$\lim_{\Lambda\to\infty}\phi_{\Lambda}=\phi^{*}$, where $\phi^{*}$ is the
unique constant solution of equation (1.17) with the Neumann boundary
condition (1.9), and

$M^{*}=\lim_{\phi\to-\infty}f^{*}\left(\phi\right)>0,\quad
m^{*}=\lim_{\phi\to\infty}f^{*}\left(\phi\right)<0.$

  2. (ii)

Suppose that $\rho_{0}\in\mathcal{C}^{\infty}(\overline{\Omega})$ is a
nonconstant function with $\int_{\Omega}\rho_{0}\,\mathrm{d}x=0$. Then the PB-
steric equation (1.12) with the Neumann boundary condition (1.9) has a unique
solution $\phi_{\Lambda}$ satisfying

$\lim_{\Lambda\to\infty}\|\phi_{\Lambda}-\phi^{*}\|_{\mathcal{C}^{m}(\Omega)}=0\quad\text{for}~{}m\in\mathbb{N},$

where $\phi^{*}$ is the unique solution of equation (1.17) with the Neumann
boundary condition (1.9).

######  Remark 1.

As $\Lambda=0$, equations (1.10) and (1.11) imply $c_{0,0}=e^{\hat{\mu}_{0}}$
and $c_{i,0}=\exp(\hat{\mu}_{i}-z_{i}\phi)$ for $i=1,\dots,N$. Here we replace
$\phi_{0}$ by $\phi$ so the PB-steric equation (1.12) becomes the conventional
Poisson–Boltzmann (PB) equation
$-\nabla\cdot(\epsilon\nabla\phi)=\rho_{0}+\sum_{i=1}^{N}z_{i}\exp(\hat{\mu}_{i}-z_{i}\phi)$
in $\Omega$. However, as $\Lambda>0$, a larger $\Lambda$ makes the steric
repulsion (between ions and solvent molecules) stronger, and eventually (as
$\Lambda\to\infty)$ the limiting equation of PB-steric equation (1.12) becomes
equation (1.17) which has the same form (up to scalar multiples) as the
modified PB (mPB) equations in [9, 10, 28, 32, 33, 34, 37] (see below). Hence
the PB-steric equation (1.12) is a new PB equation (with the steric effects of
ions and solvent molecules) which can be regarded as a generalized model of
mPB equations.

To show equation (1.17) with the same form (up to scalar multiples) as the mPB
equations (1.18) and (1.21), we use (1.4) and (A1) with $\Lambda=S_{\sigma}$
to obtain $\lambda_{i}\lambda_{j}=\epsilon_{ij}(a_{i}+a_{j})^{12}$ for
$i,j=0,1,\dots,N$ which implies
$\lambda_{i}=(\epsilon_{ii}(2a_{i})^{12})^{1/2}$ and
$\epsilon_{ij}=(\epsilon_{ii}\epsilon_{jj})^{1/2}\frac{2^{12}(a_{i}a_{j})^{6}}{(a_{i}+a_{j})^{12}}$
for $i,j=0,1,\dots,N$. Note that because of $S_{\sigma}\sim\sigma^{d-12}$ and
$1\leq d\leq 3$, $S_{\sigma}\to\infty$ iff $\sigma\to 0^{+}$, where $\sigma>0$
is the parameter of the approximate Lennard–Jones potential (cf. [35]). Hence
the nonlinear term of equation (1.17) is function $f^{*}$ which depends on
$a_{i}$’s the radii of ions and solvent molecules. Suppose that
$\lambda_{i}=a_{i}^{3}$, i.e., $\epsilon_{ii}=(2^{12}a_{i}^{6})^{-1}$ for
$i=0,1,\dots,N$. If ions and solvent molecules have the same size, i.e.,
$a_{i}=a_{0}$ for all $i$, then equations (1.14) and (1.15) imply

$c_{i}^{*}(\phi)=\frac{\tilde{\mu}_{0}}{a_{0}^{3}}\frac{\exp(\bar{\mu}_{i}-z_{i}\phi)}{\sum_{i=0}^{N}\exp(\bar{\mu}_{i}-z_{i}\phi)}\quad\text{for}~{}\phi\in\mathbb{R}~{}\text{and}~{}i=0,1,\dots,N,$

which can be plugged in equations (1.16) and (1.17) (with $\phi=\phi^{*}$) to
get the following mPB equation (same size of ions and solvent molecules)

$-\nabla\cdot(\epsilon\nabla\phi)=\rho_{0}+f_{mPB}(\phi)\quad\text{in}~{}\Omega,$
(1.18)

where

$f_{mPB}(\phi)=\frac{\tilde{\mu}_{0}}{a_{0}^{3}}\frac{\sum_{i=1}^{N}z_{i}\exp(\bar{\mu}_{i}-z_{i}\phi)}{\sum_{i=0}^{N}\exp(\bar{\mu}_{i}-z_{i}\phi)}\quad\text{for}~{}\phi\in\mathbb{R},$

see [9, 10, 28]. Conversely, if ions and solvent molecules do not have the
same size, i.e., there exists $a_{i}\neq a_{0}$ for some
$i\in\\{1,\dots,N\\}$, then equations (1.14)–(1.17) can be transformed to the
following mPB equations

$\displaystyle\sum_{i=0}^{N}v_{i}c_{i}=1,$ (1.19)
$\displaystyle\frac{v_{i}}{v_{0}}\ln(v_{0}c_{0})-\ln(v_{i}c_{i})=\beta
q_{i}\psi-\beta\mu_{i},$ (1.20)
$\displaystyle-\nabla\cdot(\epsilon\nabla\phi)=\rho_{0}+\sum_{i=0}^{N}\beta
q_{i}v_{i}c_{i}(\psi)~{}~{}\text{in}~{}\Omega,$ (1.21)

see [32, 33, 34, 37]. Here we set
$\lambda_{i}=a_{i}^{3}=v_{i}\tilde{\mu}_{0}$, $c_{i}^{*}=c_{i}$,
$\bar{\mu}_{i}=\beta\mu_{i}+\frac{v_{i}}{v_{0}}\ln v_{0}-\ln v_{i}$,
$z_{i}=\beta q_{i}$ and $\phi^{*}=\phi=\psi$. Moreover, equations
(1.14)–(1.17) have the same form (up to scalar multiples) as equations
(1.19)–(1.21) (see Appendix II). Therefore, as for Remark 1, the PB-steric
equation (equation (1.12) with the nonlinear term given by equations (1.10),
(1.11) and (1.13)) can be regarded as a generalized model of mPB equations.

In Section 5, we provide numerical simulations on equations (1.12) and (1.17)
with the Robin and Neumann boundary conditions (1.8) and (1.9), respectively.
The one-dimensional domain $\Omega=\left(-1,1\right)$ is discretized by the
Legendre–Gauss–Lobatto (LGL) points (cf. [20]). Then we use the LGL points to
discretize equations (1.12) and (1.17), and solve them numerically with the
command fsolve in Matlab. Under the Robin boundary condition (1.8), we show
the profiles of solutions $\phi_{\Lambda}$ and $\phi^{*}$ to support Theorem
1.1. Moreover, we demonstrate the profiles of solutions $\phi_{\Lambda}$ and
$\phi^{*}$ under the Neumann boundary condition (1.9) to support Theorem 1.2.

Organization. The rest of the paper is organized as follows. In Section 2, the
analysis of functions $f_{\Lambda}$ and $f^{*}$ is presented. We prove
Theorems 1.1 and 1.2 in Sections 3 and 4, respectively. Numerical results are
shown in Section 5.

##  2 Analysis of $f_{\Lambda}$ and $f^{*}$

In this section, we show how to solve equation (1.6) by (A1) and (A2) and
obtain smooth solutions $c_{i,\Lambda}=c_{i,\Lambda}(\phi)$ for
$\phi\in\mathbb{R}$ and $i=0,1,\dots,N$. Note that because of (A1), equation
(1.6) can be expressed as

$\ln
c_{i,\Lambda}+z_{i}\phi+\Lambda\sum_{j=0}^{N}\lambda_{i}\lambda_{j}c_{j,\Lambda}=\mu_{j}\quad\text{for}~{}\phi\in\mathbb{R}~{}\text{and}~{}i=0,1,\dots,N.$
(2.1)

By (A2) and Gaussian elimination, we solve equation (2.1) and get (1.10). Then
we plug (1.10) into (2.1) with $i=0$ and obtain equation (1.11). Moreover, we
apply the implicit function theorem to prove that equation (1.11) has a smooth
solution $c_{0,\Lambda}=c_{0,\Lambda}(\phi)$ for $\phi\in\mathbb{R}$ (see
Proposition 2.1). As $\Lambda$ goes to infinity, we use the implicit function
theorem on Banach space to show that $c_{i,\Lambda}$ converges to $c_{i}^{*}$
in space $\mathcal{C}^{m}[a,b]$ for any $m\in\mathbb{N}$ and $a<b$ (see
Proposition 2.6). Hence $f_{\Lambda}=\sum_{i=1}^{N}z_{i}c_{i,\Lambda}(\phi)$
also converges to $f^{*}=\sum_{i=1}^{N}z_{i}c_{i}^{*}(\phi)$ in space
$\mathcal{C}^{m}[a,b]$.

###  2.1 Analysis of $c_{i,\Lambda}$ and $f_{\Lambda}$

For functions $c_{i,\Lambda}$ and $f_{\Lambda}$, we have the following
propositions.

######  Proposition 2.1.

Assume that $z_{0}=0$, (A1) and (A2) hold true. Then equation (2.1) has smooth
and positive solutions $c_{i,\Lambda}=c_{i,\Lambda}(\phi)$ for
$\phi\in\mathbb{R}$ and $i=0,1,\cdots,N$.

###### Proof.

Multiplying the equation (2.1) for $i=0$ by $\lambda_{i}/\lambda_{0}$, we
obtain

$\frac{\lambda_{i}}{\lambda_{0}}\ln
c_{0,\Lambda}+\Lambda\sum_{j=0}^{N}\lambda_{i}\lambda_{j}c_{j,\Lambda}=\frac{\lambda_{i}}{\lambda_{0}}\mu_{0}.$
(2.2)

Then we subtract the equation (2.1) from the equation (2.2) to get

$\ln c_{i,\Lambda}-\frac{\lambda_{i}}{\lambda_{0}}\ln
c_{0,\Lambda}+z_{i}\phi=\mu_{i}-\frac{\lambda_{i}}{\lambda_{0}}\mu_{0},$

which implies

$c_{i,\Lambda}=c_{0,\Lambda}^{\lambda_{i}/\lambda_{0}}\exp\left(\bar{\mu}_{i}-z_{i}\phi\right)\quad\text{for}~{}\phi\in\mathbb{R}~{}\text{and}~{}i=1,\dots,N,$
(2.3)

where
$\bar{\mu_{i}}=\mu_{i}-\frac{\lambda_{i}}{\lambda_{0}}\mu_{0}=\hat{\mu}_{i}-\frac{\lambda_{i}}{\lambda_{0}}\hat{\mu}_{0}$
because of (A2). Note that $\bar{\mu}_{i}$ is independent of $\Lambda$.
Plugging (2.3) into (2.2), we have

$\ln
c_{0,\Lambda}+\Lambda\lambda_{0}\sum_{j=0}^{N}\,\lambda_{j}\,c_{0,\Lambda}^{\lambda_{j}/\lambda_{0}}\,\exp(\bar{\mu}_{j}-z_{j}\phi)-\mu_{0}=0\quad\text{for}~{}\phi\in\mathbb{R},$
(2.4)

which can be denoted as $h_{1}(c_{0,\Lambda},\phi)=0$. Here $h_{1}$ is defined
by

$h_{1}\left(t,\phi\right)=\ln
t+\Lambda\lambda_{0}\sum_{j=0}^{N}\lambda_{j}t^{\lambda_{j}/\lambda_{0}}\exp\left(\bar{\mu}_{j}-z_{j}\phi\right)-\mu_{0}\quad\text{for}~{}t>0~{}\text{and}~{}\phi\in\mathbb{R}.$

Notice that, for any $\phi\in\mathbb{R}$, $h_{1}$ is strictly increasing for
$t>0$ and the range of $h_{1}$ is entire space $\mathbb{R}$. Then there exists
a unique positive number $c_{0,\Lambda}(\phi)$ such that
$h_{1}\left(c_{0,\Lambda}(\phi),\phi\right)=0$ for $\phi\in\mathbb{R}$.
Moreover, since $h_{1}$ is smooth for $t>0$, $\phi\in\mathbb{R}$, and

$\frac{\partial h_{1}}{\partial
t}\left(t,\phi\right)=\frac{1}{t}+\Lambda\sum_{j=0}^{N}\lambda_{j}^{2}t^{(\lambda_{j}-\lambda_{0})/\lambda_{0}}\exp\left(\bar{\mu}_{j}-z_{j}\phi\right)>0\quad\text{for}~{}t>0~{}\text{and}~{}\phi\in\mathbb{R},$

then by the implicit function theorem (cf. [31, Theorem 3.3.1]),
$c_{0,\Lambda}=c_{0,\Lambda}(\phi)$ is a smooth and positive function on
$\mathbb{R}$. Therefore, by (2.3), each $c_{i,\Lambda}$ is smooth and positive
on $\mathbb{R}$ and we complete the proof of Proposition 2.1. ∎

######  Proposition 2.2.

Function $\displaystyle
f_{\Lambda}(\phi)=\sum_{i=0}^{N}z_{i}c_{i,\Lambda}(\phi)$ is strictly
decreasing on $\mathbb{R}$, where functions $c_{i,\Lambda}$ are obtained in
Proposition 2.1.

###### Proof.

By Proposition 2.1, we differentiate the equation (2.1) with respect to $\phi$
and obtain

$(D_{\Lambda}+G)\frac{\mathrm{d}{\bf c}_{\Lambda}}{\mathrm{d}\phi}=-{\bf
z}\quad\text{for}~{}\phi\in\mathbb{R},$ (2.5)

where matrix $D_{\Lambda}=\text{diag}(1/c_{0,\Lambda},\cdots,1/c_{N,\Lambda})$
is positive definite, matrix
$G=\Lambda[\lambda_{0}~{}\cdots~{}\lambda_{N}]^{\mathsf{T}}[\lambda_{0}~{}\cdots~{}\lambda_{N}]$,
${\bf c}_{\Lambda}=[c_{0,\Lambda}~{}\cdots~{}c_{N,\Lambda}]^{\mathsf{T}}$ and
${\bf z}=[z_{0}~{}\cdots~{}z_{N}]^{\mathsf{T}}$. It is obvious that matrix
$D_{\Lambda}+G$ is positive definite and invertible with inverse matrix
$(D_{\Lambda}+G)^{-1}$ which is also positive definite. Then equation (2.5)
gives $\mathrm{d}{\bf c}_{\Lambda}/\mathrm{d}\phi=-(D_{\Lambda}+G)^{-1}{\bf
z}$, and $z^{\mathsf{T}}(\mathrm{d}{\bf c}_{\Lambda}/\mathrm{d}\phi)$ becomes

$\frac{\mathrm{d}f_{\Lambda}}{\mathrm{d}\phi}=\sum_{i=1}^{N}z_{i}\frac{\mathrm{d}c_{i,\Lambda}}{\mathrm{d}\phi}={\bf
z}^{\mathsf{T}}\frac{\mathrm{d}{\bf c}}{\mathrm{d}\phi}=-{\bf
z}^{\mathsf{T}}\left(D+G\right)^{-1}{\bf
z}<0\quad\text{for}~{}\phi\in\mathbb{R}.$

Here we have used the fact that ${\bf z}\neq 0$ and we complete the proof of
Proposition 2.2. ∎

######  Proposition 2.3.

Assume that (A1) and (A2) hold true.

  1. (i)

If $z_{k}\geq 0$ ($z_{k}\leq 0$) and $k\in\\{0,1,\dots,N\\}$, then
$\displaystyle\sup_{\phi\geq 0}c_{k,\Lambda}(\phi)\leq e^{\mu_{k}}$
($\displaystyle\sup_{\phi\leq 0}c_{k,\Lambda}(\phi)\leq e^{\mu_{k}}$).

  2. (ii)

Suppose $z_{i}z_{j}<0$ for some $i,j\in\\{1,\dots,N\\}$. Then there exist
$i_{0},j_{0}\in\\{1,\dots,N\\}$, $z_{i_{0}}>0$ and $z_{j_{0}}<0$ such that
$\displaystyle\limsup_{\phi\to-\infty}c_{i_{0},\Lambda}(\phi)=\infty$ and
$\displaystyle\limsup_{\phi\to\infty}c_{j_{0},\Lambda}(\phi)=\infty$.

###### Proof.

Suppose that $z_{k}\geq 0$ for some $k\in\\{0,\dots,N\\}$. Then equation (2.1)
implies

$\sup_{\phi\geq 0}\ln c_{k,\Lambda}(\phi)=\sup_{\phi\geq
0}\left(\mu_{k}-z_{k}\phi-\sum_{j=0}^{N}g_{kj}c_{j,\Lambda}(\phi)\right)\leq\mu_{k}<\infty,$

and $\displaystyle\sup_{\phi\geq 0}c_{k,\Lambda}(\phi)\leq e^{\mu_{k}}$. Here
we have used the fact that $c_{i,\Lambda}(\phi)>0$ for $\phi\in\mathbb{R}$ and
$i=0,1,\cdots,N$. Similarly, if $z_{k}\leq 0$ and $k\in\\{0,\cdots,N\\}$, then
equation (2.1) gives

$\sup_{\phi\leq 0}\ln c_{k,\Lambda}(\phi)=\sup_{\phi\leq
0}\left(\mu_{k}-z_{k}\phi-\Lambda\lambda_{k}\sum_{j=0}^{N}\lambda_{j}c_{j,\Lambda}(\phi)\right)\leq\mu_{k}<\infty,$

and $\displaystyle\sup_{\phi\leq 0}c_{k,\Lambda}(\phi)\leq e^{\mu_{k}}$. Hence
the proof of Proposition 2.3(i) is complete.

It remains to prove (ii). Since $z_{i}z_{j}<0$ for some
$i,j\in\\{1,\dots,N\\}$, then box index sets $I=\\{i:z_{i}>0\\}$ and
$J=\\{j:z_{j}<0\\}$ are nonempty. Now we claim that there exists $i_{0}\in I$
such that
$\displaystyle\limsup_{\phi\to-\infty}c_{i_{0},\Lambda}(\phi)=\infty$. We
prove this by contradiction. Suppose $\displaystyle\sup_{\phi\leq
0}c_{i,\Lambda}(\phi)<\infty$ for all $i\in I$. Then there exists $K_{1}>0$
such that $0<c_{i,\Lambda}(\phi)<K_{1}$ for $\phi\leq 0$ and $i\in I$. By
equation (2.1) and Proposition (2.3)(i), we have

$\displaystyle z_{i}\phi$ $\displaystyle=\mu_{i}-\ln
c_{i,\Lambda}(\phi)-\Lambda\lambda_{i}\sum_{k=0}^{N}\lambda_{j}c_{k,\Lambda}(\phi)$
$\displaystyle=\mu_{i}-\ln
c_{i,\Lambda}(\phi)-\Lambda\lambda_{i}\left(\sum_{k\in
I}\lambda_{k}c_{k,\Lambda}(\phi)+\sum_{j\in\\{0\\}\cup
J}\lambda_{k}c_{k,\Lambda}(\phi)\right)$ $\displaystyle\geq\mu_{i}-\ln
K_{1}-\Lambda\lambda_{i}\left(\sum_{k\in
I}\lambda_{k}K_{1}+\sum_{k\in\\{0\\}\cup
J}\lambda_{k}e^{\mu_{k}}\right)\quad\text{for}~{}i\in
I~{}\text{and}~{}\phi\leq 0,$

which leads a contradiction by letting $\phi\to-\infty$. Hence there exists
$i_{0}\in I$ such that
$\displaystyle\limsup_{\phi\to-\infty}c_{i_{0},\Lambda}(\phi)=\infty$.
Similarly, we claim that there exists $j_{0}\in J$ such that
$\displaystyle\limsup_{\phi\to\infty}c_{j_{0},\Lambda}(\phi)=\infty$. We also
prove this by contradiction. Suppose $\displaystyle\sup_{\phi\geq
0}c_{j,\Lambda}(\phi)<\infty$ for all $j\in J$. Then there exists $K_{2}>0$
such that $0<c_{j,\Lambda}(\phi)\leq K_{2}$ for $\phi\geq 0$ and $j\in J$. By
equation (2.1) and Proposition 2.3(i), we have

$\displaystyle z_{j}\phi$ $\displaystyle=\mu_{j}-\ln
c_{j,\Lambda}(\phi)-\Lambda\lambda_{j}\sum_{k=0}^{N}\lambda_{k}c_{k,\Lambda}(\phi)$
$\displaystyle=\mu_{j}-\ln
c_{j,\Lambda}(\phi)-\Lambda\lambda_{j}\left(\sum_{k\in\\{0\\}\cup
I}\lambda_{k}c_{k,\Lambda}(\phi)+\sum_{k\in
J}\lambda_{k}c_{k,\Lambda}(\phi)\right)$ $\displaystyle\geq\mu_{j}-\ln
K_{2}-\Lambda\lambda_{j}\left(\sum_{k\in\\{0\\}\cup
I}\lambda_{k}e^{\mu_{k}}+\sum_{k\in
J}\lambda_{k}K_{2}\right)\quad\text{for}~{}k\in J~{}\text{and}~{}\phi\geq 0,$

which leads a contradiction by letting $\phi\to\infty$. Therefore, the proof
of Proposition 2.3(ii) is complete.∎

######  Proposition 2.4.

For each $\Lambda>0$, the range of $f_{\Lambda}$ is entire space $\mathbb{R}$,
and
$\displaystyle\lim_{\phi\to\pm\infty}f_{\Lambda}\left(\phi\right)=\mp\infty$.

###### Proof.

By (1.13), function $f_{\Lambda}$ can be denoted as

$f_{\Lambda}(\phi)=\sum_{z_{i}>0}z_{i}c_{i,\Lambda}(\phi)+\sum_{z_{j}<0}z_{j}c_{j,\Lambda}(\phi)\quad\text{for}~{}\phi\in\mathbb{R}.$
(2.6)

Proposition 2.3(i) gives

$\sum_{z_{i}>0}z_{i}c_{i,\Lambda}(\phi)<\sum_{z_{i}>0}z_{i}e^{\mu_{i}}\quad\text{for}~{}\phi\geq
0.$ (2.7)

Moreover, by Proposition 2.3(ii), there exists $j_{0}\in\\{1,\dots,N\\}$ with
$z_{j_{0}}<0$ and a sequence $\\{\phi_{n}\\}_{n=1}^{\infty}$ with
$\displaystyle\lim_{n\to\infty}\phi_{n}=\infty$ such that
$\displaystyle\lim_{n\to\infty}z_{j_{0}}c_{j_{0},\Lambda}(\phi_{n})=-\infty$.
Thus, by (2.6) and (2.7), we have

$f_{\Lambda}(\phi_{n})\leq\sum_{z_{i}>0}z_{i}e^{\mu_{i}}+z_{j_{0}}c_{j_{0},\Lambda}(\phi_{n})\to-\infty\quad\text{as}~{}n\to\infty,$

which implies that
$\displaystyle\lim_{\phi\to\infty}f_{\Lambda}(\phi)=-\infty$ because of the
monotone decreasing of $f_{\Lambda}$ (see Proposition 2.2). On the other hand,
Proposition 2.3(i), also implies

$\sum_{z_{i}<0}z_{i}c_{i,\Lambda}(\phi)\geq\sum_{z_{i}<0}z_{i}e^{\mu_{i}}\quad\text{for}~{}\phi\leq
0.$ (2.8)

Moreover, by Proposition 2.3(ii), there exist $i_{0}\in\\{1,\dots,N\\}$ with
$z_{i_{0}}>0$ and a sequence $\\{\tilde{\phi}_{n}\\}_{n=1}^{\infty}$ with
$\displaystyle\lim_{n\to\infty}\tilde{\phi}_{n}=-\infty$ such that
$\displaystyle\lim_{n\to\infty}z_{i_{0}}c_{i_{0},\Lambda}(\tilde{\phi}_{n})=\infty$.
Thus, by (2.6) and (2.8), we obtain

$f_{\Lambda}(\tilde{\phi}_{n})\geq
z_{i_{0}}c_{i_{0},\Lambda}(\tilde{\phi}_{n})+\sum_{z_{i}<0}z_{i}e^{\mu_{i}}\to\infty\quad\text{as}~{}n\to\infty,$

which implies that
$\displaystyle\lim_{\phi\to-\infty}f_{\Lambda}(\phi)=\infty$ due to the
monotone decreasing of $f_{\Lambda}$. Therefore, the proof of Proposition 2.4
is complete. ∎

###  2.2 Analysis of $c_{i}^{*}$ and $f^{*}$

Function $c_{0}^{*}$ is the limit
$\displaystyle\lim_{\Lambda\to\infty}\,c_{0,\Lambda}$ (see Proposition 2.6),
where function $c_{0,\Lambda}$ is the solution of equation (2.4) for
$\Lambda>0$. Let $\delta=\Lambda^{-1}$ and
$\tilde{c}_{0,\delta}=c_{0,\Lambda}$. Then by (A2) and equation (2.4),
function $\tilde{c}_{0,\delta}$ satisfies

$\delta\ln\tilde{c}_{0,\delta}(\phi)+\lambda_{0}\sum_{i=0}^{N}\,\lambda_{i}\,(\tilde{c}_{0,\delta}(\phi))^{\lambda_{i}/\lambda_{0}}\,\exp(\bar{\mu}_{i}-z_{i}\phi)=\lambda_{0}\tilde{\mu}_{0}+\delta\hat{\mu_{0}}\quad\text{for}~{}\phi\in\mathbb{R}.$
(2.9)

Notice that $\Lambda\to\infty$ is equivalent to $\delta\to 0^{+}$ so
$c_{0}^{*}$ is also equal to the limit $\displaystyle\lim_{\delta\to
0^{+}}\,\tilde{c}_{0,\delta}$. Moreover, $c_{0}^{*}$ satisfies equation (1.15)
which is equation (2.9) with $\delta=0$. Besides, by (2.3), functions
$c_{i}^{*}$ satisfy the equation (1.14) for $i=1,\dots,N$. The existence and
uniqueness of equations (1.14) and (1.15) is proved in Proposition 2.5. The
convergence of $\tilde{c}_{0,\delta}$ as $\delta\to 0^{+}$, i.e. the
convergence of $c_{0,\Lambda}$ as $\Lambda\to\infty$ is proved in Proposition
2.6 so by (2.3), we obtain the convergence of $c_{i,\Lambda}$ as
$\Lambda\to\infty$.

Now we state Propositions 2.5 and 2.6 as follows.

######  Proposition 2.5.

Equations (1.14) and (1.15) have a unique solution
$(c_{0}^{*},\dots,c_{N}^{*})$ and each function $c_{i}^{*}=c_{i}^{*}(\phi)$ is
a smooth and positive function for $\phi\in\mathbb{R}$ and $i=0,1,\dots,N$.

###### Proof.

Equations (1.14) and (1.15) can be solved by the following problem.

$\sum_{j=0}^{N}\,\lambda_{j}(c_{0}^{*}(\phi))^{\lambda_{j}/\lambda_{0}}\,\exp(\bar{\mu}_{j}-z_{j}\phi)-\tilde{\mu}_{0}=0\quad\text{for}~{}\phi\in\mathbb{R},$

which can be represented by $h_{2}(c_{0}^{*}(\phi),\phi)=0$ for
$\phi\in\mathbb{R}$. Here $h_{2}$ is defined by

$h_{2}(t,\phi)=\sum_{j=0}^{N}\lambda_{j}t^{\lambda_{j}/\lambda_{0}}\exp(\bar{\mu}_{j}-z_{j}\phi)-\tilde{\mu}_{0}\quad\text{for}~{}t>0~{}\text{and}~{}\phi\in\mathbb{R}.$

Note that $\displaystyle\lim_{t\to 0^{+}}h_{2}(t,\phi)=-\tilde{\mu}_{0}<0$ and
$\displaystyle\lim_{t\to\infty}h_{2}(t,\phi)=\infty$ for any
$\phi\in\mathbb{R}$. Thus, for any $\phi\in\mathbb{R}$, there exists $t>0$
such that $h_{2}(t,\phi)=0$. Since $\frac{\partial h_{2}}{\partial
t}=\frac{1}{\lambda_{0}}\sum_{j=0}^{N}\lambda_{j}^{2}t^{(\lambda_{j}-\lambda_{0})/\lambda_{0}}\exp(\bar{\mu}_{j}-z_{j}\phi)>0$
for $t>0$ and $\phi\in\mathbb{R}$, then we may use the implicit function
theorem (cf. [31, Theorem 3.3.1]) to conclude that there exists a unique
smooth and positive function $c_{0}^{*}:\mathbb{R}\to\mathbb{R}^{+}$ such that
$h_{2}(c_{0}^{*}(\phi),\phi)=0$ for $\phi\in\mathbb{R}$. Therefore by equation
(1.14), we obtain the smooth and positive functions $c_{i}^{*}$ for
$i=1,\dots,N$, and complete the proof of Proposition 2.5. ∎

######  Proposition 2.6.

$\displaystyle\lim_{\Lambda\to\infty}\|c_{i,\Lambda}-c_{i}^{*}\|_{\mathcal{C}^{m}\left[a,b\right]}=0$
for $i=0,1,\dots,N$, $a<b$ and $m\in\mathbb{N}$, where
$\displaystyle\|h\|_{\mathcal{C}^{m}\left[a,b\right]}:=\sum_{k=0}^{m}\|h^{(k)}\|_{\infty}$
for $h\in\mathcal{C}^{m}\left[a,b\right]$.

###### Proof.

Fix $m\in\mathbb{N}$, $a,b\in\mathbb{R}$ and $a<b$ arbitrarily. Let
$\|\cdot\|_{\mathcal{C}^{m}}:=\|\cdot\|_{\mathcal{C}^{m}[a,b]}$ for notation
convenience. For $\Lambda>0$, let $\delta=\Lambda^{-1}$,
$\tilde{c}_{0,\delta}=c_{0,\Lambda}$, $w_{\delta}=\ln\tilde{c}_{0,\delta}$ and
$w^{*}=\ln c_{0}^{*}$. Obviously, $\delta\to 0^{+}$ is equivalent to
$\Lambda\to\infty$. Hence by (2.3), it suffices to show that
$\displaystyle\lim_{\delta\to
0^{+}}\|e^{w_{\delta}}-e^{w^{*}}\|_{\mathcal{C}^{m}}=0$. Because
$w_{\delta}=\ln\tilde{c}_{0,\delta}$ and $\tilde{c}_{0,\delta}=c_{0,\Lambda}$,
equation (2.9) can be denoted as

$H(w_{\delta}(\phi),\delta)=0$ (2.10)

for $\delta>0$ and $\phi\in[a,b]$, where $H$ is a $\mathcal{C}^{1}$-function
on $\mathcal{C}^{m}[a,b]\times\mathbb{R}$ defined by

$H(w(\phi),\delta)=-\tilde{\mu}_{0}+\sum_{i=0}^{N}\lambda_{i}\exp\left(\frac{\lambda_{i}}{\lambda_{0}}w(\phi)+\bar{\mu}_{i}-z_{i}\phi\right)-\frac{\hat{\mu}_{0}-w(\phi)}{\lambda_{0}}\cdot\delta$
(2.11)

for all $w\in\mathcal{C}^{m}[a,b]$ and $\phi\in[a,b]$. Note that
$H(w^{*},0)=0$ by Proposition 2.5. A direct calculation for Fréchet derivative
of (2.11) gives

$D_{w}H(w(\phi),\delta)=\frac{1}{\lambda_{0}}\sum_{i=0}^{N}\lambda_{i}^{2}\exp\left(\frac{\lambda_{i}}{\lambda_{0}}w(\phi)+\bar{\mu}_{i}-z_{i}\phi\right)+\frac{\delta}{\lambda_{0}}$

for all $w\in\mathcal{C}^{m}[a,b]$ and $\phi\in[a,b]$. Then
$D_{w}H(w^{*}(\phi),0)=\frac{1}{\lambda_{0}}\sum_{i=0}^{N}\lambda_{i}^{2}(c_{0}^{*}(\phi))^{\lambda_{i}/\lambda_{0}}\exp(\bar{\mu}_{i}-z_{i}\phi)>0$
for $\phi\in[a,b]$. Here we have used the fact that $w^{*}=\ln c_{0}^{*}$.
This implies that $D_{w}H(w^{*},0)I$ is a bounded and invertible linear map on
the Banach space $\mathcal{C}^{m}[a,b]$, where $I$ is an identity map. Hence
by the implicit function theorem on Banach spaces (cf. [13, Corollary 15.1]),
there exist an open subset
$B_{\delta_{0}}(w^{*})\times(-\delta_{0},\delta_{0})\subsetneq\mathcal{C}^{m}[a,b]\times\mathbb{R}$
and a unique $\mathcal{C}^{1}$-function $\tilde{w}(\cdot,\delta)$ of
$\delta\in(-\delta_{0},\delta_{0})$ with $\tilde{w}(\cdot,\delta)\in
B_{\delta_{0}}(w^{*})\subsetneq\mathcal{C}^{m}[a,b]$ for
$\delta\in(-\delta_{0},\delta_{0})$ such that

$H(\tilde{w}(\cdot,\delta),\delta)=0~{}~{}\text{for
all}~{}\delta\in(-\delta_{0},\delta_{0}),$

which gives $\displaystyle\lim_{\delta\to
0^{+}}\|\tilde{w}(\cdot,\delta)-w^{*}\|_{\mathcal{C}^{m}}=0$. By Proposition
2.1, equation (2.10) has a unique solution $w_{\delta}$, which implies
$\tilde{w}(\cdot,\delta)=w_{\delta}(\cdot)$ for
$\delta\in(-\delta_{0},\delta_{0})$. Therefore, we obtain
$\displaystyle\lim_{\delta\to
0^{+}}\left\|w_{\delta}-w^{*}\right\|_{\mathcal{C}^{m}}=0$, i.e.,
$\displaystyle\lim_{\delta\to
0^{+}}\|e^{w_{\delta}}-e^{w^{*}}\|_{\mathcal{C}^{m}}=0$, and complete the
proof of Proposition 2.6. ∎

######  Remark 2.

Because $\displaystyle
f_{\Lambda}(\phi)=\sum_{i=1}^{N}z_{i}c_{i,\Lambda}(\phi)$ and $\displaystyle
f^{*}(\phi)=\sum_{i=1}^{N}z_{i}c_{i}^{*}(\phi)$ for $\phi\in\mathbb{R}$,
Proposition 2.6 gives
$\displaystyle\lim_{\Lambda\to\infty}\|f_{\Lambda}-f^{*}\|_{\mathcal{C}^{m}[a,b]}$
for $m\in\mathbb{N}$ and $a<b$.

For function $f^{*}$, we have the following propositions.

######  Proposition 2.7.

Function $\displaystyle f^{*}(\phi)=\sum_{i=1}^{N}z_{i}c_{i}^{*}(\phi)$ is
strictly decreasing on $\mathbb{R}$, where functions $c_{i}^{*}$ are obtained
in Proposition 2.5.

###### Proof.

Differentiating (1.14) and (1.15) with respect to $\phi$ gives

$\displaystyle\sum_{j=0}^{N}\lambda_{j}\frac{\mathrm{d}c_{j}^{*}}{\mathrm{d}\phi}=0,$
(2.12)
$\displaystyle\frac{\mathrm{d}c_{i}^{*}}{\mathrm{d}\phi}=\frac{\lambda_{i}}{\lambda_{0}}\frac{c_{i}^{*}}{c_{0}^{*}}\frac{\mathrm{d}c_{0}^{*}}{\mathrm{d}\phi}-z_{i}c_{i}^{*},\quad\text{for}~{}\phi\in\mathbb{R}~{}\text{and}~{}i=1,\dots,N.$
(2.13)

Multiply the equation (2.13) by $\lambda_{i}$ for $i=1,\dots,N$ and add them
together. Then by equation (2.12), we obtain

$\frac{\mathrm{d}c_{0}^{*}}{\mathrm{d}\phi}=\lambda_{0}c_{0}^{*}\frac{\sum_{i=1}^{N}z_{i}\lambda_{i}c_{i}^{*}}{\sum_{i=0}^{N}\lambda_{i}^{2}c_{i}^{*}}\quad\text{for}~{}\phi\in\mathbb{R}.$

Consequently, we have

$\frac{\mathrm{d}f^{*}}{\mathrm{d}\phi}=\sum_{i=1}^{N}z_{i}\frac{\mathrm{d}c_{i}^{*}}{\mathrm{d}\phi}=\frac{1}{\lambda_{0}c_{0}^{*}}\frac{\mathrm{d}c_{0}^{*}}{\mathrm{d}\phi}\sum_{i=1}^{N}z_{i}\lambda_{i}c_{i}^{*}-\sum_{i=1}^{N}z_{i}^{2}c_{i}^{*}=\frac{(\sum_{i=1}^{N}z_{i}\lambda_{i}c_{i}^{*})^{2}}{\sum_{i=0}^{N}\lambda_{i}^{2}c_{i}^{*}}-\sum_{i=1}^{N}z_{i}^{2}c_{i}^{*}<0\quad\text{for}~{}\phi\in\mathbb{R}.$

Here the last inequality comes from Cauchy’s inequality, and the fact that
$\lambda_{i}>0$, $c_{i}^{*}>0$ for $i\in\\{0,1,\dots,N\\}$, and $z_{i}z_{j}<0$
for some $i,j\in\\{1,\dots,N\\}$. Therefore, we complete the proof of
Proposition 2.7. ∎

######  Proposition 2.8.

Function $f^{*}$ is bounded and $m^{*}<f^{*}(\phi)<M^{*}$ for all
$\phi\in\mathbb{R}$, where $\displaystyle
m^{*}=\lim_{\phi\to\infty}f^{*}(\phi)<0$ and $\displaystyle
M^{*}=\lim_{\phi\to-\infty}f^{*}(\phi)>0$.

###### Proof.

To prove boundedness of $f^{*}$, we use Proposition 2.5 and equation 1.15 to
get

$0<c_{i}^{*}(\phi)<\tilde{\mu}_{0}/\lambda_{i}\quad\text{for}~{}\phi\in\mathbb{R}~{}\text{and}~{}i=0,1,\dots,N,$
(2.14)

and

$|f^{*}(\phi)|\leq\sum_{i=1}^{N}|z_{i}|c_{i}^{*}(\phi)\leq\tilde{\mu}_{0}\sum_{i=1}^{N}\frac{|z_{i}|}{\lambda_{i}}\leq\tilde{\mu}_{0}N\max_{0\leq
i\leq N}\frac{|z_{i}|}{\lambda_{i}}\quad\text{for}~{}\phi\in\mathbb{R}.$
(2.15)

On the other hand, by (1.14) and (1.16), function $f^{*}$ can be represented
as

$f^{*}(\phi)=\sum_{i=1}^{N}z_{i}(c_{0}^{*}(\phi))^{\lambda_{i}/\lambda_{0}}\exp(\bar{\mu}_{i}-z_{i}\phi)\quad\text{for}~{}\phi\in\mathbb{R}.$
(2.16)

By (2.15) and Proposition 2.7, function $f^{*}$ is strictly decreasing and
bounded on $\mathbb{R}$ so the limit
$\displaystyle\lim_{\phi\to\infty}f^{*}(\phi)$, denoted by $m^{*}$, exists and
is finite. By (1.14) and (2.14), we obtain

$\lim_{\phi\to\infty}c_{0}^{*}(\phi)=\lim_{\phi\to\infty}\left[c_{i}^{*}(\phi)\exp(-\bar{\mu}_{i}+z_{i}\phi)\right]^{\lambda_{0}/\lambda_{i}}=0\quad\text{for}~{}z_{i}<0.$
(2.17)

Moreover, we use (1.14) and (2.17) to get

$\lim_{\phi\to\infty}c_{i}^{*}(\phi)=\lim_{\phi\to\infty}\left[(c_{0}^{*}(\phi))^{\lambda_{i}/\lambda_{0}}\exp(\bar{\mu}_{i}-z_{i}\phi)\right]=0\quad\text{for}~{}z_{i}>0,$
(2.18)

and hence (2.16) and (2.18) imply $\displaystyle
m^{*}=\lim_{\phi\to\infty}f^{*}(\phi)\leq 0$. Now, we prove $m^{*}<0$ by
contradiction. Suppose $m^{*}=0$. Then (2.18) gives
$\displaystyle\lim_{\phi\to\infty}c_{i}^{*}(\phi)=0$ for $i=0,1,\dots,N$ which
contradicts with equation (1.15) (by letting $\phi\to\infty$) and
$\tilde{\mu}_{0}>0$. Similarly, we may prove
$\displaystyle\lim_{\phi\to-\infty}f^{*}(\phi)=M^{*}>0$ and complete the proof
of Proposition 2.8. ∎

##  3 Proof of Theorem 1.1

The existence and uniqueness of $\phi_{\Lambda}$ can be proved by the standard
Direct method. One may refer the proof in Appendix I. The uniform boundedness
of $\phi_{\Lambda}$ and the convergence of $\phi_{\Lambda}$ are proved as
follows.

###  3.1 Uniform boundedness of $\phi_{\Lambda}$ and
$c_{0,\Lambda}(\phi_{\Lambda})$

######  Lemma 3.1.

There exists a constant $M_{0}\geq 1$ independent of $\Lambda$ such that
$\displaystyle\max_{x\in\overline{\Omega}}c_{0,\Lambda}(\phi_{\Lambda}(x))\leq
M_{0}$ for $\Lambda\geq 1$.

###### Proof.

Let $\displaystyle
M_{0}=\max\left\\{\frac{\tilde{\mu_{0}}}{\lambda_{0}}+\frac{\hat{\mu_{0}}}{\lambda_{0}^{2}},1\right\\}\geq
1$ and $\Lambda\geq 1$. We claim that $c_{0,\Lambda}(\phi_{\Lambda}(x))\leq
M_{0}$ for $x\in\overline{\Omega}$. Suppose that
$\Omega_{1}=\\{x\in\overline{\Omega}:c_{0,\Lambda}(\phi_{\Lambda}(x))>1\\}$ is
nonempty. Otherwise, due to $M_{0}\geq 1$,
$c_{0,\Lambda}(\phi_{\Lambda}(x))\leq 1\leq M_{0}$ for $x\in\overline{\Omega}$
which is trivial. By (2.4) and (A2), we obtain

$c_{0,\Lambda}(\phi_{\Lambda}(x))\leq\frac{\tilde{\mu}_{0}}{\lambda_{0}}+\frac{\hat{\mu}_{0}}{\Lambda\lambda_{0}^{2}}\leq
M_{0}$

for $x\in\Omega_{1}$. Here we have used the fact that $\ln
c_{0,\Lambda}(\phi_{\Lambda}(x))>0$ for $x\in\Omega_{1}$, $\Lambda\geq 1$, and
$z_{0}=\bar{\mu}_{0}=0$. Therefore, we complete the proof of Lemma 3.1. ∎

######  Lemma 3.2.

There exists a positive constant $M_{1}$ independent of $\Lambda$ such that
$\left\|\phi_{\Lambda}\right\|_{L^{\infty}(\Omega)}\leq M_{1}$ for
$\Lambda\geq 1$.

###### Proof.

Let $\psi$ be the solution of the equation
$-\nabla\cdot(\epsilon\nabla\psi)=\rho_{0}$ in $\Omega$ with the Robin
boundary condition $\psi+\eta\frac{\partial\psi}{\partial\nu}=0$ on
$\partial\Omega$, and let $\bar{\phi}_{\Lambda}=\phi_{\Lambda}-\psi$. Then
function $\bar{\phi}_{\Lambda}$ satisfies

$\nabla\cdot(\epsilon\nabla\bar{\phi}_{\Lambda})=f_{\Lambda}(\phi_{\Lambda})\quad\text{in}~{}\Omega.$
(3.1)

By (2.3), $c_{i,\Lambda}(\phi_{\Lambda})$ satisfies

$c_{i,\Lambda}(\phi_{\Lambda})=\left(c_{0,\Lambda}(\phi_{\Lambda})\right)^{\lambda_{i}/\lambda_{0}}\exp\left(\bar{\mu}_{i}-z_{i}\psi-
z_{i}\bar{\phi}_{\Lambda}\right)\quad\text{for}~{}i=1,\dots,N.$ (3.2)

Since $\psi$ is independent of $\Lambda$ and is continuous on
$\overline{\Omega}$, then $\bar{\phi}_{\Lambda}$ is uniformly bounded if and
only if$\phi_{\Lambda}$ is uniformly bounded. Thus, it suffices to show that
$\displaystyle\max_{x\in\overline{\Omega}}\bar{\phi}_{\Lambda}(x)\leq M_{1}$
and $\displaystyle\min_{x\in\overline{\Omega}}\bar{\phi}_{\Lambda}\geq-M_{1}$
for $\Lambda\geq 1$, where $M_{1}$ is a positive constant independent of
$\Lambda$.

Now we prove that
$\displaystyle\max_{x\in\overline{\Omega}}\bar{\phi}_{\Lambda}(x)\leq M_{1}$
for $\Lambda\geq 1$, where $M_{1}$ is a positive constant independent of
$\Lambda$. Suppose by contradiction that there exists a sequence $\Lambda_{k}$
with $\displaystyle\lim_{k\to\infty}\Lambda_{k}=\infty$ such that
$\displaystyle\max_{\overline{\Omega}}\bar{\phi}_{\Lambda_{k}}\geq k$ for
$k\in\mathbb{N}$. Then there exists $x_{k}\in\Omega$ such that
$\displaystyle\bar{\phi}_{\Lambda_{k}}(x_{k})=\max_{\overline{\Omega}}\bar{\phi}_{\Lambda_{k}}$
which implies $\nabla\bar{\phi}_{\Lambda_{k}}\left(x_{k}\right)=0$ and
$\Delta\bar{\phi}_{\Lambda_{k}}(x_{k})\leq 0$. Note that because of the Robin
boundary condition of $\bar{\phi}_{\Lambda_{k}}$, maximum point $x_{k}$ cannot
be located on the boundary $\partial\Omega$ as $k$ sufficiently large. Hence
without loss of generality, we assume each $x_{k}\in\Omega$ for
$k\in\mathbb{N}$. For a sake of simplicity, in this proof, we set
$c_{i,k}:=c_{i,\Lambda_{k}}$, $\phi_{k}:=\phi_{\Lambda_{k}}$,
$f_{k}:=f_{\Lambda_{k}}$, and $\bar{\phi}_{k}:=\bar{\phi}_{\Lambda_{k}}$.
Hence by equation (3.1) with $\nabla\bar{\phi}_{k}(x_{k})=0$,
$\Delta\bar{\phi}_{k}(x_{k})\leq 0$ and function $\epsilon$ is positive, we
have

$0\leq-\nabla\epsilon\left(x_{k}\right)\cdot\nabla\bar{\phi}_{k}\left(x_{k}\right)-\epsilon\left(x_{k}\right)\Delta\bar{\phi}_{k}\left(x_{k}\right)=f_{k}(\phi_{k}(x_{k})).$
(3.3)

Because $\displaystyle f_{k}(\phi_{k})=\sum_{i=1}^{N}z_{i}c_{i,k}(\phi_{k})$,
we may use (3.2) and (3.3) to get

$\displaystyle 0$
$\displaystyle<\sum_{z_{i}<0}\left(-z_{i}\right)[c_{0,k}(\phi_{k}(x_{k}))]^{\lambda_{i}/\lambda_{0}}\exp\left(\bar{\mu}_{i}-z_{i}\psi\left(x_{k}\right)-z_{i}\bar{\phi}_{k}\left(x_{k}\right)\right)$
(3.4)
$\displaystyle\leq\sum_{z_{i}>0}z_{i}[c_{0,k}(\phi_{k}(x_{k}))]^{\lambda_{i}/\lambda_{0}}\exp\left(\bar{\mu}_{i}-z_{i}\psi\left(x_{k}\right)-z_{i}\bar{\phi}_{k}\left(x_{k}\right)\right).$

Applying Lemma 3.1 to (3.2) for $z_{i}>0$, we obtain
$\displaystyle\lim_{k\to\infty}c_{i,k}(\phi_{k}(x_{k}))=0$ for all $i$ with
$z_{i}>0$. Then (3.2) and (3.4) give

$\lim_{k\to\infty}\sum_{z_{i}<0}\left(-z_{i}\right)[c_{0,k}(\phi_{k}(x_{k}))]^{\lambda_{i}/\lambda_{0}}\exp\left(\bar{\mu}_{i}-z_{i}\psi\left(x_{k}\right)-z_{i}\bar{\phi}_{k}\left(x_{k}\right)\right)=0,$

which means $\displaystyle\lim_{k\to\infty}c_{i,k}(\phi_{k}(x_{k}))=0$ for all
$i$ with $z_{i}<0$. Inserting $x=x_{k}$ into (3.2) with $z_{i}<0$, we can
obtain

$c_{0,k}(\phi_{k}(x_{k}))=[(c_{i,k}(\phi_{k}(x_{k}))]^{\lambda_{0}/\lambda_{i}}\left[\exp\left(-\bar{\mu}_{i}+z_{i}\psi\left(x_{k}\right)+z_{i}\bar{\phi}_{k}\left(x_{k}\right)\right)\right]^{\lambda_{0}/\lambda_{i}}\to
0$

as $k\to\infty$. Hence
$\displaystyle\lim_{k\to\infty}c_{i,k}(\phi_{k}(x_{k}))=0$ for
$i=0,1,\dots,N$. By (2.3), (2.4) and (A2), we get the following contradiction:

$0\geq\lim_{k\to\infty}\frac{\ln
c_{0,k}(\phi_{k}(x_{k}))}{\Lambda_{k}}=\lim_{k\to\infty}\left(\lambda_{0}\tilde{\mu}_{0}+\frac{\hat{\mu}_{0}}{\Lambda_{k}}-\sum_{i=0}^{N}\lambda_{i}c_{i,k}(\phi_{k}(x_{k}))\right)=\lambda_{0}\tilde{\mu}_{0}>0.$
(3.5)

Therefore, we complete the proof to show that
$\displaystyle\max_{x\in\overline{\Omega}}\bar{\phi}_{\Lambda}(x)\leq M_{1}$
for $\Lambda\geq 1$, where $M_{1}$ is a positive constant independent of
$\Lambda$.

It remains to prove that
$\displaystyle\min_{x\in\overline{\Omega}}\bar{\phi}_{\Lambda}\geq-M_{1}$ for
$\Lambda\geq 1$, where $M_{1}$ is a positive independent of $\Lambda$. Suppose
by contradiction that there exists a sequence $\Lambda_{k}$ with
$\displaystyle\lim_{k\to\infty}\Lambda_{k}=\infty$ such that
$\displaystyle\min_{\overline{\Omega}}\bar{\phi}_{k}\leq-k$ for
$k\in\mathbb{N}$. Then there exists $x_{k}\in\Omega$ such that
$\displaystyle\bar{\phi}_{k}(x_{k})=\min_{\overline{\Omega}}\bar{\phi}_{k}$,
which implies $\nabla\bar{\phi}_{k}(x_{k})=0$ and
$\Delta\bar{\phi}_{k}(x_{k})\geq 0$. Notice that because of the Robin boundary
condition of $\bar{\phi}_{\Lambda_{k}}$, minimum point $x_{k}$ cannot be
located on the boundary $\partial\Omega$ as $k$ sufficiently large. Hence
without loss of generality, we assume each $x_{k}\in\Omega$ for
$k\in\mathbb{N}$. Thus, as for (3.4), we have

$\displaystyle 0$
$\displaystyle<\sum_{z_{i}>0}z_{i}[c_{0,k}(\phi_{k}(x_{k}))]^{\lambda_{i}/\lambda_{0}}\exp\left(\bar{\mu}_{i}-z_{i}\psi\left(x_{k}\right)-z_{i}\bar{\phi}_{k}\left(x_{k}\right)\right)$
(3.6)
$\displaystyle\leq\sum_{z_{i}<0}(-z_{i})[c_{0,k}(\phi_{k}(x_{k}))]^{\lambda_{i}/\lambda_{0}}\exp\left(\bar{\mu}_{i}-z_{i}\psi\left(x_{k}\right)-z_{i}\bar{\phi}_{k}\left(x_{k}\right)\right).$

Applying Lemma 3.1 to (3.2) for $z_{i}<0$, we obtain
$\displaystyle\lim_{k\to\infty}c_{i,k}(\phi_{k}(x_{k}))=0$ for all $i$ with
$z_{i}<0$. Then (3.2) and (3.6) give
$\displaystyle\lim_{k\to\infty}c_{i,k}(\phi_{k}(x_{k}))=0$ for all $i$ with
$z_{i}>0$. Inserting $x=x_{k}$ into (3.2) with $z_{i}>0$, we have
$\displaystyle\lim_{k\to\infty}c_{0,k}(\phi_{k}(x_{k}))=0$, and hence
$\displaystyle\lim_{k\to\infty}c_{i,k}(\phi_{k}(x_{k}))=0$ for each
$i=0,1,\dots,N$. As for (3.5), we also get a contradiction and complete the
proof to show
$\displaystyle\min_{x\in\overline{\Omega}}\bar{\phi}_{\Lambda}(x)\geq-M_{1}$
for $\Lambda\geq 1$, where $M_{1}$ is a positive constant independent of
$\Lambda$. Therefore, we complete the proof of Lemma 3.2.∎

By Lemma 3.2, we have

######  Lemma 3.3.

There exists positive constant $M_{2}$ independent of $\Lambda$ such that
$\displaystyle\min_{x\in\overline{\Omega}}c_{0,\Lambda}(\phi_{\Lambda}(x))\geq
M_{2}$ for $\Lambda\geq 1$.

###### Proof.

Suppose by contrary that there exists $\Lambda_{k}$ with
$\displaystyle\lim_{k\to\infty}\Lambda_{k}=\infty$, and
$x_{k}\in\overline{\Omega}$ is the minimum point of
$c_{0,\Lambda_{k}}(\phi_{\Lambda_{k}}(x))$ such that
$\displaystyle\lim_{k\to\infty}c_{0,\Lambda_{k}}(\phi_{\Lambda_{k}}(x_{k}))=0$.
Notice that by (2.4), $c_{0,\Lambda}(\phi_{\Lambda}(x))>0$ for $\Lambda\geq 1$
and $x\in\overline{\Omega}$. Then, there exists $N_{1}\in\mathbb{N}$ such that
$\ln c_{0,\Lambda_{k}}(\phi_{\Lambda_{k}}(x_{k}))<0$ for $k>N_{1}$. By (2.3),
(2.4) at $\phi=\phi_{\Lambda_{k}}(x_{k})$ and assumption (A2), we get

$\sum_{i=0}^{N}\lambda_{i}c_{i,\Lambda_{k}}(\phi_{\Lambda_{k}}(x_{k}))>\tilde{\mu}_{0}+\frac{\hat{\mu}_{0}}{\Lambda_{k}\lambda_{0}}\quad\text{for}~{}k>N_{1}.$
(3.7)

On the other hand, by (3.2) and Lemma 3.2, we obtain
$\displaystyle\lim_{k\to\infty}c_{i,\Lambda_{k}}(\phi_{\Lambda_{k}}(x_{k}))=0$
for $i=0,1,\dots,N$, which contradicts with (3.7) as $k\to\infty$. Hence we
complete the proof of Lemma 3.3. ∎

######  Remark 3.

By Lemmas 3.1–3.3 and (2.3), we have $M_{3}\leq
c_{i,\Lambda}(\phi_{\Lambda}(x))\leq M_{4}$ for $x\in\overline{\Omega}$,
$\Lambda\geq 1$ and $i=0,1,\dots,N$, where $M_{3}$ and $M_{4}$ are positive
constants independent of $\Lambda$. Moreover, due to $\displaystyle
f_{\Lambda}(\phi_{\Lambda})=\sum_{i=1}^{N}z_{i}c_{i,\Lambda}(\phi_{\Lambda})$,
$\|f_{\Lambda}(\phi_{\Lambda})\|_{L^{\infty}}\leq M_{5}$ for $\Lambda\geq 1$,
where $M_{5}$ is a positive constant independent of $\Lambda$.

###  3.2 Convergence of $\phi_{\Lambda}$

Since $\phi_{\Lambda}$ is the solution of equation
$-\nabla\cdot(\epsilon\nabla\phi_{\Lambda})=\rho_{0}+f_{\Lambda}(\phi_{\Lambda})$
in $\Omega$ with the Robin boundary condition
$\phi_{\Lambda}+\eta\frac{\partial\phi_{\Lambda}}{\partial\nu}=\phi_{bd}$ on
$\partial\Omega$, we use the $W^{2,p}$ estimate (cf. [2, Theorem 15.2]) to get

$\|\phi_{\Lambda}\|_{W^{2,p}(\Omega)}\leq
C(\|\rho_{0}+f_{\Lambda}(\phi_{\Lambda})\|_{L^{p}(\Omega)}+\|\phi_{bd}\|_{W^{1,p}(\Omega)})$

for all $p>1$, where $C$ is a positive constant independent of $\Lambda$.
Hence by Remark 3, we have the uniform bound estimate of $\phi_{\Lambda}$ in
$W^{2,p}$ norm. This implies that there exists a sequence of functions
$\\{\phi_{\Lambda_{k}}\\}_{k=1}^{\infty}$, with
$\displaystyle\lim_{k\to\infty}\Lambda_{k}=\infty$, such that
$\phi_{\Lambda_{k}}$ converges to $\phi^{*}$ weakly in $W^{2,p}(\Omega)$. By
Remark 2 of Proposition 2.6, we have the convergence of function
$f_{\Lambda_{k}}$ to $f^{*}$ in $\mathcal{C}^{m}[-M_{1},M_{1}]$ for
$m\in\mathbb{N}$ so $\phi^{*}$ satisfies equation (1.17) in weak sense, where
the positive constant $M_{1}$ comes from Lemma 3.2. Let
$w_{k}=\phi_{\Lambda_{k}}-\phi^{*}$, $c_{i,k}:=c_{i,\Lambda_{k}}$, and
$f_{k}:=f_{\Lambda_{k}}$. Then by Sobolev embedding,
$w_{k}\in\mathcal{C}^{1,\alpha}(\Omega)$ for $\alpha\in(0,1)$, and
$\displaystyle\lim_{k\to\infty}\|w_{k}\|_{\mathcal{C}^{1,\alpha}(\Omega)}=0$.
Moreover, $w_{k}$ satisfies

$-\nabla\cdot(\epsilon\nabla
w_{k})=f_{k}(w_{k}+\phi^{*})-f^{*}(\phi^{*})\quad\text{in}~{}\Omega$

with the boundary condition $w_{k}+\eta\frac{\partial w_{k}}{\partial\nu}=0$
on $\partial\Omega$. Using the Schauder’s estimate (cf. [22, Theorem 6.30])
with the mathematical induction, we get

$\displaystyle\|w_{k}\|_{\mathcal{C}^{m+2,\alpha}(\Omega)}$ $\displaystyle\leq
C\|f_{k}(w_{k}+\phi^{*})-f^{*}(\phi^{*})\|_{\mathcal{C}^{m,\alpha}(\Omega)}$
(3.8) $\displaystyle\leq
C^{\prime}\sum_{i=1}^{N}\left(\|c_{i,k}(w_{k}+\phi^{*})-c_{i}^{*}(w_{k}+\phi^{*})\|_{\mathcal{C}^{m,\alpha}(\Omega)}+\|c_{i}^{*}(w_{k}+\phi^{*})-c_{i}^{*}(\phi^{*})\|_{\mathcal{C}^{m,\alpha}(\Omega)}\right),$

for all $m\in\mathbb{N}$ and $\alpha\in(0,1)$, where $C$ and $C^{\prime}$ are
positive constants independent of $k$. By Proposition 2.6 and induction
hypothesis
$\displaystyle\lim_{k\to\infty}\|w_{k}\|_{\mathcal{C}^{m,\alpha}(\Omega)}=0$,
we may use (3.8) to get
$\displaystyle\lim_{k\to\infty}\|w_{k}\|_{\mathcal{C}^{m+2,\alpha}(\Omega)}=0$,
i.e.
$\displaystyle\lim_{k\to\infty}\|\phi_{\Lambda_{k}}-\phi^{*}\|_{C^{m+2,\alpha}(\Omega)}=0$
for $m\in\mathbb{N}$ and $\alpha\in(0,1)$. Therefore, $\phi^{*}$ is the
solution of equation (1.17) with the Robin boundary condition (1.8).

To complete the proof of Theorem 1.1, we need to prove

######  Claim 1.

For any $m\in\mathbb{N}$, we have
$\displaystyle\lim_{\Lambda\to\infty}\|\phi_{\Lambda}-\phi^{*}\|_{\mathcal{C}^{m}(\Omega)}=0$.

###### Proof.

Suppose that there exist the sequences $\\{\Lambda_{k}\\}$ and
$\\{\tilde{\Lambda}_{k}\\}$ tending to infinity such that sequences
$\\{\phi_{\Lambda_{k}}\\}$ and $\\{\phi_{\tilde{\Lambda}_{k}}\\}$ have limits
$\phi_{1}^{*}$ and $\phi_{2}^{*}$, respectively. It is clear that
$\phi_{1}^{*}$ and $\phi_{2}^{*}$ satisfy the equation (1.17) with the Robin
boundary condition (1.8). Now we want to prove that
$\phi_{1}^{*}\equiv\phi_{2}^{*}$.

Let $u=\phi_{1}^{*}-\phi_{2}^{*}$. Subtracting (1.17) with
$\phi^{*}=\phi_{2}^{*}$ from that with $\phi^{*}=\phi_{1}^{*}$, we obtain

$-\nabla\cdot(\epsilon\nabla
u)=f^{*}(\phi_{1}^{*})-f^{*}(\phi_{2}^{*})=c(x)u\quad\text{in}~{}\Omega,$

where function $c$ is defined by

$c(x)=\begin{cases}\displaystyle\frac{f^{*}(\phi_{1}^{*}(x))-f^{*}(\phi_{2}^{*}(x))}{\phi_{1}^{*}(x)-\phi_{2}^{*}(x)},&\text{if}~{}\phi_{1}^{*}(x)\neq\phi_{2}^{*}(x);\\\
\displaystyle\frac{\mathrm{d}f^{*}}{\mathrm{d}\phi}(\phi_{1}^{*}(x)),&\text{if}~{}\phi_{1}^{*}(x)=\phi_{2}^{*}(x).\end{cases}$

By Proposition 2.7, we have $c<0$ in $\Omega$ which comes from the fact that
if $f$ is a strictly decreasing function on $\mathbb{R}$, then
$\frac{f(\alpha)-f(\beta)}{\alpha-\beta}<0$ for $\alpha\neq\beta$. Since
$\nabla\cdot(\epsilon\nabla u)+c(x)u=0$ with $c<0$ in $\Omega$, it is obvious
that $u$ cannot be a nonzero constant. Then by the strong maximum principle,
we have that $u$ attains its nonnegative maximum value and nonnpositive
minimum value at the boundary point. Suppose $u$ has nonnegative maximum value
and attains its maximum value at $x^{*}\in\partial\Omega$. Then by the
boundary condition of $u$ which is $u+\eta\frac{\partial u}{\partial\nu}=0$ on
$\partial\Omega$, we get $u(x^{*})=-\eta\frac{\partial
u}{\partial\nu}(x^{*})\leq 0$, and hence $u\leq u(x^{*})\leq 0$ on
$\overline{\Omega}$. Similarly, we obtain $u\geq 0$ in $\overline{\Omega}$,
and hence $u\equiv 0$. Therefore, we conclude that
$\phi_{1}^{*}\equiv\phi_{2}^{*}$ and complete the proof of Theorem 1.1. ∎

##  4 Proof of Theorem 1.2

We show the asymptotic limit of $\phi_{\Lambda}$ for the following two cases:
(i) function $\rho_{0}$ is constant; (ii) function $\rho_{0}$ is nonconstant
and satisfies $\int_{\Omega}\rho_{0}\,\mathrm{d}x=0$ in Sections 4.1 and 4.2,
respectively. For the case (i), the unique constant solution $\phi_{\Lambda}$
follows from Proposition 2.2 and 2.4, and the uniform convergence of
$\phi_{\Lambda}$ can be obtained from Remark 2 of Proposition 2.6. For the
case (ii), the existence, uniqueness, and regularity of $\phi_{\Lambda}$ can
be proved by the standard Direct method, which is stated in Appendix I. Unlike
Lemma 3.2, we use Hopf’s lemma to prove the uniform boundedness of
$\phi_{\Lambda}$ with the Neumann boundary condition (See Lemma 4.1). By
Lemmas 3.1, 3.3 and 4.1, we can follow the same argument in Section 3.2 to
complete the proof of Theorem 1.2(ii).

###  4.1 Constant $\rho_{0}$

The following is the proof of Theorem 1.2(i).

###### Proof of Theorem 1.2(i).

Firstly, we fix $\Lambda>0$ and suppose that $\rho_{0}$ is constant. Then by
the strict monotonicity and unboundedness of $f_{\Lambda}$ (see Propositions
2.2 and 2.4), there exists a unique $\hat{\phi}_{\Lambda}\in\mathbb{R}$ such
that $\rho_{0}+f_{\Lambda}(\hat{\phi}_{\Lambda})=0$. Hence
$\phi_{\Lambda}\equiv\hat{\phi}_{\Lambda}$ is a constant solution of equation
(1.12) with the Neumann boundary condition (1.9). For the uniqueness of
solution $\phi_{\Lambda}$ up to a constant, suppose that $\phi_{1,\Lambda}$
and $\phi_{2,\Lambda}$ are the solutions of equation (1.12) with the Neumann
boundary condition (1.9). We subtract the equation (1.7) for
$\phi=\phi_{2,\Lambda}$ from that for $\phi=\phi_{1,\Lambda}$ to get

$-\nabla\cdot(\epsilon\nabla(\phi_{1,\Lambda}-\phi_{2,\Lambda}))=\sum_{i=1}^{N}z_{i}\left[c_{i,\Lambda}(\phi_{1,\Lambda})-c_{i,\Lambda}(\phi_{2,\Lambda})\right]=f_{\Lambda}(\phi_{1,\Lambda})-f_{\Lambda}(\phi_{2,\Lambda})\quad\text{in}~{}\Omega.$
(4.1)

Multiplying the equation (4.1) by $\phi_{1,\Lambda}-\phi_{2,\Lambda}$ and
integrating it over $\Omega$, we have

$-\int_{\Omega}(\phi_{1,\Lambda}-\phi_{2,\Lambda})\nabla\cdot(\epsilon\nabla(\phi_{1,\Lambda}-\phi_{2,\Lambda}))\mathrm{d}x=\int_{\Omega}(\phi_{1,\Lambda}-\phi_{2,\Lambda})[f_{\Lambda}(\phi_{1,\Lambda})-f_{\Lambda}(\phi_{2,\Lambda})]\mathrm{d}x.$

Then using integration by parts and the Neumann boundary condition, we obtain

$\int_{\Omega}\epsilon|\nabla(\phi_{1,\Lambda}-\phi_{2,\Lambda})|^{2}\mathrm{d}x=\int_{\Omega}(\phi_{1,\Lambda}-\phi_{2,\Lambda})[f_{\Lambda}(\phi_{1,\Lambda})-f_{\Lambda}(\phi_{2,\Lambda})]\mathrm{d}x\leq
0.$

Here we have used the monotone decreasing of $f_{\Lambda}$ (see Proposition
2.2) and function $\epsilon$ is positive. Thus,
$\phi_{1,\Lambda}-\phi_{2,\Lambda}\equiv C$ for some constant $C$, and we get
the uniqueness of $\phi_{\Lambda}$ up to a constant. This implies that
$\phi_{\Lambda}$ is a constant solution and satisfies
$\phi_{\Lambda}=\hat{\phi}_{\Lambda}+C$. Due to the uniqueness of equation
$\rho_{0}+f_{\Lambda}(\hat{\phi}_{\Lambda})=0$, we obtain the constant $C=0$
and $\phi_{\Lambda}\equiv\hat{\phi}_{\Lambda}$. By Propositions 2.7 and 2.8,
$f^{*}$ is strictly decreasing and the range of $f^{*}$ is exactly the
interval $(m^{*},M^{*})$, there exists a unique $\phi^{*}\in\mathbb{R}$ such
that $\rho_{0}+f^{*}(\phi^{*})=0$ for $\rho_{0}\in(-M^{*},-m^{*})$. By Remark
2 of Proposition 2.6, we conclude that $\phi_{\Lambda}$ is uniformly bounded
and $\displaystyle\lim_{\Lambda\to\infty}\hat{\phi}_{\Lambda}=\phi^{*}$. The
proof of Theorem 1.2(i) is complete. ∎

###  4.2 Uniform boundedness of $\phi_{\Lambda}$ and $c_{0,\Lambda}$

Under the Neumann boundary condition, Lemma 3.1 also holds true so
$\displaystyle\max_{x\in\overline{\Omega}}c_{0,\Lambda}(\phi_{\Lambda}(x))\leq
M_{0}$ for $\Lambda\geq 1$, where $M_{0}$ is a positive constant independent
of $\Lambda$. However, due to the Neumann boundary condition, we have to use
the Hopf’s Lemma to modify the proof of Lemma 3.2 for the uniform boundedness
of $\phi_{\Lambda}$ (see Lemma 4.1).Then as for Lemma 3.3, we can use Lemma
4.1 to prove
$\displaystyle\min_{x\in\overline{\Omega}}c_{0,\Lambda}(\phi_{\Lambda}(x))\geq
M_{2}$ for $\Lambda\geq 1$, where $M_{2}$ is a positive constant independent
of $\Lambda$. Now we state and prove Lemma 4.1 as follows.

######  Lemma 4.1.

Under the same hypotheses of Theorem 1.2(ii), let $\phi_{\Lambda}$ be the
solution of equation (1.12) with the Neumann boundary condition (1.9). Then
there exists a positive constant $M_{6}$ independent of $\Lambda$ such that
$\|\phi_{\Lambda}\|_{L^{\infty}(\Omega)}\leq M_{6}$ for $\Lambda\geq 1$.

###### Proof.

Based on $\int_{\Omega}\rho_{0}\,\mathrm{d}x=0$, let $\psi_{1}$ and $\psi_{2}$
be the solution of Poisson equation
$-\nabla\cdot(\epsilon\nabla\psi)=\rho_{0}$ in $\Omega$ with the Neumann
boundary condition (1.9) such that
$\displaystyle\min_{\overline{\Omega}}\psi_{1}=0$ and
$\displaystyle\max_{\overline{\Omega}}\psi_{2}=0$. Note that $\psi_{1}$ and
$\psi_{2}$ are independent of $\Lambda$, and up to a constant. Let
$\bar{\phi}_{i,\Lambda}=\phi_{\Lambda}-\psi_{i}$ for $i=1,2$. Then
$\bar{\phi}_{i,\Lambda}$ satisfies

$-\nabla\cdot(\epsilon\nabla\bar{\phi}_{i,\Lambda})=f_{\Lambda}(\phi_{\Lambda})\quad\text{in}~{}\Omega,\quad\text{for}~{}i=1,2.$
(4.2)

Since $\psi_{i}$ is continuous on $\overline{\Omega}$, then
$\bar{\phi}_{i,\Lambda}$ is uniformly bounded if and only if $\phi_{\Lambda}$
is uniformly bounded. Thus, it suffices to show that
$\displaystyle\max_{x\in\overline{\Omega}}\bar{\phi}_{1,\Lambda}(x)\leq M_{7}$
and $\displaystyle\min_{x\in\overline{\Omega}}\bar{\phi}_{2,\Lambda}(x)\geq-
M_{7}$ for $\Lambda\geq 1$, where $M_{7}$ is a positive constant independent
of $\Lambda$.

Firstly, we prove that there exists a positive constant $M_{7}$ independent of
$\Lambda$ such that
$\displaystyle\max_{x\in\overline{\Omega}}\bar{\phi}_{1,\Lambda}(x)\leq M_{7}$
for $\Lambda\geq 1$. Suppose by contradiction that there exists a sequence
$\Lambda_{k}$ with $\displaystyle\lim_{k\to\infty}\Lambda_{k}=\infty$ such
that $\displaystyle\max_{\overline{\Omega}}\bar{\phi}_{1,\Lambda_{k}}\geq k$
for $k\in\mathbb{N}$. As for the proof of Lemma 3.2, we need to prove

######  Claim 2.

$\bar{\phi}_{1,\Lambda_{k}}$ attains its maximum value at $x_{k}\in\Omega$ for
sufficient large $k$.

###### Proof.

Suppose by contradiction that the maximum point $x_{k}$ of
$\bar{\phi}_{1,\Lambda_{k}}$ is located on the boundary $\partial\Omega$,
i.e., $x_{k}\in\partial\Omega$ for $k\in\mathbb{N}$. By Theorem 1.2(i), there
exists a unique constant $\hat{\phi}_{\Lambda}\in\mathbb{R}$ such that
$f_{\Lambda}(\hat{\phi}_{\Lambda})=0$ and $|\hat{\phi}_{\Lambda}|\leq M$ for
$\Lambda\geq 1$, where $M$ is a positive constant independent of $\Lambda$.
For notation convenience, let
$\tilde{\phi}_{k}:=\bar{\phi}_{1,\Lambda_{k}}-\hat{\phi}_{\Lambda_{k}}=\phi_{\Lambda_{k}}-\psi_{1}-\hat{\phi}_{\Lambda_{k}}$,
$\phi_{k}:=\phi_{\Lambda_{k}}$, $\bar{\phi}_{k}:=\bar{\phi}_{1,\Lambda_{k}}$,
$\hat{\phi}_{k}:=\hat{\phi}_{\Lambda_{k}}$, and $f_{k}:=f_{\Lambda_{k}}$.
Since $\hat{\phi}_{k}$ is a constant with $f_{k}(\hat{\phi}_{k})=0$ and
$|\hat{\phi}_{k}|\leq M$ for $k\in\mathbb{N}$, then equation (4.2) can be
transformed to

$-\nabla\cdot(\epsilon\nabla\tilde{\phi}_{k})=f_{k}(\tilde{\phi}_{k}+\hat{\phi}_{k}+\psi_{1})-f_{k}(\hat{\phi}_{k})=c(x)(\tilde{\phi}_{k}+\psi_{1})\quad\text{in}~{}\Omega,$
(4.3)

where function $c=c(x)$ is defined by

$c(x)=\begin{cases}\displaystyle\frac{f_{k}(\tilde{\phi}_{k}(x)+\hat{\phi}_{k}+\psi_{1}(x))-f_{k}(\hat{\phi}_{k})}{(\tilde{\phi}_{k}(x)+\hat{\phi}_{k}+\psi_{1}(x))-\hat{\phi}_{k}},&\text{if~{}}\tilde{\phi}_{k}(x)+\psi_{1}(x)\neq
0;\\\
\displaystyle\frac{\mathrm{d}f_{k}}{\mathrm{d}\phi}(\hat{\phi}_{k}),&\text{if}~{}\tilde{\phi}_{k}(x)+\psi_{1}(x)=0.\end{cases}$

By equation (4.3), we have
$\nabla\cdot(\epsilon\nabla\tilde{\phi}_{k})+c(x)\tilde{\phi}_{k}=-c(x)\psi_{1}\geq
0$ in $\Omega$. Here we have used that
$\displaystyle\min_{\overline{\Omega}}\psi_{1}=0$ and $c<0$ in $\Omega$
because of Proposition 2.2. Notice that if $f$ is a strictly decreasing
function on $\mathbb{R}$, then $\frac{f(\alpha)-f(\beta)}{\alpha-\beta}<0$ for
$\alpha\neq\beta$. Since $\hat{\phi}_{k}$ is a constant, $\tilde{\phi}_{k}$
and $\bar{\phi}_{k}$ have the same maximum point $x_{k}$ for $k\in\mathbb{N}$.
Moreover,
$\tilde{\phi}_{k}(x_{k})=\bar{\phi}_{1,\Lambda_{k}}(x_{k})-\hat{\phi}_{\Lambda_{k}}\geq
k-M>0$ if $k>M$. Hence by Hopf’s lemma, we get
$\frac{\partial\tilde{\phi}_{k}}{\partial\nu}(x_{k})>0$, which contradicts
with the Neumann boundary condition. Thus, $\bar{\phi}_{k}$ attains its
maximum value at interior point $x_{k}\in\Omega$ for $k>M$, and we complete
the proof of Claim 2. ∎

Claim 2 implies $\nabla\bar{\phi}_{1,\Lambda_{k}}(x_{k})=0$ and
$\Delta\bar{\phi}_{1,\Lambda_{k}}(x_{k})\leq 0$ for $k>M$. Then as for the
proof of (3.3)–(3.5), we conclude that
$\displaystyle\max_{x\in\overline{\Omega}}\bar{\phi}_{1,\Lambda}(x)\leq M_{7}$
for $\Lambda\geq 1$ which gives
$\displaystyle\max_{x\in\overline{\Omega}}\phi_{\Lambda}(x)\leq M_{6}$ for
$\Lambda\geq 1$, where $M_{6}$ and $M_{7}$ are positive constants independent
of $\Lambda$.

It remains to prove that
$\displaystyle\min_{x\in\overline{\Omega}}\bar{\phi}_{2,\Lambda}(x)\geq-M_{7}$
for $\Lambda\geq 1$, where $M_{7}$ is a positive constant independent of
$\Lambda$. Suppose by contradiction that there exists a sequence $\Lambda_{k}$
with $\displaystyle\lim_{k\to\infty}\Lambda_{k}=\infty$ such that
$\displaystyle\min_{\overline{\Omega}}\bar{\phi}_{2,\Lambda_{k}}\leq-k$ for
$k\in\mathbb{N}$. As for the proof of Lemma 3.2, we need to prove

######  Claim 3.

$\bar{\phi}_{2,\Lambda_{k}}$ attains its minimum value at $x_{k}\in\Omega$ for
sufficient large $k$.

###### Proof.

Suppose by contradiction that the minimum point $x_{k}$ of
$\bar{\phi}_{2,\Lambda_{k}}$ is located on the boundary $\partial\Omega$,
i.e., $x_{k}\in\partial\Omega$ for $k\in\mathbb{N}$. For a sake simplicity, we
use the similar notations of Claim 2 and let
$\tilde{\phi}_{k}:=\bar{\phi}_{2,\Lambda_{k}}-\hat{\phi}_{\Lambda_{k}}=\phi_{\Lambda_{k}}-\psi_{2}-\hat{\phi}_{\Lambda_{k}}$,
$\phi_{k}:=\phi_{\Lambda_{k}}$, $\bar{\phi}_{k}:=\bar{\phi}_{2,\Lambda_{k}}$,
$\hat{\phi}_{k}:=\hat{\phi}_{\Lambda_{k}}$. Then equation (4.2) becomes

$-\nabla\cdot(\epsilon\nabla\tilde{\phi}_{k})=f_{k}(\tilde{\phi}_{k}+\hat{\phi}_{k}+\psi_{2})-f_{k}(\hat{\phi}_{k})=c(x)(\tilde{\phi}_{k}+\psi_{2})\quad\text{in}~{}\Omega,$
(4.4)

where function $c=c(x)$ is defined by

$c(x)=\begin{cases}\displaystyle\frac{f_{k}(\tilde{\phi}_{k}(x)+\hat{\phi}_{k}+\psi_{2}(x))-f_{k}(\hat{\phi}_{k})}{(\tilde{\phi}_{k}(x)+\hat{\phi}_{k}+\psi_{2}(x))-\hat{\phi}_{k}},&\text{if~{}}\tilde{\phi}_{k}(x)+\psi_{2}(x)\neq
0;\\\
\displaystyle\frac{\mathrm{d}f_{k}}{\mathrm{d}\phi}(\hat{\phi}_{k}),&\text{if}~{}\tilde{\phi}_{k}(x)+\psi_{2}(x)=0.\end{cases}$

From equation (4.4), we have
$\nabla\cdot(\epsilon\nabla\tilde{\phi}_{k})+c(x)\tilde{\phi}_{k}=-c(x)\psi_{2}\leq
0$ in $\Omega$. Here we have used that
$\displaystyle\max_{\overline{\Omega}}\psi_{2}=0$ and $c<0$ in $\Omega$
because of Proposition 2.2. Notice that if $f$ is a strictly decreasing
function on $\mathbb{R}$, then $\frac{f(\alpha)-f(\beta)}{\alpha-\beta}<0$ for
$\alpha\neq\beta$. Since $\hat{\phi}_{k}$ is a constant, $\tilde{\phi}_{k}$
and $\bar{\phi}_{k}$ have the same minimum point $x_{k}$ for $k\in\mathbb{N}$.
Moreover,
$\tilde{\phi}_{k}(x_{k})=\bar{\phi}_{2,\Lambda_{k}}(x_{k})-\hat{\phi}_{\Lambda_{k}}\leq-k+M<0$
if $k>M$. Hence by Hopf’s lemma, we get
$\frac{\partial\tilde{\phi}_{k}}{\partial\nu}(x_{k})<0$, which contradicts
with the Neumann boundary condition. Thus, $\bar{\phi}_{k}$ attains its
minimum value at interior point $x_{k}\in\Omega$ for $k>M$, and we complete
the proof of Claim 3. ∎

Claim 3 implies $\nabla\bar{\phi}_{2,\Lambda_{k}}(x_{k})=0$ and
$\Delta\bar{\phi}_{2,\Lambda_{k}}(x_{k})\geq 0$ for $k>M$. Then as for (3.6),
we conclude that
$\displaystyle\min_{x\in\overline{\Omega}}\bar{\phi}_{2,\Lambda}(x)\geq-M_{7}$
for $\Lambda\geq 1$ which gives
$\displaystyle\lim_{x\in\overline{\Omega}}\phi_{\Lambda}(x)\geq-M_{6}$ for
$\Lambda\geq 1$, where $M_{6}$ and $M_{7}$ are positive constants independent
of $\Lambda$. Therefore the proof of Lemma 4.1 is complete.∎

The following is the proof of Theorem 1.2(ii). Because $\phi_{\Lambda}$ is the
solution of
$-\nabla\cdot(\epsilon\nabla\phi_{\Lambda})=\rho_{0}+f_{\Lambda}(\phi_{\Lambda})$
in $\Omega$ with the Neumann boundary condition
$\frac{\partial\phi_{\Lambda}}{\partial\nu}=0$ on $\partial\Omega$, we use the
$W^{2,p}$ estimate (cf. [2, Theorem 15.2]) to get

$\|\phi_{\Lambda}\|_{W^{2,p}(\Omega)}\leq
C\|\rho_{0}+f_{\Lambda}(\phi_{\Lambda})\|_{L^{p}(\Omega)}$

for all $p>1$, where $C$ is a positive constant independent of $\Lambda$. Then
as for Remark 3 of Section 3.1, we use Lemma 4.1 to get the uniform bound
estimate of $\phi_{\Lambda}$ in $W^{2,p}$ norm, and hence there exists a
sequence of function $\\{\phi_{\Lambda_{k}}\\}_{k=1}^{\infty}$ such that
$\phi_{\Lambda_{k}}$ converges to $\phi^{*}$ weakly in $W^{2,p}(\Omega)$. By
Remark 2 of Proposition 2.6, $\phi^{*}$ satisfies the equation (1.17) in weak
sense. Let $w_{k}=\phi_{\Lambda_{k}}-\phi^{*}$, $c_{i,k}:=c_{i,\Lambda_{k}}$,
and $f_{k}:=f_{\Lambda_{k}}$. Then by Sobolev embedding,
$w_{k}\in\mathcal{C}^{1,\alpha}(\Omega)$ for $\alpha\in(0,1)$, and
$\displaystyle\lim_{k\to\infty}\|w_{k}\|_{\mathcal{C}^{1,\alpha}(\Omega)}=0$.
Moreover, $w_{k}$ satisfies

$-\nabla\cdot(\epsilon\nabla
w_{k})=f_{k}(w_{k}+\phi^{*})-f^{*}(\phi^{*})\quad\text{in}~{}\Omega$

with the boundary condition $\frac{\partial w_{k}}{\partial\nu}=0$ on
$\partial\Omega$. Using the Schauder’s estimate (cf. [22, Theorem 6.30]) with
the mathematical induction, we get

$\displaystyle\|w_{k}\|_{\mathcal{C}^{m+2,\alpha}(\Omega)}$ $\displaystyle\leq
C\|f_{k}(w_{k}+\phi^{*})-f^{*}(\phi^{*})\|_{\mathcal{C}^{m,\alpha}(\Omega)}$
(4.5) $\displaystyle\leq
C^{\prime}\sum_{i=1}^{N}\left(\|c_{i,k}(w_{k}+\phi^{*})-c_{i}^{*}(w_{k}+\phi^{*})\|_{\mathcal{C}^{m,\alpha}(\Omega)}+\|c_{i}^{*}(w_{k}+\phi^{*})-c_{i}^{*}(\phi^{*})\|_{\mathcal{C}^{m,\alpha}(\Omega)}\right),$

for $m\in\mathbb{N}$ and $\alpha\in(0,1)$, where $C$ and $C^{\prime}$ are
positive constants independent of $k$. By Proposition 2.6 and induction
hypothesis
$\displaystyle\lim_{k\to\infty}\|w_{k}\|_{\mathcal{C}^{m,\alpha}(\Omega)}=0$,
we may use (4.4) to get
$\displaystyle\lim_{k\to\infty}\|w_{k}\|_{\mathcal{C}^{m+2,\alpha}(\Omega)}=0$,
i.e.
$\displaystyle\lim_{k\to\infty}\|\phi_{\Lambda_{k}}-\phi^{*}\|_{C^{m+2,\alpha}(\Omega)}=0$
for $m\in\mathbb{N}$ and $\alpha\in(0,1)$. Therefore, $\phi^{*}$ is the
solution of (1.17) with the Neumann boundary condition (1.9).

To complete the proof of Theorem 1.2(ii), we need to prove

######  Claim 4.

For any $m\in\mathbb{N}$, we have
$\displaystyle\lim_{\Lambda\to\infty}\|\phi_{\Lambda}-\phi^{*}\|_{\mathcal{C}^{m}(\Omega)}=0$.

###### Proof.

Suppose that there exist sequences $\\{\Lambda_{k}\\}$ and
$\\{\tilde{\Lambda}_{k}\\}$ tending to infinity such that sequences
$\\{\phi_{\Lambda_{k}}\\}$ and $\\{\phi_{\tilde{\Lambda}_{k}}\\}$ have limits
$\phi_{1}^{*}$ and $\phi_{2}^{*}$, respectively. It is clear that
$\phi_{1}^{*}$ and $\phi_{2}^{*}$ satisfy the equation (1.17) with the Neumann
boundary condition (1.9). Now we want to prove that
$\phi_{1}^{*}\equiv\phi_{2}^{*}$. Let $u=\phi_{1}^{*}-\phi_{2}^{*}$.
Subtracting the equation (1.17) with $\phi^{*}=\phi_{2}^{*}$ from that with
$\phi^{*}=\phi_{1}^{*}$, we obtain

$-\nabla\cdot(\epsilon\nabla
u)=f^{*}(\phi_{1}^{*})-f^{*}(\phi_{2}^{*})=c(x)u\quad\text{in}~{}\Omega,$

where function $c$ is defined by

$c(x)=\begin{cases}\displaystyle\frac{f^{*}(\phi_{1}^{*}(x))-f^{*}(\phi_{2}^{*}(x))}{\phi_{1}^{*}(x)-\phi_{2}^{*}(x)},&\text{if}~{}\phi_{1}^{*}(x)\neq\phi_{2}^{*}(x);\\\
\displaystyle\frac{\mathrm{d}f^{*}}{\mathrm{d}\phi}(\phi_{1}^{*}(x)),&\text{if}~{}\phi_{1}^{*}(x)=\phi_{2}^{*}(x).\end{cases}$

By Proposition 2.7, we have $c<0$ in $\Omega$. Here we have used the fact that
$\frac{f(\alpha)-f(\beta)}{\alpha-\beta}<0$ for $\alpha\neq\beta$ if $f$ is
strictly decreasing on $\mathbb{R}$. Since $\nabla\cdot(\epsilon\nabla
u)+c(x)u=0$ with $c<0$ in $\Omega$, then it is obvious that $u$ cannot be a
nonzero constant. Then by the strong maximum principle, we have that $u$
attains its nonnegative maximum value and nonnpositive minimum value at the
boundary point. Suppose $u$ has positive maximum value which is attained at
$x^{*}\in\partial\Omega$. Then by Hopf’s lemma, we get $\frac{\partial
u}{\partial\nu}(x^{*})>0$, which contradicts with the Neumann boundary
condition. Thus $u\leq 0$ in $\Omega$. Similarly, $u\geq 0$ in $\Omega$.
Therefore, $u\equiv 0$, i.e. $\phi^{*}_{1}\equiv\phi^{*}_{2}$ and we complete
the proof of Claim 4 and Theorem 1.2(ii).∎

##  5 Numerical results

Here we present numerical results to support Theorems 1.1 and 1.2. Throughout
this section, we assume the one-dimensional domain $\Omega=(-1,1)$, $N=3$,
$z_{0}=0$, $z_{1}=1$, $z_{2}=-1$, $z_{3}=2$, $\eta=\epsilon=0.1$, and
$\bar{\mu}_{i}=\tilde{\mu}_{0}=1$ for all $i=0,1,\dots,N$. By the method of
[20], we employ the Legendre–Gauss–Lobatto (LGL) points
$\\{x_{k}\\}_{k=0}^{L}$ as the partition of the interval $[-1,1]$ to discrtize
equations (1.11) and (1.12) as the following algebraic equations.

$\displaystyle\ln
c_{0,\Lambda}(x_{k})+\Lambda\lambda_{0}\sum_{j=0}^{N}\lambda_{j}(c_{0,\Lambda}(x_{k}))^{\lambda_{j}/\lambda_{0}}\exp(\bar{\mu}_{j}-z_{j}\phi_{\Lambda}(x_{k}))=\mu_{0},\quad\text{for}~{}k=0,\dots,L,$
(5.1)
$\displaystyle-\epsilon\text{D}_{L+1}^{2}\boldsymbol{\phi}_{\Lambda}=\boldsymbol{\rho}_{0}+\sum_{i=1}^{N}z_{i}{\bf
c}_{0,\Lambda}^{\lambda_{i}/\lambda_{0}}\circ\exp(\bar{\mu}_{i}-z_{i}\boldsymbol{\phi}_{\Lambda}),$
(5.2)

where
$\boldsymbol{\rho}_{0}=[\rho_{0}(x_{0}),\dots,\rho_{0}(x_{L})]^{\mathsf{T}}$,
${\bf
c}_{0,\Lambda}=[c_{0,\Lambda}(x_{0}),\dots,c_{0,\Lambda}(x_{L})]^{\mathsf{T}}$,
$\boldsymbol{\phi}_{\Lambda}=[\phi_{\Lambda}(x_{0}),\dots,\phi_{\Lambda}(x_{L})]^{\mathsf{T}}$
as column vectors, $\circ$ is the Hadamard product, and
$D_{L+1}=[d_{ij}]_{0\leq i,j\leq L}$ is the differentiation matrix which
satisfies
$D_{L+1}\boldsymbol{\phi}_{\Lambda}=[\phi_{\Lambda}^{\prime}(x_{0}),\dots,\phi_{\Lambda}^{\prime}(x_{L})]^{\mathsf{T}}$.
Besides, the Robin boundary condition (1.8) and the Neumann boundary condition
(1.9) are also discretized as

$\displaystyle(\text{Robin})\qquad~{}\phi(x_{0})-\eta\sum_{k=0}^{L}d_{0k}\phi(x_{k})=\phi_{bd}\left(-1\right),\quad\phi(x_{L})+\eta\sum_{k=0}^{L}d_{Lk}\phi(x_{k})=\phi_{bd}\left(1\right),$
(5.3)
$\displaystyle(\text{Neumann})\quad\sum_{k=0}^{L}d_{0k}\phi(x_{k})=0,\qquad\sum_{k=0}^{L}d_{Lk}\phi(x_{k})=0.$
(5.4)

For the Robin and Neumann boundary conditions, we replace the first and last
equations of (5.2) by equations (5.3) and(5.4), respectively. As for (5.1) and
(5.2), we discrtize the equations (1.15) and (1.17) as follows.

$\displaystyle\sum_{i=0}^{N}\lambda_{i}(c_{0}^{*}(x_{k}))^{\lambda_{i}/\lambda_{0}}\exp\left(\bar{\mu}_{i}-z_{i}\phi^{*}(x_{k})\right)=\tilde{\mu}_{0},\quad\text{for}~{}k=0,\dots,L,$
(5.5)
$\displaystyle-\epsilon\,\text{D}_{L+1}^{2}\boldsymbol{\phi}^{*}=\boldsymbol{\rho}_{0}+\sum_{i=1}^{N}z_{i}({\bf
c}_{0}^{*})^{\lambda_{i}/\lambda_{0}}\circ\exp\left(\bar{\mu}_{i}-z_{i}\boldsymbol{\phi}^{*}\right),$
(5.6)

where ${\bf c}_{0}^{*}=[c_{0}^{*}(x_{0}),\dots,c_{0}^{*}(x_{L})]^{\mathsf{T}}$
and
$\boldsymbol{\phi}^{*}=[\phi^{*}(x_{0}),\dots,\phi^{*}(x_{L})]^{\mathsf{T}}$.
Similarly, for the Robin and Neumann boundary conditions, we replace the first
and last equations of (5.6) by equations (5.3) and (5.4), respectively.

$\Lambda$ | (a) $\lambda_{i}=1$ for all $i$ | (b) $\lambda_{1}=2$, $\lambda_{2}=1.5$, $\lambda_{3}=1$
---|---|---
$1$ |  $5.9950$e-$01$ |  $6.4144$e-$01$
$10$ |  $1.4501$e-$01$ |  $1.3755$e-$01$
$10^{2}$ |  $1.7715$e-$02$ |  $1.5933$e-$02$
$10^{3}$ |  $1.8144$e-$03$ |  $1.6200$e-$03$
$10^{4}$ |  $1.8189$e-$04$ |  $1.6227$e-$04$
$10^{5}$ |  $1.8193$e-$05$ |  $1.6230$e-$05$

Table 1: The maximum norms $\left\|\phi_{\Lambda}-\phi^{*}\right\|_{\infty}$
with the Robin boundary condition (5.3), where $\eta=0.1$,
$\phi_{bd}\left(-1\right)=-2$ and $\phi_{bd}\left(1\right)=3$.

$\Lambda$ | (a) $\lambda_{i}=1$ for all $i$ | (b) $\lambda_{1}=2$, $\lambda_{2}=1.5$, $\lambda_{3}=1$
---|---|---
$1$ |  $6.2243$e-$01$ |  $6.4429$e-$01$
$10$ |  $1.4654$e-$01$ |  $1.3285$e-$01$
$10^{2}$ |  $1.7727$e-$02$ |  $1.5170$e-$02$
$10^{3}$ |  $1.8133$e-$03$ |  $1.5393$e-$03$
$10^{4}$ |  $1.8175$e-$04$ |  $1.5416$e-$04$
$10^{5}$ |  $1.8179$e-$05$ |  $1.5419$e-$05$

Table 2: The maximum norms $\left\|\phi_{\Lambda}-\phi^{*}\right\|_{\infty}$
with the Dirichlet boundary condition (5.3), where $\eta=0$,
$\phi_{bd}\left(-1\right)=-2$ and $\phi_{bd}\left(1\right)=3$.

$\Lambda$ | (a) $\lambda_{i}=1$ for all $i$ | (b) $\lambda_{1}=2$, $\lambda_{2}=1.5$, $\lambda_{3}=1$
---|---|---
$1$ |  $3.0027$e-$01$ |  $5.3706$e-$01$
$10$ |  $9.4618$e-$02$ |  $1.6046$e-$01$
$10^{2}$ |  $1.2565$e-$02$ |  $2.0737$e-$02$
$10^{3}$ |  $1.3010$e-$03$ |  $2.1385$e-$03$
$10^{4}$ |  $1.3057$e-$04$ |  $2.1452$e-$04$
$10^{5}$ |  $1.3062$e-$05$ |  $2.1459$e-$05$

Table 3: The maximum norms $\left\|\phi_{\Lambda}-\phi^{*}\right\|_{\infty}$
with the Neumann boundary condition (5.4), and $\rho_{0}(x)=x^{3}$.

Using the command fsolve in Matlab, we numerically solve the equations
(5.1)–(5.6) and have the following results. For various $\Lambda$ with
$\rho_{0}=0$ and the Robin boundary condition with $\eta=0.1$, the profiles of
$\phi_{\Lambda}$ and $\phi^{*}$ are represented in Figure 2, and the maximum
norms $\|\phi_{\Lambda}-\phi^{*}\|_{\infty}$ are listed in Table 3. For
various $\Lambda$ with $\rho_{0}=0$ and the Dirichlet boundary condition (the
Robin boundary condition with $\eta=0$), the profiles of $\phi_{\Lambda}$ and
$\phi^{*}$ are sketched in Figure 4, and the maximum norms
$\|\phi_{\Lambda}-\phi^{*}\|_{\infty}$ are shown in Table 3. These results are
consistent with Theorem 1.1. In addition, for the Neumann boundary condition,
we choose $\rho_{0}(x)=x^{3}$ to fulfill the constraint
$\int_{-1}^{1}\rho_{0}\mathrm{d}x=0$. For various $\Lambda$, the profiles of
$\phi_{\Lambda}$ and $\phi^{*}$ are presented in Figure 4 and the maximum norm
$\left\|\phi_{\Lambda}-\phi^{*}\right\|_{\infty}$ are expressed in Table 3
which support Theorem 1.2.

Figure 2: The numerical profiles of $\phi_{\Lambda}$ and $\phi^{*}$ with the
Robin boundary condition, where the curves 1–3 are the profiles of
$\phi_{\Lambda}$ with $\Lambda=0.5,2,8$ and the curve 4 is the profile of
$\phi^{*}$. In Figure 2 (a), $\lambda_{i}=1$ for all $i$. In Figure 2 (b),
$\lambda_{1}=2$, $\lambda_{2}=1.5$, $\lambda_{3}=1$.

Figure 3: The numerical profiles of $\phi_{\Lambda}$ and $\phi^{*}$ with the
Dirichlet boundary condition, where curves 1–3 are the profiles of
$\phi_{\Lambda}$ with $\Lambda=0.5,2,8$, and curve 4 is the profile of
$\phi^{*}$.

Figure 4: The numerical profiles of $\phi_{\Lambda}$ and $\phi^{*}$ with the
Neumann boundary condition, where curves 1–3 are the profiles of
$\phi_{\Lambda}$ with $\Lambda=0.5,2,8$, and curve 4 is the profile of
$\phi^{*}$.

Final Remark. (Poisson–Boltzmann approach)

Instead of steady state PNP-steric equations, we provide another approach to
obtain equations (1.6) and (1.7) which can be derived by the following energy
functional.

$E[c_{0},c_{1},\dots,c_{N},\phi]=E_{PB}[c_{0},c_{1},\dots,c_{N},\phi]+E_{LJ}[c_{0},c_{1},\dots,c_{N}],$

where

$\displaystyle
E_{PB}[c_{0},c_{1},\cdots,c_{N},\phi]=\int_{\Omega}\left[-\frac{1}{2}\epsilon\left|\nabla\phi\right|^{2}+\sum_{i=0}^{N}c_{i}(\ln
c_{i}-1)+\left(\rho_{0}+\sum_{i=1}^{N}z_{i}c_{i}\right)\phi\right]\mathrm{d}x,$
$\displaystyle
E_{LJ}[c_{0},c_{1},\dots,c_{N}]=\frac{1}{2}\sum_{i,j=0}^{N}g_{ij}\int_{\Omega}c_{i}(x)c_{j}(x)\mathrm{d}x.$

Notice that $E_{PB}$ is the energy functional of conventional
Poisson–Boltzmann equation with the form
$-\nabla\cdot(\epsilon\nabla\phi)=\rho_{0}+\sum_{i=1}^{N}z_{i}e^{\mu_{i}-z_{i}\phi}$
which can be obtained by $\delta E_{PB}/\delta\phi=0$ and $\delta
E_{PB}/\delta c_{i}=\mu_{i}$ for $i=0,1,\dots,N$. Besides, $E_{LJ}$ is the
energy functional of the approximate Lennard–Jones potentials (cf. [35]), and
one may derive the equations (1.6) and (1.7) by $\delta E/\delta\phi=0$ and
$\delta E/\delta c_{i}=\mu_{i}$ for $i=0,1,\dots,N$.

Conclusion

Modified Poisson-Boltzmann (mPB) equations play an important role to
understand the steric effects of ion and solvent molecules. To get such
equations, we study steady state Poisson–Nernst–Planck equations with steric
effects (PNP-steric equatons) and derive the Poisson–Boltzmann equation with
steric effects (PB-steric equation) by the assumptions of steric effects.
Under the Robin (or Neumann) boundary condition, the PB-steric equation has a
unique solution $\phi_{\Lambda}$, a parameter $\Lambda$, and positive
constants $\lambda_{i}$’s which depend on the radii of ions and solvent
moleclues. As $\Lambda$ goes to infinity, the limiting equation of PB-steric
equation becomes a mPB equation. We firstly use the implicit function theorem
on Banach spaces to prove the convergence of nonlinear term of PB-steric
equation as $\Lambda$ tends to infinity. Then we prove the uniform bound
estimate of solution $\phi_{\Lambda}$ and obtain its convergence. Numerical
simulations are provided to support theoretical results. Further work is
needed on the application of PB-steric equation.

Acknowledgement

The research of T.C. Lin is partially supported by the National Center for
Theoretical Sciences (NCTS) and MOST grant 109-2115-M-002 -003 of Taiwan.

##  6 Appendix I. Existence and uniqueness of $\phi_{\Lambda}$

In this section, we fix $\Lambda>0$ and prove the existence and uniqueness of
solution $\phi_{\Lambda}$ of equation (1.12) with the Robin boundary condition
and Neumann boundary condition (1.8) and (1.9), respectively. For the
existence of solution $\phi_{\Lambda}$, we study the following energy
minimization problems (N) and (R) for the Neumann and Robin boundary
conditions, respectively.

  * (N)

Minimize $E_{eq}[\phi]$ subject to $\phi\in H^{1}(\Omega)$,

  * (R)

Minimize $E[\phi]:=E_{eq}[\phi]+B_{\eta}[\phi]$ subject to
$\phi\in\mathbb{H}_{\eta}(\Omega)$,

where the functionals

$E_{eq}[\phi]=\frac{1}{2}\int_{\Omega}\epsilon|\nabla\phi|^{2}\,\mathrm{d}x-\int_{\Omega}\rho_{0}\phi\,\mathrm{d}x-\int_{\Omega}F_{\Lambda}(\phi)\,\mathrm{d}x,\qquad
B_{\eta}[\phi]=\begin{cases}\frac{1}{2\eta}\int_{\partial\Omega}\epsilon(\phi-\phi_{bd})^{2}\,\mathrm{d}S_{x}&\text{if}~{}\eta>0,\\\
0&\text{if}~{}\eta=0,\end{cases}$

and defined on the space

$\mathbb{H}_{\eta}(\Omega):=\begin{cases}H^{1}(\Omega)&\text{if}~{}\eta>0,\\\
\\{u\in H^{1}(\Omega):u-\phi_{bd}\in
H_{0}^{1}(\Omega)\\}&\text{if}~{}\eta=0.\end{cases}$

Here the function
$F_{\Lambda}(\phi)=\int_{0}^{\phi}f_{\Lambda}(s)\,\mathrm{d}s$, and the
function $f_{\Lambda}=f_{\Lambda}(\phi)$ is defined in Proposition 2.2 which
gives $F_{\Lambda}^{\prime\prime}(\phi)=f_{\Lambda}^{\prime}(\phi)<0$ for
$\phi\in\mathbb{R}$. Hence we may apply the Direct method (cf. [40]) to solve
problems (N) and (R). The argument of problem (N) is similar to that of
problem (R) so we omit here and only provide the argument of problem (R) in
the rest of this section.

To apply the Direct method on problem (R), we need the following lemma.

######  Lemma 6.1.

Functional $E$ is coercive on $H^{1}(\Omega)$ for $\eta>0$.

###### Proof.

By Young’s inequality, we have

$B_{\eta}\left[\phi\right]\geq\frac{1}{4\eta}\int_{\partial\Omega}\epsilon\phi^{2}\,\mathrm{d}S_{x}-\frac{1}{2\eta}\int_{\partial\Omega}\epsilon\phi_{bd}^{2}\,\mathrm{d}S_{x}.$
(6.1)

By Propositions 2.2 and 2.4, function $F_{\Lambda}$ is strictly concave and
has an absolute maximum denoted by $M_{F}$. Hence by (6.1), we obtain

$E[\phi]\geq
C_{\eta}\left(\int_{\Omega}\left|\nabla\phi\right|^{2}\mathrm{d}x+\int_{\partial\Omega}\phi^{2}\,\mathrm{d}S_{x}\right)-\int_{\Omega}|\rho_{0}\phi|\,\mathrm{d}x-\frac{1}{2\eta}\int_{\partial\Omega}\epsilon\phi_{bd}^{2}\,\mathrm{d}S_{x}-M_{F}\left|\Omega\right|,$
(6.2)

where
$C_{\eta}=\frac{1}{4\eta}\min\left\\{1,2\eta\right\\}\min_{\overline{\Omega}}\epsilon>0$
and $\left|\Omega\right|$ denotes the Lebesuge measure of $\Omega$. Moreover,
for any $\phi\in H^{1}(\Omega)$, Friedrichs’ inequality gives

$\int_{\Omega}\left|\nabla\phi\right|^{2}\mathrm{d}x+\int_{\partial\Omega}\phi^{2}\,\mathrm{d}S_{x}\geq
C_{1}\int_{\Omega}\phi^{2}\mathrm{d}x,$ (6.3)

where $C_{1}$ is a positive constant depending only on the dimension $d$ and
the measures of $\Omega$ and $\partial\Omega$. Besides, Cauchy–Schwarz
inequality gives

$\int_{\Omega}\left|\rho_{0}\phi\right|\,\mathrm{d}x\leq\left(\int_{\Omega}\rho_{0}^{2}\,\mathrm{d}x\right)^{1/2}\left(\int_{\Omega}\phi^{2}\,\mathrm{d}x\right)^{1/2}$
(6.4)

and by (6.2)–(6.4), Lemma 6.1 follows. ∎

By Lemma 6.1, we may prove the existence of minimizer of problem (R) as
follows.

######  Proposition 6.2.

Functional $E$ has a minimizer $\phi\in\mathbb{H}_{\eta}$ for any $\eta\geq
0$.

###### Proof.

Suppose $\eta>0$ and $\mathbb{H}_{\eta}(\Omega)=H^{1}(\Omega)$. By
(6.2)–(6.4), $\displaystyle\inf_{\phi\in H^{1}(\Omega)}E[\phi]$ exists so
there exists a minimizing sequence $\\{\phi_{n}\\}_{n=1}^{\infty}\subsetneq
H^{1}(\Omega)$ such that

$\lim_{n\to\infty}E[\phi_{n}]=\inf_{\phi\in H^{1}(\Omega)}E[\phi]:=m_{E}.$

Due to the coerciveness of $E$ (see Lemma 6.1), we have
$\displaystyle\sup_{n\in\mathbb{N}}\left\|\phi_{n}\right\|_{H^{1}(\Omega)}<\infty$.
Along with (6.2), we can get
$\displaystyle\sup_{n\in\mathbb{N}}\left\|\phi_{n}\right\|_{L^{2}(\partial\Omega)}<\infty$.
Thus, there exists a subsequence $\\{\phi_{n_{k}}\\}$ of $\\{\phi_{n}\\}$ and
$\phi\in H^{1}(\Omega)$ such that $\phi_{n_{k}}\rightharpoonup\phi$ weakly in
$H^{1}(\Omega)$ and $\phi_{n_{k}}\rightharpoonup\Gamma\phi$ weakly in
$L^{2}(\partial\Omega)$ as $k\to\infty$, where $\Gamma\phi$ is the trace of
$\phi$ on $\partial\Omega$. Since $\phi_{n_{k}}\rightharpoonup\phi$ weakly in
$H^{1}(\Omega)$ implies $\nabla\phi_{n_{k}}\rightharpoonup\nabla\phi$ weakly
in $L^{2}(\Omega)$, then we obtain

$\liminf_{k\to\infty}\int_{\Omega}\left|\nabla\phi_{n_{k}}\right|^{2}\mathrm{d}x\geq\int_{\Omega}\left|\nabla\phi\right|^{2}\mathrm{d}x,\quad\liminf_{k\to\infty}\int_{\partial\Omega}\left|\phi_{n_{k}}-\phi_{bd}\right|^{2}\mathrm{d}S_{x}\geq\int_{\partial\Omega}\left|\Gamma\phi-\phi_{bd}\right|^{2}\mathrm{d}S_{x},$

and

$\lim_{k\to\infty}\phi_{n_{k}}=\phi~{}~{}\text{a.e. in}~{}\Omega.$

Along with the Fatou’s lemma, we have

$\liminf_{k\to\infty}\int_{\Omega}-F_{\Lambda}(\phi_{n_{k}})\,\mathrm{d}x\geq\int_{\Omega}-F_{\Lambda}(\phi)\,\mathrm{d}x,$

which gives

$m_{E}=\lim_{k\to\infty}E[\phi_{n_{k}}]\geq E[\phi]\geq m_{E},$

and $E$ attains its minimum $m_{E}$ at $\phi\in
H^{1}(\Omega)=\mathbb{H}_{\eta}(\Omega)$ for $\eta>0$. As $\eta=0$ (i.e., the
Dirichlet boundary condition), $B_{\eta}[\phi]=0$ so we may use the similar
argument to prove the existence of minimizer and complete the proof.∎

The minimizer of Proposition 6.2 satisfies equation (1.12) with the Robin
boundary condition (1.8) in weak sense. Hereafter, for simplicity, we denote
the solution as $\phi$ instead of $\phi_{\Lambda}$ before. Now we prove the
regularity of solution $\phi$ as follows. Because $\phi$ is a minimizer of the
functional $E$ on $H^{1}(\Omega)$, $\phi$ satisfies

$\int_{\Omega}\left[\epsilon\nabla\phi\cdot\nabla
v-v\rho_{0}-vf_{\Lambda}(\phi)\right]\mathrm{d}x+\frac{1}{\eta}\int_{\partial\Omega}\epsilon\left(\phi-\phi_{bd}\right)v\,\mathrm{d}S_{x}=0$
(6.5)

for any $v\in H^{1}(\Omega)$. Let $\Phi$ be the solution of the auxiliary
Poisson equation

$-\nabla\cdot\left(\epsilon\nabla\Phi\right)=\rho_{0}+f_{\Lambda}(\phi)\quad\text{in}~{}\Omega$
(6.6)

with the Robin boundary condition
$\Phi+\eta\frac{\partial\Phi}{\partial\nu}=\phi_{bd}$ on $\partial\Omega$. We
multiply equation (6.6) by $(\Phi-\phi)$ and integrate it over $\Omega$. Then
we may use integration by parts to obtain

$\int_{\Omega}\epsilon\nabla\Phi\cdot\nabla(\Phi-\phi)\,\mathrm{d}x-\int_{\Omega}(\Phi-\phi)\left[\rho_{0}+f_{\Lambda}(\phi)\right]\mathrm{d}x+\frac{1}{\eta}\int_{\partial\Omega}\epsilon(\Phi-\phi_{bd})(\Phi-\phi)\,\mathrm{d}S_{x}=0.$
(6.7)

Here we have used the fact
$\frac{\partial\Phi}{\partial\nu}=\frac{1}{\eta}(\phi_{bd}-\Phi)$ on
$\partial\Omega$. Besides, we subtract (6.5) with $v=\Phi-\phi$ from (6.7) and
get

$\int_{\Omega}\epsilon\left|\nabla\left(\Phi-\phi\right)\right|^{2}\mathrm{d}x+\frac{1}{\eta}\int_{\partial\Omega}\epsilon(\Phi-\phi)^{2}\mathrm{d}S_{x}=0,$

which implies that $\Phi\equiv\phi$ a.e. in $\Omega$. Thus, by a standard
bootstrap argument on equation (6.6), solution $\phi$ becomes a classical
solution. On the other hand, as $\eta=0$, we can use a similar argument to
prove the regularity of solution $\phi$.

Now we prove the uniqueness of equation (1.12) with the Robin boundary
condition (1.8). Suppose that
$\phi_{1},\phi_{2}\in\mathcal{C}^{\infty}\left(\Omega\right)\cap\mathcal{C}^{2}\left(\overline{\Omega}\right)$
are solutions of equation (1.12) with the Robin boundary condition (1.8). We
subtract equation (1.12) with $\phi=\phi_{2}$ from that with $\phi=\phi_{1}$.
Then we have

$-\nabla\cdot(\epsilon\nabla(\phi_{1}-\phi_{2}))=f_{\Lambda}(\phi_{1})-f_{\Lambda}(\phi_{2})\quad\text{in}~{}\Omega.$
(6.8)

Multiplying the equation (6.8) by $\phi_{1}-\phi_{2}$ and integrating it over
$\Omega$, we get

$-\int_{\partial\Omega}(\phi_{1}-\phi_{2})\nabla\cdot(\epsilon\nabla(\phi_{1}-\phi_{2}))\mathrm{d}x=\int_{\Omega}\left(\phi_{1}-\phi_{2}\right)[f_{\Lambda}(\phi_{1})-f_{\Lambda}(\phi_{2})]\,\mathrm{d}x.$

Then using integration by parts, we obtain

$\int_{\Omega}\epsilon|\nabla(\phi_{1}-\phi_{2})|^{2}\mathrm{d}x+\eta\in_{\partial\Omega}\epsilon\left|\frac{\partial(\phi_{1}-\phi_{2})}{\partial\nu}\right|^{2}\mathrm{d}S_{x}=\int_{\Omega}(\phi_{1}-\phi_{2})[f_{\Lambda}(\phi_{1})-f_{\Lambda}(\phi_{2})]\mathrm{d}x.$
(6.9)

Here we have used the fact
$\phi_{1}-\phi_{2}=-\eta\frac{\partial(\phi_{1}-\phi_{2})}{\partial\nu}$ on
$\partial\Omega$ which comes from the Robin boundary condition (1.8). By
Proposition 2.2,

$\left(\phi_{1}-\phi_{2}\right)[f_{\Lambda}(\phi_{1})-f_{\Lambda}(\phi_{2})]\leq
0\quad\text{in}~{}\Omega.$ (6.10)

Thus, by (6.9) and (6.10), we have

$\int_{\Omega}\epsilon\left|\nabla\left(\phi_{1}-\phi_{2}\right)\right|^{2}\mathrm{d}x+\eta\int_{\partial\Omega}\epsilon\left|\frac{\partial\left(\phi_{1}-\phi_{2}\right)}{\partial\nu}\right|^{2}\mathrm{d}S_{x}\leq
0,$

which implies that $\phi_{1}=\phi_{2}$ in $\overline{\Omega}$ so we complete
the proof of uniqueness and conclude the following theorem.

######  Theorem 6.3.

Under the same hypotheses in Theorem 1.1, equation (1.12) with the Robin
boundary condition (1.8) has a unique solution
$\phi_{\Lambda}\in\mathcal{C}^{\infty}\left(\Omega\right)\cap\mathcal{C}^{2}\left(\overline{\Omega}\right)$
for $\Lambda>0$.

##  7 Appendix II.

In this section, we prove that equations (1.14)–(1.17) and equations
(1.19)–(1.21) have the same form (up to scalar multiples). Firstly, we recall
these equations as follows.

$\displaystyle\sum_{i=0}^{N}\lambda_{i}c_{i}^{*}(\phi^{*})=\tilde{\mu}_{0},$
(7.1) $\displaystyle
c_{i}^{*}(\phi^{*})=(c_{0}^{*}(\phi^{*}))^{\lambda_{i}/\lambda_{0}}\exp(\bar{\mu}_{i}-z_{i}\phi^{*})\quad\text{for}~{}i=1,\dots,N,$
(7.2)
$\displaystyle-\nabla\cdot(\epsilon\nabla\phi^{*})=\rho_{0}+\sum_{i=1}^{N}z_{i}c_{i}^{*}(\phi^{*})\quad\text{in}~{}\Omega,$
(7.3)

and

$\displaystyle\sum_{i=0}^{N}v_{i}c_{i}(\psi)=1,$ (7.4)
$\displaystyle\frac{v_{i}}{v_{0}}\ln(v_{0}c_{0})-\ln(v_{i}c_{i})=\beta
q_{i}\psi-\beta\mu_{i}\quad\text{for}~{}i=1,\dots,N,$ (7.5)
$\displaystyle-\nabla\cdot(\epsilon\nabla\psi)=\rho_{0}+\sum_{i=1}^{N}\beta
q_{i}c_{i}(\psi)\quad\text{in}~{}\Omega.$ (7.6)

From (7.1) and (7.4), we may assume that $\tilde{\mu}_{0}=1$. Then we have
$\lambda_{i}=v_{i}$, $c_{i}^{*}=c_{i}$ and $\phi^{*}=\psi$ for
$i=0,1,\dots,N$. Moreover, by taking the natural logarithm on both sides of
(7.2), we have

$\frac{\lambda_{i}}{\lambda_{0}}\ln c_{0}^{*}(\phi^{*})-\ln
c_{i}^{*}(\phi^{*})=z_{i}\phi^{*}-\bar{\mu}_{i}\quad\text{for}~{}i=1,\dots,N.$
(7.7)

Comparing (7.5) and (7.7) with $\phi^{*}=\psi^{*}$, we get

$z_{i}=\beta q_{i},\qquad\bar{\mu}_{i}=\beta\mu_{i}+\frac{v_{i}}{v_{0}}\ln
v_{0}-\ln v_{i}\quad\text{for}~{}i=1,\dots,N.$

Hence it is clear that (7.3) and (7.6) are identical. Therefore, equations
(7.1)–(7.3) and equations (7.4)–(7.6) have the same form (up to scalar
multiples).

## References

  * [1] Z. Abbas, M. Gunnarsson, E. Ahlberg, S. Nordholm, Corrected Debye–Hückel Theory of Salt Solutions: Size Asymmetry and Effective Diameters, J. Phys. Chem. B 106 (2002), 1403–1420.
  * [2] S. Agmon, A. Douglis, L. Nirenberg, Estimates Near the Boundary for Solutions of Elliptic Partial Differential Equations Satisfying General Boundary Conditions. I, Commun. Pure Appl. 12 (1959), 623–727.
  * [3] R. Aitbayev, P. W. Bates, H. Lu, L. Zhang, M. Zhang, Mathematical studies of Poisson–Nernst–Planck model for membrance channels: Finite ion size effects without electroneutrality boundary conditions, J. Comput. Appl. Math. 362 (2019), 510–527.
  * [4] D. Andelman, Electrostatic Properties of Membranes: The Poisson–Boltzmann Theory, Handbook of Biological Physics, R. Lipowsky and E. Sackmann 1 (1995), 603–642.
  * [5] V. Barcilon, D.-P. Chen, R. S. Eisenberg, J. W. Jerome, Qualitative Properties of Steady-State Poisson–Nernst–Planck Systems: Perturbation and Simulation Study, SIAM J. Appl. Math. 57 (1997), 631–648.
  * [6] M. Z. Bazant, B. D. Storey, A. A. Kornshev, Double Layer in Ionic Liquids: Overscreening versus Crowding, Phys. Rev. Lett. 106 (2011), 046102.
  * [7] P. W. Bates, W. Liu, H. Lu, M. Zhang, Ion Size and Valence Effects on Ionic Flows via Poisson–Nernst–Planck Models, Commun. Math. Sci. 15 (2017), 881–901.
  * [8] G. M.-Baude, B. Stadlbauer, S. Howorka, C. Heitzinger Protein Transport through Nanopores Illuminated by Long-Time-Scale Simulations, ACS Nano 15 (2021), 9900–9912.
  * [9] J. J. Bikerman, XXXIX. Structure and capacity of electrical double layer, Philos. Mag. 7 (1942), 384–397.
  * [10] I. Borukhov, D. Andelman, Steric Effects in Electrolytes: A Modified Poisson–Boltzmann Equation, Phys. Rev. Lett. 79 (1997), 435–438.
  * [11] D. Chen, A New Poisson–Nernst–Planck Model with Ion-Water Interactions for Charge Transport in Ion Channels, Bull. Math. Biol. 78 (2016), 1703–1726.
  * [12] R. D. Coalson, M. G. Kurnikova, Poisson–Nernst–Planck Theory Approach to the Calculation of Current through Biological Ion Channels, IEEE Trans. Nanobioscience 4 (2005), 81–93.
  * [13] K. Deimling, Nonlinear Functional Analysis, (1973), Springer–Verlag.
  * [14] J. Ding, Z. Wang, S. Zhou, Positivity preserving finite difference methods for Poisson–Nernst–Planck equations with steric interactions: Application to slit-shaped nanopore conductance, J. Comput. Phys. 397 (2019), 108864.
  * [15] J. Ding, Z. Wang, S. Zhou, Structure-preserving and efficient numerical methods for ion transport, J. Comput. Phys. 418 (2020), 109597.
  * [16] J. Ding, H. Sun, S. Zhou Hysteresis and linear stability analysis on multiple steady-state solutions to the Poisson–Nernst–Planck equations with steric interactions, Phys. Rev. E 102 (2020), 053301.
  * [17] R. Eisenberg, D. Chen, Poisson–Nernst–Planck (PNP) theory of an open ionic channel, Biophys. J., 64(2 Pt. 2) (1993), A22.
  * [18] B. Eisenberg, W. Liu, Poisson–Nernst–Planck Systems for Ion Channels with Permanent Charges, SIAM J. Math. Anal., 38 (2007), 1932–1966.
  * [19] B. Eisenberg, Crowded Charges in Ion Channels. In Advanced in Chemical Physics 148 (2012); S. A. Rice, A. R. Dinner, Eds.; John Wiley & Sons: Hoboken, N. J., 77–223.
  * [20] G. I. El–Baghdady, M. S. El–Azab, W. S. El–Beshbeshy, Legendre–Gauss–Lobatto Pseudo-spectral Method for One-Dimensional Advection-Diffusion Equation, Sohag J. Math. 2 (2015), 29–35.
  * [21] N. Gavish, Poisson–Nernst–Planck equations with steric effects—non-convexity and multiple stationary solutions, Physica D 368 (2018), 50–65.
  * [22] D. Gilbarg, N. S. Trudinger, Elliptic Partial Differential Equations of Second Order, (1977), Springer–Verlag.
  * [23] P. Grochowski, J. Trylska, Continuum Molecular Electrostatics, Salt Effects, and Counterion Binding—A Review of the Poisson–Boltzmann Theory and Its Modifications, Biopolymers 89 (2007), 93–113.
  * [24] T.-L. Horng, T.-C. Lin, C. Liu, B. Eisenberg, PNP Equations with Steric Effects: A Model of Ion Flow through Channels, J. Phys. Chem. B 116 (2012), 11422–11441.
  * [25] T.-L. Horng, R. S. Eisenberg, C. Liu, F. Bezanilla, Continuum Gating Current Models Computed with Consistent Interactions, Biophys. J. 116 (2019), 270–282.
  * [26] B. Huang, S. Maset, K. Bohinc, Interaction between Charged Cylinders in Electrolyte Solution; Excluded Volume Effect, J. Phys. Chem. B 121 (2017), 9013–9023.
  * [27] I. Kalcher, J. C. F. Schulz, J. Dzubiella, Ion-Specific Excluded-Volume Correlations and Solvation Forces, Phys. Rev. Lett., 104 (2010), 097802.
  * [28] M. S. Kilic, M. Z. Bazant, Steric effects in the dynamics of electrolytes at large applied voltages. I. Double-layer charging, Phys. Rev. E 75 (2007), 021502.
  * [29] C. Köhn, D. v. Laethem, J. Deconinck, A. Hubin, A simulation study of steric effects on the anodic dissolution at high current densities, Materials and Corrosion 72 (2020), 610–619.
  * [30] A. A. Kornyshev, Double-Layer in Ioninc Liquids: Paradigm Change, J. Phys. Chem. B 111 (2007), 5545–5557.
  * [31] S. G. Krantz, H. R. Parks, The Implicit Function Theorem: History, Theory, and Applications, (2002), Birkhäuser.
  * [32] B. Li, Continuum electrostatics for ionic solutions with non-uniform ionic sizes, Nonlinearity 22 (2009), 811–833.
  * [33] B. Li, Minimization of electrostatic free energy and the Poisson–Boltzmann equation for molecular solvation with implicit solvent, SIAM J. Math. Anal. 40 (2009), 2536–2566.
  * [34] B. Li, P. Liu, Z. Xu, S. Zhou, Ionic size effects: generalized Boltzmann distributions, counterion stratification and modified Debye length, Nonlinearity 26 (2013), 2899–2922.
  * [35] T.-C. Lin, B. Eisenberg, A new approach to the Lennard–Jones potential and a new model: PNP-steric equations, Commun. Math. 12 (2014), 149–173.
  * [36] T.-C. Lin, B. Eisenberg, Multiple solutions of steady-state Poisson–Nernst–Planck equations with steric effects, Nonlinearity 28 (2015), 2053–2080.
  * [37] B. Lu, Y. C. Zhou, Poisson–Nernst–Planck Equations for Simulating Biomolecular Diffusion-Reaction Processes II: Size Effects on Ionic Distributions and Diffusion-Reaction Rates, Biophys. J. 100 (2011), 2475–2485.
  * [38] Y. S. Park, I. S. Kang, Multi-ionic effects on the equilibrium and dynamic properties of electric double layers based on the Bikerman correction, J. Electroanal. Chem. 880 (2021), 114923.
  * [39] F. Siddiqua, Z. Wang, S. Zhou, A modified Poisson–Nernest–Planck model with excluded volume effect: theory and numerical implementation, Commun. Math. Sci. 16 (2018), 251–271.
  * [40] M. Struwe, Variational Methods: Applications to Nonlinear Partial Differential Equations and Hamiltonian Systems, (2008), Springer–Verlag.