# PB-steric equations: A general model of Poisson–Boltzmann equations

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 (which can
produce oscillations in charge density profiles) becomes significant and the
conventional Poisson–Boltzmann (PB) equation should be modified. Several
modified PB equations were developed but the associated total ionic charge
density has no oscillation. This motivates us to derive a general model of PB
equations called the PB-steric equations with a parameter $\Lambda$, which not
only include the conventional and modified PB equations but also have
oscillatory total ionic charge density under different assumptions of steric
effects and chemical potentials. As $\Lambda=0$, the PB-steric equation
becomes the conventional PB equation, but as $\Lambda>0$, the concentrations
of ions and solvent molecules are determined by the Lambert type functions. To
approach the modified PB equations, we study the asymptotic limit of PB-steric
equations with the Robin boundary condition as $\Lambda$ goes to infinity. Our
theoretical results show that the PB-steric equations (for
$0\leq\Lambda\leq\infty$) may include the conventional and modified PB
equations. On the other hand, we use the PB-steric equations to find
oscillatory total ionic charge density which cannot be obtained in the
conventional and modified PB equations.
Key words. PB-steric equations, total ionic charge density, oscillation
AMS subject classifications. 35C20, 35J60, 35Q92

##  1 Introduction

Understanding the distribution of ions in the electrolyte is one of the most
crucial problems in many physical and electrochemical fields. As a well-known
mathematical model, the Poisson–Boltzmann (PB) equations play a central role
in the study of ionic distribution. The conventional PB equations [24, 31]
treat ions as point particles without size and can be denoted as

$\displaystyle-\nabla\cdot(\varepsilon\nabla\phi)=4\pi\rho_{0}+4\pi\sum_{i=1}^{N}z_{i}e_{0}c_{i}(\phi),$
(1.1) $\displaystyle
k_{B}T\ln(c_{i})+z_{i}e_{0}\phi=\mu_{i}\quad\text{for}~{}i=1,\dots,N,$ (1.2)

where $\phi$ is the electrostatic potential, $\varepsilon$ is the dielectric
function, and $\rho_{0}$ is the permanent charge density function. In
addition, $N$ is the number of ion species, $c_{i}$ is the concentration of
$i$th ion species with the valence $z_{i}\neq 0$ for $i=1,\dots,N$. Moreover,
$c_{i}=c_{i}^{\text{b}}\exp(-z_{i}e_{0}\phi/k_{B}T)$, where $c_{i}^{\text{b}}$
is the concentration of the $i$th ion species in the bulk, $k_{B}$ is the
Boltzmann constant, $T$ is the absolute temperature, $e_{0}$ is the elementary
charge, and $\mu_{i}=k_{B}T\ln(c_{i}^{\text{b}})$ is the chemical potential.
However, when ions are crowded, steric repulsions may appear due to ion sizes,
so the conventional PB equations should be modified (cf. [8, 9, 20, 21]).

Under the hypothesis of volume exclusion, one may study the case of two-
species ions with the same radius and obtain the following modified PB
equations (cf. [4, 5, 6, 17]):

$\displaystyle-\nabla\cdot(\varepsilon\nabla\phi)=4\pi(z_{1}e_{0}c_{1}+z_{2}e_{0}c_{2}),$
(1.3) $\displaystyle
c_{i}=c_{i}^{\text{b}}\frac{\exp(-z_{i}e_{0}\phi/k_{B}T)}{\displaystyle
1-\gamma+\frac{\gamma}{z_{1}-z_{2}}(c_{1}^{\text{b}}\exp(-z_{1}e_{0}\phi/k_{B}T)+c_{2}^{\text{b}}\exp(-z_{2}e_{0}\phi/k_{B}T))},$
(1.4)

for $i=1,2$, where $\gamma$ is the total bulk volume fraction of ions, and
$z_{1}>0$ and $z_{2}<0$ are the valence of cations and anions, respectively.
Moreover, the concentrations of ion species in the bulk are given by
$c_{1}^{\text{b}}=-z_{2}$ and $c_{2}^{\text{b}}=z_{1}$. For the background and
development of (1.3)–(1.4), we refer the interested reader to [14, 22].
Moreover, when ions and solvent molecules may have different sizes, the
associated PB equations become

$\displaystyle-\nabla\cdot(\varepsilon\nabla\phi)=4\pi\rho_{0}+4\pi\sum_{i=1}^{N}z_{i}e_{0}c_{i}(\phi),$
(1.5) $\displaystyle
k_{B}T\ln(c_{i})-k_{B}T\frac{v_{i}}{v_{0}}\ln(c_{0})+z_{i}e_{0}\phi=\bar{\mu}_{i}\quad\text{for}~{}i=1,\dots,N,$
(1.6) $\displaystyle\sum_{i=0}^{N}v_{i}c_{i}=1,$ (1.7)

where $\displaystyle c_{0}=\frac{1}{v_{0}}(1-\sum\limits_{i=1}^{N}v_{i}c_{i})$
is the concentration of solvent molecules with the volume $v_{0}$ and the
valence $z_{0}=0$, $v_{i}$ is the volume of $i$th ion species, and
$\bar{\mu}_{i}$ is the associated chemical potential for $i=1,\dots,N$ (cf.
[27, 28, 29, 32]). Equations (1.1)–(1.7) can be denoted as

$\displaystyle-\nabla\cdot(\varepsilon\nabla\phi)=4\pi\rho_{0}+4\pi\sum_{i=1}^{N}z_{i}e_{0}c_{i},$
(1.8)
$\displaystyle\mu_{i}=k_{B}T\ln(c_{i})+z_{i}e_{0}\phi+\mu_{i}^{\text{ex}}+U_{i}^{\text{wall}}\quad\text{for}~{}i=0,1,\dots,N,$
(1.9)

with different $\mu_{i}$ and $\mu_{i}^{\text{ex}}$, where $\mu_{i}$ is the
chemical potential, $\mu_{i}^{\text{ex}}$ is the excess chemical potential
which describes the interaction potential of ions and solvent molecules.
Besides, $U_{i}^{\text{wall}}$ is the potential which comes from the
interactions of ions and solvent molecules with the wall. Under
$U_{i}^{\text{wall}}=0$ and conditions of $\mu_{i}$ and $\mu_{i}^{\text{ex}}$,
equations (1.8)–(1.9) become the following equations.

  * •

(1.1)–(1.2) if $\mu_{i}=k_{B}T\ln(c_{i}^{\text{b}})$ and
$\mu_{i}^{\text{ex}}=0$ for $i=1,\cdots,N$.

  * •

(1.3)–(1.4) if $\mu_{i}=k_{B}T\ln(c_{i}^{\text{b}})$ and
$\displaystyle\mu_{i}^{\text{ex}}=k_{B}T\ln\left(1-\gamma+\frac{\gamma}{z_{1}-z_{2}}(c_{1}^{\text{b}}e^{-z_{1}e_{0}\phi/k_{B}T}+c_{2}^{\text{b}}e^{-z_{2}e_{0}\phi/k_{B}T})\right)$
for $i=1,2$.

  * •

(1.5)–(1.7) if $\mu_{i}=\bar{\mu}_{i}$,
$\mu_{i}^{\text{ex}}=-k_{B}T\frac{v_{i}}{v_{0}}\ln(c_{0})$ for $i=1,\cdots,N$.

The steric interactions of ions can produce oscillations in concentration
(density) profiles but such oscillations cannot be obtained by (1.1)–(1.4)
(cf. [13, 18, 36, 38]). In [40], the oscillatory concentrations can be found
in (1.5)–(1.7), so it is natural to ask if there exists oscillatory total
ionic charge density. However, the total ionic charge density
$\sum\limits_{i=1}^{N}z_{i}c_{i}(\phi)$ of (1.1)–(1.7) is decreasing to $\phi$
(see [29] and Propositions 2.6 and 2.10 below). This motivates us to derive a
general model of PB equations, called the Poisson–Boltzmann equations with
steric effects (PB-steric equations), which not only include (1.1)–(1.7) but
also have oscillatory total ionic charge density
$\sum\limits_{i=1}^{N}z_{i}c_{i}(\phi)$ under different assumptions of steric
effects and chemical potentials.

To get the PB-steric equations, we consider to use the Lennard-Jones (LJ)
potential which is a well-known model for the interaction between a pair of
ions (and solvent molecules) and is important in the density functional theory
and molecular dynamics simulation (cf. [2, 19, 25, 26, 37]). We set the energy
functional
$\iint_{\mathbb{R}^{d}}\psi(x-y)c_{i}(x)c_{j}(y)\,\mathrm{d}x\,\mathrm{d}y$ to
describe the energy of steric repulsion of $c_{i}$ and $c_{j}$ for
$i,j=0,1,\dots,N$. Here $\psi(z)=|z|^{-12}$ comes from the repulsive part of
LJ potential with strong singularity at the origin which makes the energy
functional
$\iint_{\mathbb{R}^{d}}\\!\psi(x-y)c_{i}(x)c_{j}(y)\,\mathrm{d}x\,\mathrm{d}y$
hard to study. Hence we replace the LJ potential by the approximate LJ
potentials (cf. [16, 30]) to describe the steric repulsion of ions and solvent
molecules and we obtain the following PB-steric equations:

$\displaystyle-\nabla\cdot(\varepsilon\nabla\phi)=4\pi\rho_{0}+4\pi\sum_{i=1}^{N}z_{i}e_{0}c_{i},$
(1.10) $\displaystyle
k_{B}T\ln(c_{i})+z_{i}e_{0}\phi+\Lambda\sum_{j=0}^{N}g_{ij}c_{j}=\Lambda\tilde{\mu}_{i}+\hat{\mu}_{i}\quad\text{for}~{}i=0,1,\dots,N,$
(1.11)

where $c_{0}$ is the concentration of solvent molecules, $c_{i}$ is the
concentration of $i$th ion species, and $\Lambda g_{ij}\geq 0$ represents the
strength of steric repulsion between the $i$th and $j$th ions (or solvent
molecules). One may see Appendix A for the derivation of (1.10)–(1.11) when
matrix $(g_{ij})$ is symmetric.

In this paper, we generalize (1.10)–(1.11) to all matrices $(g_{ij})$ (which
may include nonsymmetric matrices) such that (1.11) has a unique solution
$c_{i}=c_{i}(\phi)$ for $\phi\in\mathbb{R}$ and $i=0,1,\dots,N$. The
nonsymmetry of matrix $(g_{ij})$ may come from the potential
$U_{i}^{\text{wall}}$ which describes how the ions and solvent molecules
interact with the wall. Note that (1.10)–(1.11) can be expressed by
(1.8)–(1.9) with the chemical potential
$\mu_{i}=\Lambda\tilde{\mu}_{i}+\hat{\mu}_{i}$ and the excess chemical
potential $\mu_{i}^{\text{ex}}$ and the potential $U_{i}^{\text{wall}}$
satisfying
$\mu_{i}^{\text{ex}}+U_{i}^{\text{wall}}=\Lambda\sum\limits_{j=0}^{N}g_{ij}c_{j}$
for $i=0,1,\dots,N$. For the sake of simplify, we set $k_{B}T=4\pi e_{0}=1$
and write $\rho_{0}$ instead of $4\pi\rho_{0}$ so (1.10)–(1.11) become

$\displaystyle-\nabla\cdot\left(\varepsilon\nabla\phi\right)=\rho_{0}+\sum_{i=1}^{N}z_{i}c_{i}\quad\text{in}~{}\Omega,$
(1.12)
$\displaystyle\ln(c_{i})+z_{i}\phi+\Lambda\sum_{j=0}^{N}g_{ij}c_{j}=\Lambda\tilde{\mu}_{i}+\hat{\mu}_{i}\quad\text{for}~{}i=0,1,\dots,N,$
(1.13)

where $\Omega$ is a bounded smooth domain in $\mathbb{R}^{d}$ ($d\geq 2$).
Under suitable conditions of $g_{ij}$ such that system (1.13) has a unique
solution $c_{i}=c_{i}(\phi)$ (which may depend on parameter $\Lambda$) for
$\phi\in\mathbb{R}$ and $i=0,1,\dots,N$, (1.12) becomes a nonlinear elliptic
equation

$-\nabla\cdot(\varepsilon\nabla\phi)=\rho_{0}+\sum_{i=1}^{N}z_{i}c_{i}(\phi)\quad\text{in}~{}\Omega.$
(1.14)

Then for all possible $g_{ij}\geq 0$, equation (1.14) and its limiting
equation (as $\Lambda\to\infty$) are called as the PB-steric equations.
Hereafter, the boundary condition of (1.14) is considered as the Robin
boundary condition

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

where $\phi_{bd}\in\mathcal{C}^{2}(\partial\Omega)$ is the extra electrostatic
potential and $\eta\geq 0$ is a constant related to the surface dielectric
constant (cf. [3, 35, 39]).

The strength of steric effect is determined by $\Lambda g_{ij}$’s in the PB-
steric equations (1.13)–(1.14). As $\Lambda=0$ (or $g_{ij}=0$), (1.13) has the
same form as (1.2) and the PB-steric equations (1.13)–(1.14) become the
conventional PB equations (1.1)–(1.2). On the other hand, when $g_{ij}$’s are
positive constants, a larger $\Lambda$ produces stronger steric repulsion.
This motivates us to expect that as $\Lambda$ tends to infinity, the steric
repulsion becomes extremely strong so that volume exclusion holds true and
ions can be described by the lattice gas model as for the mean-field
approximation of the modified PB equations (1.3)–(1.4) (cf. [6, 14]). To
justify this, we prove that as $\Lambda\to\infty$, the PB-steric equations
(1.13)–(1.14) may approach to the modified PB equations (1.3)–(1.4) (see
Theorem 1.1 and Remark 1.2). Moreover, we prove that as $\Lambda\to\infty$,
the PB-steric equations (1.13)–(1.14) may approach to the modified PB
equations (1.5)–(1.7) (see Theorem 1.6 and Remark 1.7).

In order to obtain equations (1.3)–(1.4), we need the assumptions of steric
effects ($g_{ij}$) and chemical potentials ($\tilde{\mu}_{i}$) given by

  * (A1)

$g_{i0}=g_{00}=1-\gamma$ and $g_{ij}=g_{0j}=\gamma/Z$ for $i,j=1,\dots,N$,

  * (A2)

$\tilde{\mu}_{i}=\tilde{\mu}_{0}$ for $i=1,\dots,N$,

where $0\leq\gamma\leq 1$, $\tilde{\mu}_{0}>0$, $\hat{\mu}_{0}$ (in (1.13))
and $Z=\sum\limits_{i=1}^{N}|z_{i}|$ are constants independent of $\Lambda$.
Then by (A1) and (A2), system (1.13) has a unique solution
$c_{i,\Lambda}=c_{i,\Lambda}(\phi)$ for $\phi\in\mathbb{R}$ and
$i=0,1,\dots,N$, which satisfies

$\displaystyle
c_{i,\Lambda}=c_{0,\Lambda}\exp(\bar{\mu}_{i}-z_{i}\phi)\quad\text{for}~{}\phi\in\mathbb{R}~{}\text{and}~{}i=1,\dots,N,$
(1.16) $\displaystyle\ln(c_{0,\Lambda})+\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)c_{0,\Lambda}=\Lambda\tilde{\mu}_{0}+\hat{\mu}_{0}\quad\text{for}~{}\phi\in\mathbb{R},$
(1.17) $\displaystyle
H(\gamma,z_{1},\dots,z_{N},\phi)=1-\gamma+\frac{\gamma}{Z}\sum_{j=1}^{N}e^{\bar{\mu}_{j}-z_{j}\phi},$
(1.18)

where $\bar{\mu}_{i}=\hat{\mu}_{i}-\hat{\mu}_{0}$ is a constant independent of
$\Lambda$ for $i=1,\dots,N$. Note that (1.17) can be solved in terms of the
principal branch of Lambert W function (cf. [34]), which implies that
$c_{i,\Lambda}$ are positive smooth functions for $i=0,1,\dots,N$ (see
Proposition 2.1). Recall that the Lambert W function $W_{0}(x)$ can be defined
by $W_{0}(x)e^{W_{0}(x)}=x$ for all $x\geq e^{-1}$, and the range of
$W_{0}(x)$ is $[-1,\infty)$. Hence the PB-steric equation (1.14) can be
expressed as

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

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.20)

By the standard method of a nonlinear elliptic equation, we may obtain the
existence and uniqueness of (1.19) with the Robin boundary condition (1.15).

Using the implicit function theorem on Banach spaces (cf. [7, Theorem 15.1]),
we prove that $c_{i,\Lambda}$ converges to $c_{i}^{*}$ in space
$\mathcal{C}^{m}[a,b]$ as $\Lambda$ tends to infinity for $m\in\mathbb{N}$ and
$a<b$, where $c_{i}^{*}$ satisfies

$c_{i}^{*}(\phi)=\frac{\tilde{\mu}_{0}e^{\bar{\mu}_{i}-z_{i}\phi}}{H(\gamma,z_{1},\dots,z_{N},\phi)}\quad\text{for}~{}\phi\in\mathbb{R}~{}\text{and}~{}i=0,1,\dots,N.$
(1.21)

Thus function $f_{\Lambda}$ also converges to $f^{*}$ in space
$\mathcal{C}^{m}[a,b]$ 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)\quad\text{for}~{}\phi\in\mathbb{R}$
(1.22)

(see Corollary 2.5). Moreover, for the asymptotic limit of (1.19), we obtain
the following theorem.

######  Theorem 1.1.

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded smooth domain,
$\varepsilon\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 solution $\phi_{\Lambda}$ of PB-steric equation (1.19) with the Robin
boundary condition (1.15) satisfies

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

where $\phi^{*}$ is the solution of

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

with the Robin boundary condition (1.15).

######  Remark 1.2.

When $N=2$, (1.21)–(1.23) become (1.3)–(1.4). Besides, (1.23) is the limiting
equation of (1.19) as $\Lambda\to\infty$. This shows the PB-steric equations
include the modified PB equations (1.3)–(1.4).

######  Remark 1.3.

When $N=2$, we replace the assumption (A1) by

  * (A1)

′ $g_{i0}=g_{00}=1-\gamma$, $g_{ij}=g_{0j}=\gamma_{j}$ for $i,j=0,1,2$, and
$\gamma=\gamma_{1}+\gamma_{2}$.

Then by (A1) ′ and (A2), system (1.13) has a unique solution
$c_{i,\Lambda}=c_{i,\Lambda}(\phi)$ for $\phi\in\mathbb{R}$ and $i=0,1,2$. Due
to $z_{1}z_{2}<0$, we can calculate directly to get
$\mathrm{d}f_{\Lambda}/\mathrm{d}\phi<0$, which implies function $f_{\Lambda}$
is decreasing to $\phi$. Hence one can follow the similar argument in sections
2 and 3 to obtain the same conclusion of Theorem 1.1.

######  Remark 1.4.

When $N\geq 2$, we replace the assumption (A1) by

  * (A1)

′′ $g_{i0}=g_{00}=1-\gamma$, and $g_{ij}=g_{0j}=\gamma|z_{j}|/Z$ for
$i,j=1,\dots,N$. Then by (A1)′′ and (A2), system (1.13) has a unique solution
$c_{i,\Lambda}=c_{i,\Lambda}(\phi)$ for $\phi\in\mathbb{R}$ and
$i=0,1,\dots,N$. Since $z_{i}z_{j}<0$ for some $i,j\in\\{1,\dots,N\\}$, we can
calculate directly to get $\mathrm{d}f_{\Lambda}/\mathrm{d}\phi<0$, which
implies function $f_{\Lambda}$ is decreasing to $\phi$. Hence one can follow
the similar argument in sections 2 and 3 to obtain the same conclusion of
Theorem 1.1.

######  Remark 1.5.

As $N=3$, we replace the assumption (A1) by

  * (A1)

′′′ $g_{i0}=g_{00}=1-\gamma$, $g_{ij}=g_{0j}=\gamma_{j}$ for $i,j=0,1,2,3$
and $\gamma=\gamma_{1}+\gamma_{2}+\gamma_{3}$.

Then by (A1) ′′′ and (A2), system (1.13) has a unique solution
$c_{i,\Lambda}=c_{i,\Lambda}(\phi)$ for $\phi\in\mathbb{R}$ and $i=0,1,2,3$.
But for some $\Lambda$, function $f_{\Lambda}$ and
$f_{\Lambda}\circ\phi_{\Lambda}$ (total ionic charge density) may become
oscillatory (see Figure 1). One can see the proof in Appendix B and numerical
methods in section 4. Such oscillatory total ionic charge density
$f_{\Lambda}=f_{\Lambda}(\phi)$ cannot be obtained in the conventional and
modified PB equations (1.1)–(1.7).

Figure 1: The profiles of $f_{\Lambda}(\phi)$ and
$f_{\Lambda}(\phi_{\Lambda}(x))$ under assumptions (A1)′′′ and (A2). In (a),
curves 1–4 are profiles of function
$f_{\Lambda}=\sum_{i=1}^{3}z_{i}c_{i,\Lambda}$ with $\Lambda=0.5$, $1$, $2$,
$4$. In (b), curves 1–4 are the profiles of function
$f_{\Lambda}\circ\phi_{\Lambda}$ with $\Lambda=0.5$, $1$, $2$, $4$.

To obtain (1.5)–(1.7), we need the assumptions of $g_{ij}$ and
$\tilde{\mu}_{i}$ given by

  * (A3)

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

  * (A4)

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

where $\lambda_{i}>0$, $\tilde{\mu}_{0}>0$ and $\hat{\mu}_{i}$ (in (1.13)) are
constants independent of $\Lambda$. By (A3) and (A4), system (1.13) has a
unique solution $c_{i,\Lambda}=c_{i,\Lambda}(\phi)$ for $\phi\in\mathbb{R}$
and $i=0,1,\dots,N$, which satisfies

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

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$. Note that as
$\lambda_{j}/\lambda_{0}=1$ for $j=1,\dots,N$, (1.25) can be solved by the
Lambert W function, but here some $\lambda_{j}/\lambda_{0}$ may not be equal
to one so we may call $c_{0,\Lambda}$ as a Lambert-type function. Then we
apply the implicit function theorem on (1.25) and obtain that
$c_{0,\Lambda}(\phi)$ is a positive smooth function (see Proposition 2.7).
Hence the PB-steric equation (1.14) can be expressed as

$-\nabla\cdot(\varepsilon\nabla\phi_{\Lambda})=\rho_{0}+\tilde{f}_{\Lambda}(\phi_{\Lambda})\quad\text{in}~{}\Omega,$
(1.26)

where function $\tilde{f}_{\Lambda}=\tilde{f}_{\Lambda}(\phi)$ is denoted as

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

Note that the existence and uniqueness of (1.26) with the Robin boundary
condition (1.15) can be obtain by the standard method of the nonlinear
elliptic equation, but (1.19) and (1.26) are different because they have
different nonlinear terms $f_{\Lambda}$ and $\tilde{f}_{\Lambda}$.

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

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

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

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

(see Corollary 2.9). Moreover, for the asymptotic limit of (1.26), we have the
following theorem.

######  Theorem 1.6.

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded smooth domain,
$\varepsilon\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 (A3) and (A4) hold true. Then
the solution $\phi_{\Lambda}$ of PB-steric equation (1.26) with the Robin
boundary condition (1.15) satisfies

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

where $\phi^{*}$ is the solution of

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

with the Robin boundary condition (1.15).

######  Remark 1.7.

As $\lambda_{i}=v_{i}/v_{0}$, $\tilde{\mu}_{0}=1/v_{0}$, $c_{i}^{*}=c_{i}$ and
$\phi^{*}=\phi$, (1.28)–(1.31) become (1.5)–(1.7) (up to scalar multiples).
Besides, (1.31) is the limiting equation of (1.26) as $\Lambda\to\infty$. This
shows the PB-steric equations include the modified PB equations (1.5)–(1.7).

The rest of the paper is organized as follows. In section 2, we analyze of
functions $f_{\Lambda}$, $f^{*}$, $\tilde{f}_{\Lambda}$ and $\tilde{f}^{*}$
under the assumptions (A1)–(A2) and (A3)–(A4), respectively. The proof of
Theorems 1.1 and 1.6 are stated in section 3. Numerical schemes are shown in
section 4.

##  2 Analysis of nonlinear terms under (A1)–(A2) and (A3)–(A4)

###  2.1 Analysis of $c_{i,\Lambda}$ and $f_{\Lambda}$ under (A1) and (A2)

In this section, we firstly use assumptions (A1) and (A2) and apply Gaussian
elimination to solve system (1.13) and obtain a unique positive smooth
solution $c_{i,\Lambda}=c_{i,\Lambda}(\phi)$ for $\phi\in\mathbb{R}$ and
$i=0,1,\dots,N$. Then we establish the asymptotic behavior of $c_{i,\lambda}$
(as $\phi$ tends to infinity) and the strict decrease and unboundedness of
$f_{\Lambda}=\sum\limits_{i=1}^{N}z_{i}c_{i,\Lambda}$ (see Proposition 2.1).

For functions $c_{i,\Lambda}$ and $f_{\Lambda}$ under the assumptions (A1) and
(A2), we have the following.

######  Proposition 2.1.

Assume that $z_{0}=0$, (A1) and (A2) hold true. Suppose $z_{i}z_{j}<0$ for
some $i,j\in\\{1,\dots,N\\}$.

  * (i)

System (1.13) can be solved uniquely by

$c_{i,\Lambda}(\phi)=\frac{W_{0}\left(\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)e^{\Lambda\tilde{\mu}_{0}+\hat{\mu}_{0}}\right)}{\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)}e^{\bar{\mu}_{i}-z_{i}\phi}$ (2.1)

for $\phi\in\mathbb{R}$ and $i=0,1,\dots,N$, where $W_{0}$ denotes the
principal branch of Lambert W function (cf. [34]).

  * (ii)

Function $f_{\Lambda}=\sum\limits_{i=1}^{N}z_{i}c_{i,\Lambda}$ is strictly
decreasing on $\mathbb{R}$.

  * (iii)

Let $I=\\{i:z_{i}=\max\limits_{0\leq k\leq N}z_{k}\\}$ and
$J=\\{j:z_{j}=\min\limits_{0\leq k\leq N}z_{k}\\}$. Then

    * (a)

$\lim\limits_{\phi\to-\infty}c_{i_{0},\Lambda}(\phi)=\infty$ for $i_{0}\in I$,
and $\lim\limits_{\phi\to-\infty}c_{k,\Lambda}(\phi)=0$ for $k\notin I$;

    * (b)

$\lim\limits_{\phi\to\infty}c_{j_{0},\Lambda}(\phi)=\infty$ for $j_{0}\in J$,
and $\lim\limits_{\phi\to\infty}c_{k,\Lambda}(\phi)=0$ for $k\notin J$.

  * (iv)

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

###### Proof.

By (A1) and (A2), system (1.13) can be expressed as

$\ln(c_{i,\Lambda})+z_{i}\phi+\Lambda\left((1-\gamma)c_{0,\Lambda}+\frac{\gamma}{Z}\sum_{j=1}^{N}c_{j,\Lambda}\right)=\Lambda\tilde{\mu}_{0}+\hat{\mu}_{i}\quad\text{for}~{}i=0,1,\dots,N.$
(2.2)

Subtracting (2.2) for $i\neq 0$ from that for $i=0$, we get
$\ln(c_{i,\Lambda}/c_{0,\Lambda})+z_{i}\phi=\hat{\mu}_{i}-\hat{\mu}_{0}:=\bar{\mu}_{i}$,
which implies (1.16) under the assumption $z_{0}=0$. Then we plug (1.16) into
(2.2) for $i=0$ to get (1.17), which implies

$\left[\Lambda H(\gamma,z_{1},\dots,z_{N},\phi)c_{0,\Lambda}\right]e^{\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)c_{0,\Lambda}}=\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)e^{\Lambda\tilde{\mu}_{0}+\hat{\mu}_{0}}$
(2.3)

for $\phi\in\mathbb{R}$. From theorems of Lambert W function in [34], (2.3)
has a unique smooth positive solution $c_{0,\Lambda}$ denoted by

$c_{0,\Lambda}(\phi)=\frac{W_{0}\left(\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)e^{\Lambda\tilde{\mu}_{0}+\hat{\mu}_{0}}\right)}{\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)}\quad\text{for}~{}\phi\in\mathbb{R}.$ (2.4)

By (1.16) and (2.4), we arrive at (2.1) and complete the proof of (i).

Next we state the proof of (ii). We differentiate (1.17) with respect to
$\phi$ and get

$\frac{1}{c_{0,\Lambda}}\frac{\mathrm{d}c_{0,\Lambda}}{\mathrm{d}\phi}+\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)\frac{\mathrm{d}c_{0,\Lambda}}{\mathrm{d}\phi}-\Lambda
c_{0,\Lambda}\frac{\gamma}{Z}\sum_{j=1}^{N}z_{j}e^{\bar{\mu}_{j}-z_{j}\phi}=0,$

which gives

$\frac{\mathrm{d}c_{0,\Lambda}}{\mathrm{d}\phi}=\frac{\displaystyle\frac{\gamma}{Z}\sum_{j=1}^{N}z_{j}e^{\bar{\mu}_{j}-z_{j}\phi}}{(\Lambda
c_{0,\Lambda})^{-1}+H(\gamma,z_{1},\dots,z_{N},\phi)}c_{0,\Lambda}.$ (2.5)

On the other hand, by (1.16), we write $f_{\Lambda}$ as
$f_{\Lambda}(\phi)=\sum\limits_{i=1}^{N}z_{i}c_{0,\Lambda}e^{\bar{\mu}_{i}-z_{i}\phi}$
for $\phi\in\mathbb{R}$. Then differentiating $f_{\Lambda}$ with respect to
$\phi$ and using (2.5), we get

$\frac{\mathrm{d}f_{\Lambda}}{\mathrm{d}\phi}=\left(\frac{\displaystyle\frac{\gamma}{Z}\left(\sum_{i=1}^{N}z_{i}e^{\bar{\mu}_{i}-z_{i}\phi}\right)\left(\sum_{i=1}^{N}z_{i}e^{\bar{\mu}_{i}-z_{i}\phi}\right)}{(\Lambda
c_{0,\Lambda})^{-1}+H(\gamma,z_{1},\dots,z_{N},\phi)}-\sum_{i=1}^{N}z_{i}^{2}e^{\bar{\mu}_{i}-z_{i}\phi})\right)c_{0,\Lambda}.$
(2.6)

Note the inequality

$\displaystyle~{}~{}~{}\left(\sum_{i=1}^{N}z_{i}^{2}e^{\bar{\mu}_{i}-z_{i}\phi}\right)\left[(\Lambda
c_{0,\Lambda})^{-1}+H(\gamma,z_{1},\dots,z_{N},\phi)\right]>\frac{\gamma}{Z}\left(\sum_{i=1}^{N}z_{i}^{2}e^{\bar{\mu}_{i}-z_{i}\phi}\right)\left(\sum_{i=1}^{N}e^{\bar{\mu}_{i}-z_{i}\phi}\right)>\frac{\gamma}{Z}\left(\sum_{i=1}^{N}z_{i}e^{\bar{\mu}_{i}-z_{i}\phi}\right)^{2},$
(2.7)

where we have used the fact that $z_{i}z_{j}<0$ for some
$i,j\in\\{1,\dots,N\\}$. Therefore, by (2.6) and (2.7), we have
$\mathrm{d}f_{\Lambda}/\mathrm{d}\phi<0$ on $\mathbb{R}$ and the proof of (ii)
is complete.

To prove (iii), we need the following claim.

######  Claim 2.2.

There holds that

$\lim_{\phi\to\pm\infty}\frac{W_{0}(\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)e^{\Lambda\tilde{\mu}_{0}+\hat{\mu}_{0}})}{\ln\Lambda+\ln(H(\gamma,z_{1},\dots,z_{N},\phi))+\Lambda\tilde{\mu}_{0}+\hat{\mu}_{0}}=1.$
(2.8)

###### Proof of Claim 2.2.

Recall the asymptotic behavior of Lambert W function in [15, 34]

$\ln x-\ln\ln x+\frac{\ln\ln x}{2\ln x}\leq W_{0}(x)\leq\ln x-\ln\ln
x+\frac{e}{e-1}\frac{\ln\ln x}{\ln x}\quad\text{for}~{}x\geq e.$ (2.9)

Since $z_{i}z_{j}<0$ for some $i,j\in\\{1,\dots,N\\}$, it is clear that
$\lim\limits_{\phi\to\pm\infty}\sum\limits_{i=1}^{N}e^{\bar{\mu}_{i}-z_{i}\phi}=\infty$,
which implies that

$\lim\limits_{\phi\to\pm\infty}\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)e^{\Lambda\tilde{\mu}_{0}+\hat{\mu}_{0}}=\infty.$

Then by (2.9), we obtain (2.8) and complete the proof of Claim 2.2. ∎

Now we state the proof of (iii). Since $z_{i}z_{j}<0$ for some
$i,j\in\\{1,\dots,N\\}$, we know $z_{i_{0}}>0$ for $i_{0}\in I$ and
$z_{j_{0}}<0$ for $j_{0}\in J$. Then for $i_{0}\in I$, we may use (1.16),
(2.4) and Claim 2.2 to get

$\displaystyle\lim_{\phi\to-\infty}c_{i_{0},\Lambda}(\phi)$
$\displaystyle=\lim_{\phi\to-\infty}\frac{e^{\bar{\mu}_{i_{0}}-z_{i_{0}}\phi}}{\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)}W_{0}(\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)e^{\Lambda\tilde{\mu}_{0}+\hat{\mu}_{0}})$
$\displaystyle=\frac{Z}{\Lambda\gamma\sum\limits_{i\in
I}e^{\bar{\mu}_{i}-\bar{\mu}_{i_{0}}}}\lim_{\phi\to-\infty}W_{0}(\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)e^{\Lambda\tilde{\mu}_{0}+\hat{\mu}_{0}})=\infty.$

On the other hand, for $k\notin I$, we have

$\lim_{\phi\to-\infty}c_{k,\Lambda}(\phi)=\lim_{\phi\to-\infty}\frac{e^{\bar{\mu}_{k}-z_{k}\phi}[\left(\ln\Lambda+\ln(H(\gamma,z_{1},\dots,z_{N},\phi))+\Lambda\tilde{\mu}_{0}+\hat{\mu}_{0}\right)]}{\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)}=0.$

Thus the proof of (iii)(a) is complete. The proof of (iii)(b) is similar to
(iii)(a) by using (1.16), (1.17), and Claim 2.2. Therefore, we complete the
proof of (iii).

Finally, we give the proof of (iv). By (1.21), function $f_{\Lambda}$ can be
denoted as $\displaystyle f_{\Lambda}(\phi)=\sum_{i\in
I}z_{i}c_{i,\Lambda}(\phi)+\sum_{i\notin I}z_{i}c_{i,\Lambda}(\phi)$ for
$\phi\in\mathbb{R}$, where $I$ is defined in (iii). Then by (iii), we have
$\displaystyle\lim_{\phi\to-\infty}f_{\Lambda}(\phi)=\infty$ and
$\displaystyle\lim_{\phi\to\infty}f_{\Lambda}(\phi)=-\infty$. Therefore, we
conclude (iv) and complete the proof of Proposition 2.1. ∎

###  2.2 Analysis of $c_{i}^{*}$ and $f^{*}$ under (A1) and (A2)

Function $c_{0}^{*}$ is the limit
$\lim\limits_{\Lambda\to\infty}c_{0,\Lambda}$, where function $c_{0,\Lambda}$
is the solution of (1.17) for $\Lambda>0$. Let $\delta=\Lambda^{-1}$ and
$\tilde{c}_{0,\delta}=c_{0,\Lambda}$. Then by (1.17), $\tilde{c}_{0,\delta}$
satisfies

$\delta\ln(\tilde{c}_{0,\delta}(\phi))+H(\gamma,z_{1},\dots,z_{N},\phi)\tilde{c}_{0,\delta}=\tilde{\mu}_{0}+\delta\hat{\mu}_{0}\quad\text{for}~{}\phi\in\mathbb{R}.$
(2.10)

Notice that $\Lambda\to\infty$ is equivalent to $\delta\to 0^{+}$ so
$c_{0}^{*}$ also equals the limit $\lim\limits_{\delta\to
0^{+}}\tilde{c}_{0,\delta}$. Moreover $c_{0}^{*}$ is defined by (1.21), which
is (2.10) with $\delta=0$. 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 Propositions 2.3, so by (1.16), we obtain the
convergence of $c_{i,\Lambda}$ as $\Lambda\to\infty$.

######  Proposition 2.3.

Let $c_{i}^{*}$ be defined in (1.21).

  * (i)

For $\phi\in\mathbb{R}$,
$\lim\limits_{\Lambda\to\infty}c_{i,\Lambda}(\phi)=c_{i}^{*}(\phi)$ for
$i=0,1,\dots,N$.

  * (ii)

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

###### Proof.

To prove (i), we need the following claim.

######  Claim 2.4.

Assume that (A1) and (A2) hold ture. Then there holds that

$\lim_{\Lambda\to\infty}\frac{W_{0}(\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)e^{\Lambda\tilde{\mu}_{0}+\hat{\mu}_{0}})}{\Lambda\tilde{\mu}_{0}}=1\quad\text{for}~{}\phi\in\mathbb{R}.$

###### Proof of Claim 2.4.

The proof is similar to Claim 2.2 so we omit it here.∎

By Proposition 2.1(i) and Claim 2.4, (i) follows from

$\lim_{\Lambda\to\infty}c_{i,\Lambda}(\phi)=e^{\bar{\mu}_{i}-z_{i}\phi}\lim_{\Lambda\to\infty}\frac{\Lambda\tilde{\mu}_{0}}{\Lambda
H(\gamma,z_{1},\dots,z_{N},\phi)}=\frac{\tilde{\mu}_{0}e^{\bar{\mu}_{i}-z_{i}\phi}}{H(\gamma,z_{1},\dots,z_{N},\phi)}=c_{i}^{*}(\phi).$

To prove (ii), we 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 notational
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 (1.16), it suffices to show that
$\lim\limits_{\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}$, (2.10) can be denoted as
$K_{1}(w_{\delta}(\phi),\delta)=0$ for $\delta>0$ and $\phi\in[a,b]$, where
$K_{1}$ is a $\mathcal{C}^{1}$-function on
$\mathcal{C}^{m}[a,b]\times\mathbb{R}$ defined by

$K_{1}(w(\phi),\delta)=\delta
w(\phi)+H(\gamma,z_{1},\dots,z_{N},\phi)e^{w(\phi)}-\tilde{\mu}_{0}-\delta\hat{\mu}_{0}$
(2.11)

for all $w\in\mathcal{C}^{m}[a,b]$ and $\phi\in[a,b]$. Note that
$K_{1}(w^{*},0)=0$ by (i). A direct calculation for Fréchet derivative of
(2.11) gives
$D_{w}K_{1}(w(\phi),\delta)=\delta+H(\gamma,z_{1},\dots,z_{N},\phi)e^{w(\phi)}$
for all $w\in\mathcal{C}^{m}[a,b]$ and $\phi\in[a,b]$. Then by (i), we get
$D_{w}K_{1}(w^{*}(\phi),0)=\tilde{\mu}_{0}>0$ for all $\phi\in[a,b]$. This
implies that $D_{w}K_{1}(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. [7, Corollary 15.1]),
there exist an open subset
$B_{\delta_{0}}(w^{*})\times(-\delta_{0},\delta_{0})\subset\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^{*})\subset\mathcal{C}^{m}[a,b]$ for
$\delta\in(-\delta_{0},\delta_{0})$ such that
$K_{1}(\tilde{w}(\cdot,\delta),\delta)=0$ 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(i), equation $K_{1}(w,\delta)=0$ has a unique solution $w=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^{+}}\|w_{\delta}-w^{*}\|_{\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.3(ii). ∎

######  Corollary 2.5.

$\displaystyle\lim_{\Lambda\to\infty}\|f_{\Lambda}-f^{*}\|_{\mathcal{C}^{m}[a,b]}=0$
for $m\in\mathbb{N}$ and $a<b$.

###### Proof.

It follows from Proposition 2.3(ii) and the fact that $\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}$. ∎

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

######  Proposition 2.6.

Let $f^{*}$ be the function defined in (1.22).

  * (i)

Function $f^{*}$ is strictly decreasing on $\mathbb{R}$.

  * (ii)

Function $f^{*}$ satisfies $m^{*}<f^{*}(\phi)<M^{*}$ for all
$\phi\in\mathbb{R}$, where

$m^{*}=\lim_{\phi\to\infty}f^{*}(\phi)<0~{}~{}\text{and}~{}~{}M^{*}=\lim_{\phi\to-\infty}f^{*}(\phi)>0.$

###### Proof.

By (1.22), $f^{*}$ can be expressed as

$f^{*}(\phi)=\frac{\tilde{\mu}_{0}}{H(\gamma,z_{1},\dots,z_{N},\phi)}\sum_{i=1}^{N}z_{i}e^{\bar{\mu}_{i}-z_{i}\phi}\quad\text{for}~{}\phi\in\mathbb{R}.$
(2.12)

Then differentiating (2.12) with respect to $\phi$ gives

$\displaystyle\frac{\mathrm{d}f^{*}}{\mathrm{d}\phi}$
$\displaystyle=\tilde{\mu}_{0}\left(\frac{\displaystyle\frac{\gamma}{Z}\left(\sum_{i=1}^{N}z_{i}e^{\bar{\mu}_{i}-z_{i}\phi}\right)^{2}}{\left(H(\gamma,z_{1},\dots,z_{N},\phi)\right)^{2}}-\frac{\sum_{i=1}^{N}z_{i}^{2}e^{\bar{\mu}_{i}-z_{i}\phi}}{H(\gamma,z_{1},\dots,z_{N},\phi)}\right)\leq\frac{\tilde{\mu}_{0}\gamma}{Z}\frac{\displaystyle\left(\sum_{i=1}^{N}z_{i}e^{\bar{\mu}_{i}-z_{i}\phi}\right)^{2}-\sum_{i=1}^{N}e^{\bar{\mu}_{i}-z_{i}\phi}\sum_{i=1}^{N}z_{i}^{2}e^{\bar{\mu}_{i}-z_{i}\phi}}{\left(H(\gamma,z_{1},\dots,z_{N},\phi)\right)^{2}}<0$

for $\phi\in\mathbb{R}$. Here the last inequality comes from the fact that
$z_{i}z_{j}<0$ for some $i,j\in\\{1,\dots,N\\}$. Therefore, we complete the
proof of (i).

To prove (ii), we note that

$m^{*}:=\tilde{\mu}_{0}\lim_{\phi\to\infty}\frac{\displaystyle\left(\sum_{j\in
J}+\sum_{j\notin J}\right)z_{j}e^{\bar{\mu}_{j}-z_{j}\phi}}{\displaystyle
1-\gamma+\frac{\gamma}{Z}\left(\sum_{j\in J}+\sum_{j\notin
J}\right)e^{\bar{\mu}_{j}-z_{j}\phi}}=\frac{\tilde{\mu}_{0}Z}{\gamma}\frac{\displaystyle\sum_{j\in
J}z_{j}e^{\bar{\mu}_{j}}}{\displaystyle\sum_{j\in J}e^{\bar{\mu}_{j}}}<0$

and

$M^{*}:=\tilde{\mu}_{0}\lim_{\phi\to-\infty}\frac{\displaystyle\left(\sum_{i\in
I}+\sum_{i\notin I}\right)z_{i}e^{\bar{\mu}_{i}-z_{i}\phi}}{\displaystyle
1-\gamma+\frac{\gamma}{Z}\left(\sum_{i\in I}+\sum_{i\notin
I}\right)e^{\bar{\mu}_{i}-z_{i}\phi}}=\frac{\tilde{\mu}_{0}Z}{\gamma}\frac{\displaystyle\sum_{i\in
I}z_{i}e^{\bar{\mu}_{i}}}{\displaystyle\sum_{i\in I}e^{\bar{\mu}_{i}}}>0,$

where $I$ and $J$ are defined in Proposition 2.1(iii). By (i), $f^{*}$ is
strictly decreasing so we have $m^{*}<f^{*}(\phi)<M^{*}$ for all
$\phi\in\mathbb{R}$ and complete the proof of Proposition 2.6. ∎

###  2.3 Analysis of $c_{i,\Lambda}$ and $\tilde{f}_{\Lambda}$ under (A3) and
(A4)

In this section, we use assumptions (A3) and (A4) and apply Gaussian
elimination to solve system (1.13) and obtain a unique positive smooth
solution $c_{i,\Lambda}=c_{i,\Lambda}(\phi)$ for $\phi\in\mathbb{R}$ and
$i=0,1,\dots,N$. Then we establish the asymptotic behavior of $c_{i,\Lambda}$
(as $\phi$ tends to infinity), and the strict decrease and unboundedness of
$\displaystyle\tilde{f}_{\Lambda}=\sum_{i=1}^{N}z_{i}c_{i,\Lambda}$ (see
Proposition 2.7).

For functions $c_{i,\Lambda}$ and $\tilde{f}_{\Lambda}$, we have the following
proposition.

######  Proposition 2.7.

Assume that $z_{0}=0$, (A3) and (A4) hold true. Suppose $z_{i}z_{j}<0$ for
some $i,j\in\\{1,\dots,N\\}$.

  * (i)

System (1.13) has a unique, smooth, and positive solution
$c_{i,\Lambda}=c_{i,\Lambda}(\phi)$ for $\phi\in\mathbb{R}$ and
$i=0,1,\dots,N$.

  * (ii)

Function $\tilde{f}_{\Lambda}=\sum\limits_{i=1}^{N}z_{i}c_{i,\Lambda}$ is
strictly decreasing on $\mathbb{R}$.

  * (iii)

    1. (a)

If $z_{k}\geq 0$ ($z_{k}\leq 0$), then $\sup\limits_{\phi\geq
0}c_{k,\Lambda}(\phi)\leq e^{\mu_{k}}$ ($\sup\limits_{\phi\leq
0}c_{k,\Lambda}(\phi)\leq e^{\mu_{k}}$).

    2. (b)

There exist $i_{0},j_{0}\in\\{1,\dots,N\\}$, $z_{i_{0}}>0$ and $z_{j_{0}}<0$
such that $\limsup\limits_{\phi\to-\infty}c_{i_{0},\Lambda}(\phi)=\infty$ and
$\limsup\limits_{\phi\to\infty}c_{j_{0},\Lambda}(\phi)=\infty$.

  * (iv)

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

###### Proof.

By (A3) and (A4), system (1.13) can be expressed as

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

Following a similar process in the proof of Proposition 2.1(i), we can use
(2.13) to obtain (1.24)–(1.25). Then we denote (1.25) as
$k_{1}(c_{0,\Lambda},\phi)=0$, where $k_{1}$ is defined by

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

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

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

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

To prove (ii), we differentiate (2.13) with respect to $\phi$ and obtain

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

where
$D_{\Lambda}=\text{diag}(\frac{1}{c_{0,\Lambda}},\cdots,\frac{1}{c_{N,\Lambda}})$
is positive definite,
$G=[\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
$D_{\Lambda}+\Lambda G$ is positive definite and invertible with inverse
matrix $(D_{\Lambda}+\Lambda G)^{-1}$ which is also positive definite. Then
(2.14) gives $\mathrm{d}{\bf c}_{\Lambda}/\mathrm{d}\phi=-(D_{\Lambda}+\Lambda
G)^{-1}{\bf z}$, and $\mathrm{d}\tilde{f}_{\Lambda}/\mathrm{d}\phi$ becomes

$\frac{\mathrm{d}\tilde{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}_{\Lambda}}{\mathrm{d}\phi}=-{\bf
z}^{\mathsf{T}}\left(D_{\Lambda}+\Lambda 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
(ii).

To prove (iii), we firstly suppose that $z_{k}\geq 0$ for some
$k\in\\{0,\dots,N\\}$. Then (2.13) implies $\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$. The
proof for the case of $z_{k}\leq 0$ is similar to the case of $z_{k}\geq 0$ so
we omit it here. Hence the proof of (iii)(a) is complete. On the other hand,
for (iii)(b), since $z_{i}z_{j}<0$ for some $i,j\in\\{1,\dots,N\\}$, then both
index sets $I^{\prime}=\\{i:z_{i}>0\\}$ and $J^{\prime}=\\{j:z_{j}<0\\}$ are
nonempty. Now we claim that there exists $i_{0}\in I^{\prime}$ such that
$\displaystyle\limsup_{\phi\to-\infty}c_{i_{0},\Lambda}(\phi)=\infty$. We
prove this by contradiction. Suppose $\sup_{\phi\leq
0}c_{i,\Lambda}(\phi)<\infty$ for all $i\in I^{\prime}$. Then there exists
$K_{1}>0$ such that $0<c_{i,\Lambda}(\phi)<K_{1}$ for $\phi\leq 0$ and $i\in
I^{\prime}$. By equation (2.13) and (iii)(a), we have

$z_{i}\phi=\mu_{i}-\ln(c_{i,\Lambda}(\phi))-\Lambda\lambda_{i}\sum_{k=0}^{N}\lambda_{j}c_{k,\Lambda}(\phi)\geq\mu_{i}-\ln(K_{1})-\Lambda\lambda_{i}\left(\sum_{k\in
I^{\prime}}\lambda_{k}K_{1}+\sum_{k\notin
I^{\prime}}\lambda_{k}e^{\mu_{k}}\right),$

for $i\in I^{\prime}$ and $\phi\leq 0$, which leads a contradiction by letting
$\phi\to-\infty$. Hence there exists $i_{0}\in I^{\prime}$ such that
$\displaystyle\limsup_{\phi\to-\infty}c_{i_{0},\Lambda}(\phi)=\infty$. The
proof for the case that there exists $j_{0}\in J^{\prime}$ such that
$\displaystyle\limsup_{\phi\to\infty}c_{j_{0},\Lambda}(\phi)=\infty$ is
similar to the case that
$\displaystyle\limsup_{\phi\to-\infty}c_{i_{0},\Lambda}(\phi)=\infty$ for some
$i_{0}\in I^{\prime}$ so we omit it here. Therefore, the proof of (iii)(b) is
complete.

The proof of (iv) is similar to the proof of Proposition 2.1(iv) so we omit it
here. The proof of Proposition 2.7 is complete. ∎

###  2.4 Analysis of $c_{i}^{*}$ and $\tilde{f}^{*}$ under (A3) and (A4)

In this section, we apply the implicit function theorem to show the existence
of $c_{i}^{*}$ and use the implicit function theorem on Banach space to show
the convergence of $c_{0,\Lambda}$ as $\Lambda\to\infty$ (see Proposition
2.8). Then by (1.24) and (1.27), we can obtain the convergence of
$c_{i,\Lambda}$ and $f_{\Lambda}$ as $\Lambda\to\infty$.

######  Proposition 2.8.

Let $c_{i}^{*}$ be defined in (1.28)–(1.29).

  * (i)

Equations (1.28)–(1.29) 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$.

  * (ii)

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

###### Proof.

The proof of the proposition is similar to that of Proposition 2.3 so we omit
it here. ∎

######  Corollary 2.9.

$\displaystyle\lim_{\Lambda\to\infty}\|\tilde{f}_{\Lambda}-\tilde{f}^{*}\|_{\mathcal{C}^{m}[a,b]}=0$
for $m\in\mathbb{N}$ and $a<b$.

###### Proof.

The proof of the corollary is similar to that of Corollary 2.5 so we omit it
here. ∎

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

######  Proposition 2.10.

Let $\tilde{f}^{*}$ be defined in (1.30).

  * (i)

Function $\tilde{f}^{*}$ is strictly decreasing on $\mathbb{R}$.

  * (ii)

Function $\tilde{f}^{*}$ satisfies $m^{*}<\tilde{f}^{*}(\phi)<M^{*}$ for all
$\phi\in\mathbb{R}$, where

$m^{*}=\lim_{\phi\to\infty}\tilde{f}^{*}(\phi)<0~{}\text{and}~{}M^{*}=\lim_{\phi\to-\infty}\tilde{f}^{*}(\phi)>0.$

###### Proof.

Differentiating (1.28)–(1.29) with respect to $\phi$ gives

$\sum_{j=0}^{N}\lambda_{j}\frac{\mathrm{d}c_{j}^{*}}{\mathrm{d}\phi}=0,\quad\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,$

which implies
$\displaystyle\frac{\mathrm{d}c_{0}^{*}}{\mathrm{d}\phi}=\lambda_{0}c_{0}^{*}\frac{\displaystyle\sum_{i=1}^{N}z_{i}\lambda_{i}c_{i}^{*}}{\displaystyle\sum_{i=0}^{N}\lambda_{i}^{2}c_{i}^{*}}$
for $\phi\in\mathbb{R}$. Consequently, we have

$\frac{\mathrm{d}\tilde{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{\displaystyle\left(\sum_{i=1}^{N}z_{i}\lambda_{i}c_{i}^{*}\right)^{2}}{\displaystyle\sum_{i=0}^{N}\lambda_{i}^{2}c_{i}^{*}}-\sum_{i=1}^{N}z_{i}^{2}c_{i}^{*}<0$

for $\phi\in\mathbb{R}$. Here the last inequality comes from the
Cauchy–Schwarz 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 (i).

To prove boundedness of $\tilde{f}^{*}$, we note that (1.29) implies
$0<c_{i}^{*}(\phi)<\tilde{\mu}_{0}/\lambda_{i}$ for $\phi\in\mathbb{R}$ and
$i=0,1,\dots,N$, and

$|\tilde{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.29) and (1.30), function $\tilde{f}^{*}$ can be
represented as

$\tilde{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 (i) and (2.15), function $\tilde{f}^{*}$ is strictly decreasing and bounded
on $\mathbb{R}$ so the limit
$\displaystyle\lim_{\phi\to\infty}\tilde{f}^{*}(\phi)$, denoted by $m^{*}$,
exists and is finite. Since $c_{i}^{*}<\tilde{\mu}_{0}/\lambda_{i}$ for
$i=0,1\dots,N$, (1.28) implies

$\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.28) 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}\tilde{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 (1.29) (by letting $\phi\to\infty$) and
$\tilde{\mu}_{0}>0$. Similarly, we may prove
$\displaystyle\lim_{\phi\to-\infty}\tilde{f}^{*}(\phi)=M^{*}>0$ and complete
the proof of Proposition 2.10. ∎

##  3 Proof of Theorems 1.1 and 1.6

Since $f_{\Lambda}$ and $\tilde{f}_{\Lambda}$ are unbounded on $\mathbb{R}$
but $f^{*}$ and $\tilde{f}^{*}$ are bounded on $\mathbb{R}$ so we cannot
obtain the uniform convergence of $f_{\Lambda}$ and $\tilde{f}_{\Lambda}$ on
$\mathbb{R}$ as $\Lambda$ goes to infinity (see Propositions 2.1, 2.6, 2.7 and
2.10). In order to use the convergence of $f_{\Lambda}$ and
$\tilde{f}_{\Lambda}$ in space $\mathcal{C}^{m}[a,b]$ for $m\in\mathbb{N}$ and
$a<b$, we firstly have to prove the uniform boundedness of $\phi_{\Lambda}$
(the solution of (1.21) and (1.26)) with respect to $\Lambda$ (see Lemmas 3.1
and 3.4). Here we notice that $\rho_{0}$ may be any nonzero smooth function
and the boundary condition may be the Robin boundary condition but not
Dirichlet boundary condition so one cannot simply use the maximum principle on
(1.19) and (1.20) to obtain the uniform boundedness of $\phi_{\Lambda}$.

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

######  Lemma 3.1.

There exist positive constants $M_{1}\geq 1$, $M_{2}$, and $M_{3}$ independent
of $\Lambda$ such that

  * (i)

$\displaystyle\max_{x\in\overline{\Omega}}c_{0,\Lambda}(\phi_{\Lambda}(x))\leq
M_{1}$ for $\Lambda\geq 1$.

  * (ii)

$\|\phi_{\Lambda}\|_{L^{\infty}(\Omega)}\leq M_{2}$ for $\Lambda\geq 1$.

  * (iii)

$\displaystyle\min_{x\in\overline{\Omega}}c_{0,\Lambda}(\phi_{\Lambda}(x))\geq
M_{3}$ for $\Lambda\geq 1$.

###### Proof.

Since $z_{i}z_{j}<0$ for some $i,j\in\\{1,\dots,N\\}$, then
$\displaystyle\lim_{\phi\to\pm\infty}\sum_{i=1}^{N}e^{\bar{\mu}_{i}-z_{i}\phi}=\infty$,
and

$K:=(\tilde{\mu}_{0}+\hat{\mu}_{0})/\left(1-\gamma+\frac{\gamma}{Z}\min_{\phi\in\mathbb{R}}\sum_{i=1}^{N}\exp(\bar{\mu}_{i}-z_{i}\phi)\right)<\infty.$

Let $M_{1}=\max\\{K,1\\}$ and $\Lambda\geq 1$. We claim that
$c_{0,\Lambda}(\phi_{\Lambda}(x))\leq M_{1}$. Suppose that
$\Omega_{1}=\\{x\in\overline{\Omega}:c_{0,\Lambda}(\phi_{\Lambda}(x))>1\\}$ is
nonempty. Otherwise, due to $M_{1}\geq 1$,
$c_{0,\Lambda}(\phi_{\Lambda}(x))\leq 1\leq M_{1}$ for
$x\in\overline{\Omega}$, which is trivial. By (1.17), we have
$c_{0,\Lambda}(\phi_{\Lambda}(x))\leq\frac{\tilde{\mu}_{0}+\Lambda^{-1}\hat{\mu}_{0}}{H(\gamma,z_{1},\dots,z_{N},\phi_{\Lambda}(x))}\leq
K\leq M_{1}$ 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}$ and
$\Lambda\geq 1$. Therefore, we complete the proof of (i).

To prove (ii), let $\psi$ be the solution of equation
$-\nabla\cdot(\varepsilon\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(\varepsilon\nabla\bar{\phi}_{\Lambda})=f_{\Lambda}(\phi_{\Lambda})\quad\text{in}~{}\Omega.$
(3.1)

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, it suffices to show that
$\displaystyle\max_{\overline{\Omega}}\bar{\phi}_{\Lambda}\leq M_{2}$ and
$\displaystyle\min_{\overline{\Omega}}\bar{\phi}_{\Lambda}\geq-M_{2}$ for
$\Lambda\geq 1$, where $M_{2}$ is a positive constant independent of
$\Lambda$.

Now we prove that
$\displaystyle\max_{\overline{\Omega}}\bar{\phi}_{\Lambda}\leq M_{2}$ for
$\Lambda\geq 1$, where $M_{2}$ 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 the sake of simplicity, 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 (3.1) with $\nabla\bar{\phi}_{k}(x_{k})=0$,
$\Delta\bar{\phi}_{k}(x_{k})\leq 0$ and function $\varepsilon$ is positive, we
have

$0\leq-\nabla\varepsilon(x_{k})\cdot\nabla\bar{\phi}_{k}(x_{k})-\varepsilon(x_{k})\Delta\bar{\phi}_{k}(x_{k})=f_{k}(\phi_{k}(x_{k})).$
(3.2)

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

$0<\sum_{z_{i}<0}(-z_{i})c_{i,k}(\psi(x_{k})+\bar{\phi}_{k}(x_{k}))\leq\sum_{z_{i}>0}z_{i}c_{i,k}(\psi(x_{k})+\bar{\phi}_{k}(x_{k})).$
(3.3)

By Proposition 2.1(iii)(b) and (3.3),
$\displaystyle\lim_{k\to\infty}\sum_{z_{i}<0}(-z_{i})c_{i,k}(\psi(x_{k})+\bar{\phi}_{k}(x_{k}))=0$.
This leads a contradiction with Proposition 2.1(iii)(b). Hence we complete the
proof to show that
$\displaystyle\max_{\overline{\Omega}}\bar{\phi}_{\Lambda}\leq M_{2}$ for
$\Lambda\geq 1$, where $M_{2}$ is a positive constant independent of
$\Lambda$. The proof for the case that
$\displaystyle\min_{\overline{\Omega}}\bar{\phi}_{\Lambda}\geq-M_{2}$ for
$\Lambda\geq 1$ is similar to the case that
$\displaystyle\max_{\overline{\Omega}}\bar{\phi}_{\Lambda}\leq M_{2}$ for
$\Lambda\geq 1$ so we omit it here. Thus, the proof of (ii) is complete.

Finally, we state the proof of (iii). Suppose to the 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 (1.17), $c_{0,\Lambda}(\phi_{\Lambda})>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 (1.17)
with $\phi=\phi_{\Lambda_{k}}(x_{k})$, we get

$H(\gamma,z_{1},\dots,z_{N},\phi_{\Lambda_{k}}(x_{k}))c_{0,\Lambda_{k}}(\phi_{\Lambda_{k}}(x_{k}))>\tilde{\mu}_{0}+\frac{\hat{\mu}_{0}}{\Lambda_{k}}\quad\text{for}~{}k>N_{1}.$
(3.4)

Taking $k\to\infty$, (ii) and (3.4) imply $0\geq\tilde{\mu}_{0}$, which
contradicts assumption (A2). Therefore, we complete the proof of (iii). The
proof of Lemma 3.1 is complete. ∎

######  Remark 3.2.

By Lemma 3.1 and (1.16), we have $M_{4}\leq
c_{i,\Lambda}(\phi_{\Lambda}(x))\leq M_{5}$ for $x\in\overline{\Omega}$,
$\Lambda\geq 1$ and $i=0,1,\dots,N$, where $M_{4}$ and $M_{5}$ 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})$,
there exists a positive constant $M_{6}$ independent of $\Lambda$ such that
$\|f_{\Lambda}(\phi_{\Lambda})\|_{L^{\infty}(\overline{\Omega})}\leq M_{6}$
for $\Lambda\geq 1$.

###  3.2 Convergence of $\phi_{\Lambda}$ under (A1) and (A2)

Since $\phi_{\Lambda}$ is the solution of equation
$-\nabla\cdot(\varepsilon\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. [1, 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.2, 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
Corollary 2.5, we have the convergence of function $f_{\Lambda_{k}}$ to
$f^{*}$ in $\mathcal{C}^{m}[-M_{2},M_{2}]$ for $m\in\mathbb{N}$, so $\phi^{*}$
satisfies (1.20) in a weak sense, where the positive constant $M_{2}$ comes
from Lemma 3.1(i). 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}(\overline{\Omega})$ for
$\alpha\in(0,1)$, and
$\displaystyle\lim_{k\to\infty}\|w_{k}\|_{\mathcal{C}^{1,\alpha}(\overline{\Omega})}=0$.
Moreover, $w_{k}$ satisfies $-\nabla\cdot(\varepsilon\nabla
w_{k})=f_{k}(w_{k}+\phi^{*})-f^{*}(\phi^{*})$ 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. [12, Theorem 6.30]) with
the mathematical induction, we get

$\displaystyle{\color[rgb]{1,1,1}\leq}\,\|w_{k}\|_{\mathcal{C}^{m+2,\alpha}(\overline{\Omega})}\leq
C\|f_{k}(w_{k}+\phi^{*})-f^{*}(\phi^{*})\|_{\mathcal{C}^{m,\alpha}(\overline{\Omega})}$
(3.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}(\overline{\Omega})}+\|c_{i}^{*}(w_{k}+\phi^{*})-c_{i}^{*}(\phi^{*})\|_{\mathcal{C}^{m,\alpha}(\overline{\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(i) and induction
hypothesis
$\displaystyle\lim_{k\to\infty}\|w_{k}\|_{\mathcal{C}^{m,\alpha}(\overline{\Omega})}=0$,
we may use (3.5) to get
$\displaystyle\lim_{k\to\infty}\|w_{k}\|_{\mathcal{C}^{m+2,\alpha}(\overline{\Omega})}=0$,
i.e.
$\displaystyle\lim_{k\to\infty}\|\phi_{\Lambda_{k}}-\phi^{*}\|_{C^{m+2,\alpha}(\overline{\Omega})}=0$
for $m\in\mathbb{N}$ and $\alpha\in(0,1)$. Therefore, $\phi^{*}$ is the
solution of (1.26) with the Robin boundary condition (1.15).

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

######  Claim 3.3.

For any $m\in\mathbb{N}$, we have
$\displaystyle\lim_{\Lambda\to\infty}\|\phi_{\Lambda}-\phi^{*}\|_{\mathcal{C}^{m}(\overline{\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 (1.31) with the Robin boundary
condition (1.15). Now we want to prove that $\phi_{1}^{*}\equiv\phi_{2}^{*}$.
Let $u=\phi_{1}^{*}-\phi_{2}^{*}$. Subtracting (1.31) with
$\phi^{*}=\phi_{2}^{*}$ from that with $\phi^{*}=\phi_{1}^{*}$, we obtain
$-\nabla\cdot(\varepsilon\nabla u)=c_{1}(x)u$ in $\Omega$, where
$c_{1}(x)=\frac{f^{*}(\phi_{1}^{*}(x))-f^{*}(\phi_{2}^{*}(x))}{\phi_{1}^{*}(x)-\phi_{2}^{*}(x)}$
if $\phi_{1}^{*}(x)\neq\phi_{2}^{*}(x)$;
$\frac{\mathrm{d}f^{*}}{\mathrm{d}\phi}(\phi_{1}^{*}(x))$ if
$\phi_{1}^{*}(x)=\phi_{2}^{*}(x)$. By Proposition 2.6, we have $c_{1}<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(\varepsilon\nabla u)+c_{1}(x)u=0$ with
$c_{1}<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. ∎

###  3.3 Proof of Theorem 1.6

######  Lemma 3.4.

There exist positive constants $M_{7}\geq 1$, $M_{8}$, and $M_{9}$ independent
of $\Lambda$ such that

  * (i)

$\displaystyle\max_{x\in\overline{\Omega}}c_{0,\Lambda}(\phi_{\Lambda}(x))\leq
M_{7}$ for $\Lambda\geq 1$.

  * (ii)

$\|\phi_{\Lambda}\|_{L^{\infty}(\Omega)}\leq M_{8}$ for $\Lambda\geq 1$.

  * (iii)

$\displaystyle\min_{x\in\overline{\Omega}}c_{0,\Lambda}(\phi_{\Lambda}(x))\geq
M_{9}$ for $\Lambda\geq 1$.

###### Proof.

The proof is identical to the proofs of Lemmas 3.1, 3.2, and 3.3 in [33] so we
omit it here. ∎

######  Remark 3.5.

By Lemma 3.4 and (1.24), we have $M_{10}\leq
c_{i,\Lambda}(\phi_{\Lambda}(x))\leq M_{11}$ for $x\in\overline{\Omega}$,
$\Lambda\geq 1$ and $i=0,1,\dots,N$, where $M_{10}$ and $M_{11}$ are positive
constants independent of $\Lambda$. Moreover, due to
$\displaystyle\tilde{f}_{\Lambda}(\phi_{\Lambda})=\sum_{i=1}^{N}z_{i}c_{i,\Lambda}(\phi_{\Lambda})$,
$\|\tilde{f}_{\Lambda}(\phi_{\Lambda})\|_{L^{\infty}(\overline{\Omega})}\leq
M_{12}$ for $\Lambda\geq 1$, where $M_{12}$ is a positive constant independent
of $\Lambda$.

For the proof of convergence of $\phi_{\Lambda}$ in
$\mathcal{C}^{m}(\overline{\Omega})$ for $m\in\mathbb{N}$, we can use Remark
3.5 and the same method in section 3.2 to complete the proof of Theorem 1.6.

##  4 Numerical methods

In this section, we introduce numerical methods to solve the PB-steric
equations (1.13)–(1.14) with the Robin boundary condition (1.15) to get the
oscillation of total ionic charge densities $f_{\Lambda}$ for
$\Lambda=0.5,1,2,4$ (see Figure 1). Throughout this section, we consider one
dimensional domain $\Omega=(-1,1)$ and assume $N=3$, $z_{0}=0$, $z_{1}=1$,
$z_{2}=-1$, $z_{3}=2$, $\tilde{\mu}_{0}=1$, $\hat{\mu}_{0}=0$, $\rho_{0}\equiv
0$, $\eta=\varepsilon=0.1$, and $\phi_{bd}(\pm 1)=\pm 10$. The assumptions
(A1)′′′ and (A2) are satisfied when the parameters
$1-\gamma=\gamma_{1}=\gamma_{2}=0.01$, $\gamma_{3}=0.97$,
$\bar{\mu}_{0}=\bar{\mu}_{1}=0$, and $\bar{\mu}_{2}=\bar{\mu}_{3}=-5$.

To solve (1.13), we use the command fsolve in Matlab to perform the trust-
region algorithm (cf. [11]) on the following discretized nonlinear equation:

$\displaystyle\ln(\boldsymbol{c}_{i,\Lambda})+z_{i}\boldsymbol{\phi}+\Lambda\sum_{j=0}^{N}g_{ij}\boldsymbol{c}_{j,\Lambda}=\Lambda\tilde{\mu}_{i}+\hat{\mu}_{i}\quad\text{for}~{}i=0,1,\dots,N,$
(4.1)

where $\boldsymbol{\phi}=[\phi_{0},\dots,\phi_{L}]^{\mathsf{T}}$ is the
regular partition of the interval $[-10,10]$ with $L=1024$, and
$\boldsymbol{c}_{i,\Lambda}=[c_{i,\Lambda}(\phi_{0}),\dots,c_{i,\Lambda}(\phi_{L})]^{\mathsf{T}}$.
Moreover, in order to get the solution $\phi_{\Lambda}$ of (1.13)–(1.14), we
employ the Legendre–Gauss–Lobatto (LGL) points $\\{x_{k}\\}_{k=0}^{L}$ (cf.
[10]) as the partition of the interval $[-1,1]$ with $L=256$ to discrtize
(1.13)–(1.14) as the following algebraic equations:

$\displaystyle-\text{D}_{L+1}\boldsymbol{\varepsilon}\text{D}_{L+1}\boldsymbol{\phi}_{\Lambda}=\boldsymbol{\rho}_{0}+\sum_{i=1}^{N}z_{i}\boldsymbol{c}_{i,\Lambda},$
(4.2)
$\displaystyle\ln(\boldsymbol{c}_{i,\Lambda})+z_{i}\boldsymbol{\phi}_{\Lambda}+\Lambda\sum_{j=0}^{N}g_{ij}\boldsymbol{c}_{j,\Lambda}=\Lambda\tilde{\mu}_{i}+\hat{\mu}_{i}\quad\text{for}~{}i=0,1,\dots,N,$
(4.3)

where $\text{D}_{L+1}=[d_{ij}]_{0\leq i,j\leq L}$ is the matrix satisfying
$\text{D}_{L+1}\boldsymbol{\psi}=[\psi^{\prime}(x_{0}),\dots,\psi^{\prime}(x_{L})]^{\mathsf{T}}$
for $\boldsymbol{\psi}=[\psi(x_{0}),\dots,\psi(x_{L})]^{\mathsf{T}}$,
$\boldsymbol{\varepsilon}=[\varepsilon(x_{0}),\dots,\varepsilon(x_{L})]^{\mathsf{T}}$,
$\boldsymbol{\rho}_{0}=[\rho_{0}(x_{0}),\dots,\rho_{0}(x_{L})]^{\mathsf{T}}$,
$\boldsymbol{\phi}_{\Lambda}=[\phi_{\Lambda}(x_{0}),\dots,\phi_{\Lambda}(x_{L})]^{\mathsf{T}}$,
and
$\boldsymbol{c}_{i,\Lambda}=[c_{i,\Lambda}(\phi_{\Lambda}(x_{0})),\dots,c_{i,\Lambda}(\phi_{\Lambda}(x_{L}))]^{\mathsf{T}}$
for $i=0,1,\dots,N$.

For (A1)′′′ and (A2), equations (4.2) and (4.3) can be denoted as

$\displaystyle-\text{D}_{L+1}\boldsymbol{\varepsilon}\text{D}_{L+1}\boldsymbol{\phi}_{\Lambda}=\boldsymbol{\rho}_{0}+\sum_{i=1}^{N}z_{i}\boldsymbol{c}_{0,\Lambda}\circ\exp(\bar{\mu}_{i}-z_{i}\boldsymbol{\phi}_{\Lambda}),$
(4.4)
$\displaystyle\ln(\boldsymbol{c}_{0,\Lambda})+\Lambda\left(1-\gamma+\sum_{j=1}^{N}\gamma_{j}\exp(\bar{\mu}_{j}-z_{j}\boldsymbol{\phi}_{\Lambda})\right)\circ\boldsymbol{c}_{0,\Lambda}=\Lambda\tilde{\mu}_{0}+\hat{\mu}_{0}.$
(4.5)

In addition, we discretize the Robin boundary condition (1.15) as

$\phi(x_{0})-\eta\sum_{k=0}^{L}d_{0k}\phi(x_{k})=\phi_{bd}(-1),\quad\phi(x_{L})+\eta\sum_{k=0}^{L}d_{Lk}\phi(x_{k})=\phi_{bd}(1),$
(4.6)

which replaces the first and last equations of (4.4) and (4.6). Then using the
command fsolve with the initial data $\boldsymbol{\phi}_{\Lambda}^{(0)}\equiv
0$ and $\boldsymbol{c}_{0,\Lambda}^{(0)}\equiv 1$, we obtain the profiles of
$f_{\Lambda}(\phi)$ and $f_{\Lambda}(\phi_{\Lambda}(x))$ in Figure 1.

##  5 Conclusions

In this paper, we derive the PB-steric equations with a parameter $\Lambda$,
and prove that as $\Lambda\to\infty$, the solution of PB-steric equations
(1.13)–(1.14) converges to the solution of modified PB equations (1.3)–(1.4)
and (1.5)–(1.7) under different assumptions of steric effects and chemical
potentials. Equations (1.1)–(1.2), (1.3)–(1.4), and (1.5)–(1.7) were derived
based on mean-field approximation (cf. [6, 27]), whereas equations
(1.13)–(1.14) are derived via adding the approximate Lennard–Jones potential.
We compare equations (1.1)–(1.7) and (1.13)–(1.14) in the following table to
show that the PB-steric equations (1.13)–(1.14) can be regarded as a general
model of PB equations.

Model | Total ionic charge density $\sum_{i=1}^{N}z_{i}c_{i}(\phi)$ |  $c_{i}(\phi)$ ($i=0,1,\dots,N$)
---|---|---
(1.1)–(1.2) | decreasing |  $c_{1}$, $c_{2}$ are monotonic,
$c_{0}\equiv 0$ ($N=2$).
(1.3)–(1.4) | decreasing |  $c_{1}$, $c_{2}$ are monotonic,
$c_{0}\equiv 0$ ($N=2$).
(1.5)–(1.7) | decreasing | Some of $c_{i}$’s are oscillatory;
the others are monotonic.
($N\geq 2$)
(1.13)–(1.14) |  oscillatory (Remark 1.5) | Some of $c_{i}$’s are oscillatory;
the others are monotonic.
($N=3$)
decreasing (Propositions 2.6 and 2.10) | Some of $c_{i}$’s are oscillatory;
the others are monotonic.
($N\geq 2$)

##  Appendix A Derivation of (1.10) and (1.11)

In [30], the approximate LJ potential $\Psi_{ij,\sigma}$ between the $i$th ion
and the $j$th ion can be expressed as
$\Psi_{ij,\sigma}(x)=(\Psi_{ij}\chi_{\sigma}\ast\varphi_{\sigma})(x)$ for
$x\in\mathbb{R}^{d}$ and $i,j=0,1,\dots,N$, where $\ast$ denotes the standard
convolution, $\Psi_{ij}$ is the conventional LJ potential between $i$th and
$j$th ions defined by
$\displaystyle\Psi_{ij}(x)=\frac{\epsilon_{ij}(a_{i}+a_{j})^{12}}{2|x|^{12}}$
for $x\in\mathbb{R}^{d}$ and $i,j=0,1,\dots,N$, and $\varphi_{\sigma}(x)$ is
the spatially band-limited function defined by
$\varphi_{\sigma}(x)=\mathcal{F}^{-1}(1-\chi_{\sigma^{-\gamma}}(\xi))$. Here
$\epsilon_{ij}$ denotes energy coupling constant between the $i$th and the
$j$th ions, $a_{0}$ is the radius of solvent molecules, $a_{i}$ is the radius
of the $i$th ion. In addition, $\chi_{\sigma}$ and $\chi_{\sigma^{-\gamma}}$
are the characteristic function of the exterior ball
$\mathbb{R}^{d}-\overline{B_{\sigma}(0)}$ and
$\mathbb{R}^{d}-\overline{B_{\sigma^{-\gamma}}(0)}$ for $0<\gamma<1$,
respectively, and $\mathcal{F}^{-1}$ is the inverse Fourier transform. More
precisely, the characteristic function $\chi_{s}$ is defined by
$\chi_{s}(z)=0~{}\text{if}~{}|z|\leq s;1~{}\text{if}~{}|z|>s$. With the
approximate LJ potential, we may define the approximate energy functional
$E_{LJ,\sigma}$ and show that (cf. [30])

$E_{LJ,\sigma}[c_{i},c_{j}]:=\Lambda
g_{ij}\int_{\mathbb{R}^{d}}\\!c_{i}(x)c_{j}(x)\,\mathrm{d}x\sim\int_{\mathbb{R}^{d}}\\!c_{i}(x)(\Psi_{ij}\ast
c_{j})(x)\,\mathrm{d}x=E_{LJ}[c_{i},c_{j}],$ (A.1)

where $\Lambda=a^{12}\sigma^{d-12}$ satisfies $\displaystyle\lim_{\sigma\to
0^{+}}\Lambda=\infty$, and $\displaystyle
g_{ij}=\frac{\omega_{d}\epsilon_{ij}}{12-d}\frac{(a_{i}+a_{j})^{12}}{a^{12}}$
denotes a dimensionless quantity. Here $a=\min\\{a_{i}:i=0,1\dots,N\\}$ and
$\omega_{d}$ is the surface area of a $d$ dimensional unit ball. To derive
(1.10) and (1.11), we let energy functional

$E[c_{0},c_{1},\dots,c_{N},\phi]=E_{PB}[c_{0},c_{1},\dots,c_{N},\phi]+\frac{1}{2}\sum_{i,j=0}^{N}E_{LJ,\sigma}[c_{i},c_{j}],$

where $E_{LJ,\sigma}$ is defined in (A.1) and

$\displaystyle
E_{PB}[c_{0},c_{1},\cdots,c_{N},\phi]=\int_{\mathbb{R}^{d}}\left[-\frac{1}{2}\varepsilon\left|\nabla\phi\right|^{2}+k_{B}T\sum_{i=0}^{N}c_{i}(\ln
c_{i}-1)+4\pi\left(\rho_{0}+\sum_{i=1}^{N}z_{i}e_{0}c_{i}\right)\phi\right]\mathrm{d}x.$

Notice that $E_{PB}$ is the energy functional of conventional PB equation with
the form
$\displaystyle-\nabla\cdot(\varepsilon\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$. Here $\mu_{i}$ is the
chemical potential. Besides, $E_{LJ,\sigma}$ is the energy functional of the
approximate LJ potentials (cf. [30]), and we may derive the equations

$\displaystyle-\nabla\cdot(\varepsilon\nabla\phi)=4\pi\rho_{0}+4\pi\sum_{i=1}^{N}z_{i}e_{0}c_{i},$
(A.2) $\displaystyle
k_{B}T\ln(c_{i})+z_{i}e_{0}\phi+\Lambda\sum_{j=0}^{N}g_{ij}c_{j}=\mu_{i}\quad\text{for}~{}i=0,1,\dots,N,$
(A.3)

by $\delta E/\delta\phi=0$ and $\delta E/\delta c_{i}=\mu_{i}$ for
$i=0,1,\dots,N$ when matrix $(g_{ij})$ is symmetric. From (A.3), constant
$\mu_{i}$ can be expressed by

$\mu_{i}=\Lambda\tilde{\mu}_{i}+\hat{\mu}_{i}\quad\text{for}~{}i=0,1,\dots,N.$
(A.4)

Therefore, by (A.2)–(A.4), we have (1.10)–(1.11).

##  Appendix B An example of oscillatory $f_{\Lambda}$ for Remark 1.5

Under the assumptions (A1)′′′ and (A2), system (1.13) can be solved by

$c_{i,\Lambda}(\phi)=\frac{\displaystyle
W_{0}\left(\Lambda\left(1-\gamma+\sum_{j=1}^{3}\gamma_{j}e^{\bar{\mu}_{j}-z_{j}\phi}\right)e^{\Lambda\tilde{\mu}_{0}+\hat{\mu}_{0}}\right)}{\displaystyle\Lambda\left(1-\gamma+\sum_{j=1}^{3}\gamma_{j}e^{\bar{\mu}_{j}-z_{j}\phi}\right)}e^{\bar{\mu}_{i}-z_{i}\phi}$

for $\phi\in\mathbb{R}$ and $i=0,1,2,3$. Suppose that $\Lambda=4$, $z_{0}=0$,
$z_{1}=-z_{2}=1$ $z_{3}=2$, $\bar{\mu}_{0}=\bar{\mu}_{1}=0$,
$\bar{\mu}_{2}=\bar{\mu}_{3}=-5$, $\tilde{\mu}_{0}=1$, $\hat{\mu}_{0}=0$,
$\gamma=0.99$, $\gamma_{1}=\gamma_{2}=0.01$, and $\gamma_{3}=0.97$. Recall the
derivative of the Lambert W function (cf. [15, 34]) is
$\displaystyle\frac{\mathrm{d}W_{0}}{\mathrm{d}x}=\frac{W_{0}(x)}{x(1+W_{0}(x))}$
for $x>-1/e$. Differentiating $c_{i,4}$ with respect to $\phi$, we obtain

$\displaystyle\frac{\mathrm{d}c_{i,4}}{\mathrm{d}\phi}$
$\displaystyle=\frac{W_{0}(0.04e^{4}(1+e^{-\phi}+e^{-5+\phi}+97e^{-5-2\phi}))}{0.04(1+e^{-\phi}+e^{-5+\phi}+97e^{-5-2\phi})}\exp(\bar{\mu}_{i}-z_{i}\phi)$
(B.1)
$\displaystyle~{}~{}~{}\times\left(\frac{(e^{-\phi}-e^{-5+\phi}+194e^{-5-2\phi})W_{0}(0.04e^{4}(1+e^{-\phi}+e^{-5+\phi}+97e^{-5-2\phi}))}{(1+e^{-\phi}+e^{-5+\phi}+97e^{-5-2\phi})(1+W_{0}(0.04e^{4}(1+e^{-\phi}+e^{-5+\phi}+97e^{-5-2\phi})))}-z_{i}\right)$

for $\phi\in\mathbb{R}$ and $i=0,1,2,3$. Let
$K_{1}=0.04(1+e^{2}+e^{-7}+97e^{-1})$ and $K_{2}=e^{2}-e^{-7}+194e^{-1}$. Then
we have

$\frac{\mathrm{d}f_{4}}{\mathrm{d}\phi}(-2)=\frac{W_{0}(e^{4}K_{1})}{K_{1}}\left[\frac{K_{2}W_{0}(e^{4}K_{1})}{25K_{1}(1+W_{0}(e^{4}K_{1}))}(e^{2}-e^{-7}+2e^{-1})-(e^{2}+e^{-7}+4e^{-1})\right].$
(B.2)

It is easy to check that $K_{2}/K_{1}>44$ and
$3.3e^{3.3}<e^{4}K_{1}<3.4e^{3.4}$. Due to the monotonic increase of $W_{0}$,
we have $3.3<W_{0}(e^{4}K_{1})<3.4$. Along with (B.2), we arrive at the fact
$\displaystyle\frac{\mathrm{d}f_{4}}{\mathrm{d}\phi}(-2)>0$. On the other
hand, from (B.1), we apply L’Hôspital’s law to get
$\displaystyle\lim_{\phi\to\infty}\frac{\mathrm{d}c_{1,4}}{\mathrm{d}\phi}(\phi)=0$,
$\displaystyle\lim_{\phi\to\infty}\frac{\mathrm{d}c_{2,4}}{\mathrm{d}\phi}(\phi)=25$,
$\displaystyle\lim_{\phi\to\infty}\frac{\mathrm{d}c_{3,4}}{\mathrm{d}\phi}(\phi)=0$,
which implies
$\displaystyle\lim_{\phi\to\infty}\frac{\mathrm{d}f_{4}}{\mathrm{d}\phi}(\phi)=-25$.
This shows $\displaystyle\frac{\mathrm{d}f_{4}}{\mathrm{d}\phi}(\phi)<0$ when
$\phi$ is sufficiently large. Therefore, $f_{4}$ is oscillatory and the
profile of $f_{4}$ is in Figure 1(a).

## Acknowledgments

The authors express their thanks to reviewers for the careful and insightful
review.

## References

  * [1] S. Agmon, A. Douglis, and L. Nirenberg, Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions. I, Commun. Pure Appl., 12 (1959), pp. 623–727.
  * [2] P. B. Balbuena and K. E. Gubbins, Theoretical interpretation of adsorption behavior of simple fluids in slit pores, Langmuir, 9 (1993), pp. 1801–1814.
  * [3] M. Z. Bazant, K. T. Chu, and B. J. Bayly, Current-voltage relations for electrochemical thin films, SIAM J. Appl. Math., 65 (2005), pp. 1463–1484.
  * [4] M. Z. Bazant, B. D. Storey, and A. A. Kornshev, Double layer in ionic liquids: Overscreening versus crowding, Phys. Rev. Lett., 106 (2011), 046102.
  * [5] J. J. Bikerman, Structure and capacity of electrical double layer, Philos. Mag., 7 (1942), pp. 384–397.
  * [6] I. Borukhov and D. Andelman, Steric effects in electrolytes: A Modified Poisson-Boltzmann equation, Phys. Rev. Lett., 79 (1997), pp. 435–438.
  * [7] K. Deimling, Nonlinear Functional Analysis, Springer-Verlag, Berlin, 1973.
  * [8] B. Eisenberg, Proteins, channels and crowded ions, Biophys. Chem., 100 (2003), pp. 507–517.
  * [9] B. Eisenberg, Interacting ions in biophysics: Real is not ideal, Biophys. J., 104 (2013), pp. 1849–1866.
  * [10] G. I. El-Baghdady, M. S. El-Azab, and W. S. El-Beshbeshy, Legendre-Gauss-Lobatto pseudo-spectral method for one-dimensional advection-diffusion equation, Sohag J. Math., 2 (2015), pp. 29–35.
  * [11] J. Fan and J. Pan, An improved trust region algorithm for nonlinear equations, Comput. Optim. Appl., 48 (2011), pp. 59–70.
  * [12] D. Gilbarg and N. S. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer-Verlag, Berlin, 1977.
  * [13] D. Gillespie, A review of steric interactions of ions: Why some theories succeed and others fail to account for ion size, Microfluid Nanofluid, 18 (2015), pp. 717–738.
  * [14] P. Grochowski and J. Trylska, Continuum molecular electrostatics, salt effects, and counterion binding—A Review of the Poisson–Boltzmann theory and its modifications, Biopolymers, 89 (2007), pp. 93–113.
  * [15] A. Hoorfar and M. Hassani, Inequalities on the Lambert W function and hyperpower function, JIPAM, 9 (2008), pp. 93–113.
  * [16] T.-L. Horng, T.-C. Lin, C. Liu, and B. Eisenberg, PNP equations with steric effects: A model of ion flow through channels, J. Phys. Chem. B, 116 (2012), pp. 11422–11441.
  * [17] B. Huang, S. Maset, and K. Bohinc, Interaction between charged cylinders in electrolyte solution; Excluded volume effect, J. Phys. Chem. B, 121 (2017), pp. 9013–9023.
  * [18] I. Kalcher, J. C. F. Schulz, and J. Dzubiella, Ion-specific excluded-volume correlations and solvation forces, Phys. Rev. Lett., 104 (2010), 097802.
  * [19] E. Kierlik and M. L. Rosinberg. Density-functional theory for inhomogeneous fluids: Adsorption of binary mixtures, Phys. Rev. A 44, (1991), pp. 5025–5037.
  * [20] B. K. Klein, C. F. Anderson, and M. T. Record, Comparison of Poisson-Boltzmann and condensation model expressions for the colligative properties of cylindrical polyions, Biopolymers, 20 (1981), pp. 2263–2280.
  * [21] C. Köhn, D. van Laethem, J. Deconinck, and A. Hubin, A simulation study of steric effects on the anodic dissolution at high current densities, Mater. Corros., 72 (2020), pp. 610–619.
  * [22] A. A. Kornyshev, Double-layer in ioninc liquids: Paradigm change, J. Phys. Chem. B, 111 (2007), pp. 5545–5557.
  * [23] S. G. Krantz and H. R. Parks, The Implicit Function Theorem: History, Theory, and Applications, Birkhäuser, Basel, 2002.
  * [24] G. Lamm, The Poisson-Boltzmann Equation, In Rveiews in Computational Chemistry, Rev. Compuy. Chem. 19, John Wiley & Sons, New York, 2003.
  * [25] J. W. Lee, R. H. Nilson, J. A. Templeton, S. K. Griffiths, A. Kung, and B. M. Wong, Comparison of molecular dynamics with classical density functional and Poisson-Boltzmann theories of the electric double layer in nanochannels, J. Chem. Theory Comput., 8 (2012), pp. 2012–2022.
  * [26] B. Leimkuhler and C. Matthews, Molecular Dynamics with Deterministic and Stochastic Numerical Methods, Springer, New York, 2015.
  * [27] B. Li, Continuum electrostatics for ionic solutions with non-uniform ionic sizes, Nonlinearity, 22 (2009), pp. 811–833.
  * [28] B. Li, Minimization of electrostatic free energy and the Poisson-Boltzmann equation for molecular solvation with implicit solvent, SIAM J. Math. Anal., 40 (2009), pp. 2536–2566.
  * [29] B. Li, P. Liu, Z. Xu, and S. Zhou, Ionic size effects: Generalized Boltzmann distributions, counterion stratification and modified Debye length, Nonlinearity, 26 (2013), pp. 2899–2922.
  * [30] T.-C. Lin and B. Eisenberg, A new approach to the Lennard-Jones potential and a new model: PNP-steric equations, Commun. Math., 12 (2014), pp. 149–173.
  * [31] B. Z. Lu, Y. C. Zhou, M. J. Holst, and J. A. McCammon, Recent progress in numerical methods for the Poisson-Boltzmann equation in biophysical applications, Commun. Comput. Phys., 3 (2008), pp. 973–1009.
  * [32] B. Lu and 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), pp. 2475–2485.
  * [33] J.-H. Lyu and T.-C. Lin, On the Asymptotic Limit of Steady State Poisson-Nernst-Planck Equations with Steric Effects, preprint, https://arxiv.org/abs/2205.12513v1.2022.
  * [34] I. Mező, The Lambert W Function: Its Generalizations and Applications, Chapman and Hall, London, 2022.
  * [35] Y. Mori, J. W. Jerome, and C. S. Peskin, A three-dimensional model of cellular electrical activity, Bull. Inst. Math. Acad. Sin., 2 (2007), pp. 367–390.
  * [36] R. H. Nilson and S.K. Griffiths, Influence of atomistic physics on electro-osmotic flow: An analysis based on density functional theory, J. Chem. Phys., 125 (2006), 164510.
  * [37] D. C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, Cambridge, 1995.
  * [38] Z. Tang, L. E. Scriven, and H. T. Davis, A three-component model of the electrical double layer, J. Chem. Phys., 97 (1992), pp. 494–503.
  * [39] F. Ziebert, M. Z. Bazant, and D. Lacoste, Effective zero-thickness model for a conductive membrane driven by an electric field, Phys. Rev. E, 81 (2010), 031912.
  * [40] S. Zhou, Z. Wang, and B. Li, Mean-field description of ionic size effects with nonuniform ionic sizes: A numerical approach, Phys. Rev. E, 84 (2011), 021901.