# A dynamic mass transport method for Poisson-Nernst-Planck equations

Hailiang Liu and Wumaier Maimaitiyiming  † Department of Mathematics, Iowa
State University, USA<EMAIL_ADDRESS>Department of Mathematics, University
of California Los Angeles, USA<EMAIL_ADDRESS>
###### Abstract.

A dynamic mass-transport method is proposed for approximately solving the
Poisson–Nernst–Planck (PNP) equations. The semi-discrete scheme based on the
JKO type variational formulation naturally enforces solution positivity and
the energy law as for the continuous PNP system. The fully discrete scheme is
further formulated as a constrained minimization problem, shown to be
solvable, and satisfy all three solution properties (mass conservation,
positivity and energy dissipation) independent of time step size or the
spatial mesh size. Numerical experiments are conducted to validate convergence
of the computed solutions and verify the structure preserving property of the
proposed scheme.

###### Key words and phrases:

PNP equations, Optimal transport, Wasserstein distance, Positivity, Energy
dissipation

###### 1991 Mathematics Subject Classification:

35Q92, 65N08, 65N12

##  1\. Introduction

In this paper, we consider a time-dependent system of Poisson-Nernst-Planck
(PNP) equations. Such system has been widely used to describe charge transport
in diverse applications such as biological membrane channels [37, 6, 9],
electrochemical systems [2], and semiconductor devices [29, 34].

PNP equations consist of Nernst–Planck (NP) equations that describe the drift
and diffusion of ion species, and the Poisson equation that describes the
electrostatic interaction. Such mean field approximation of diffusive ions
admits several variants, and in non-dimensional form we consider the following

$\displaystyle\partial_{t}\rho_{i}=\nabla\cdot\left[D_{i}(x)\left(\nabla\rho_{i}+z_{i}\rho_{i}\nabla\phi\right)\right],\quad
x\in\Omega\subset\mathbb{R}^{d},\quad t>0,$ (1.1a)
$\displaystyle-\nabla\cdot(\epsilon(x)\nabla\phi)=f(x)+\sum_{i=1}^{s}z_{i}\rho_{i},$
(1.1b)

subject to initial data $\rho_{i}(x,0)=\rho_{i}^{in}(x)\geq 0$
($i=1,\cdots,s$) and appropriate boundary conditions to be specified in
section 2. The equations are valid in a bounded domain $\Omega$ with boundary
$\partial\Omega$ and for time $t\geq 0$. Here $\rho_{i}=\rho_{i}(x,t)$ is the
charge carrier density for the $i$-th species, and $\phi=\phi(x,t)$ the
electrostatic potential. $D_{i}(x)$ is the diffusion coefficient, $z_{i}$ is
the rescaled charge. In the Poisson equation, $\epsilon(x)$ is the
permittivity, $f(x)$ is the permanent (fixed) charge density of the system,
$s$ is the number of species.

Due to the wide variety of devices modeled by the PNP equations, computer
simulation for this system of differential equations is of great interest.
However, the PNP system is a strongly coupled system of nonlinear equations,
also, the PNP system as a gradient flow can take very long time evolution to
reach steady states. Hence, designing efficient and stable numerical methods
for the PNP system remains an active area of research (see, e.g., [22, 28, 25,
35, 36, 8]).

PNP system possesses two immediate properties: it preserves the non-negativity
of $\rho_{i}$ and conserves total mass. Therefore, we can consider non-
negative initial data with mass one, so that the density is in the set of
probability measures $\mathcal{P}(\Omega)$ on $\Omega$. The third property is
the dissipation of the total energy, which can be expressed as follows. Given
energy

$E=\int_{\Omega}\bigg{(}\sum_{i=1}^{s}\rho_{i}\log\rho_{i}+\frac{1}{2}(f+\sum_{i=1}^{s}z_{i}\rho_{i})\phi\bigg{)}dx+B,$
(1.2)

with boundary correction term $B$, the NP equation (1.1a) can be written as
the gradient flow

$\partial_{t}\rho_{i}=\nabla\cdot(\rho_{i}D_{i}(x)\nabla(\delta_{\rho_{i}}E)).$

Differentiating the energy along solutions of the PNP system, one formally
obtains the energy dissipation along the gradient flow

$\frac{dE}{dt}=-\int_{\Omega}\sum_{i=1}^{s}D_{i}(x)\rho_{i}|\nabla(\log\rho_{i}+z_{i}\phi)|^{2}dx\leq
0,$

which indicates that the solution evolves in the direction of steepest descent
of the energy. This property entails a characterization of the set of
stationary states, and provides a useful tool to study its stability.
Numerical methods for (1.1) are desired to attain all three properties at the
discrete level, which are rather challenging.

###  1.1. Related work

The most common numerical approach is the direct discretization of (1.1) using
classical finite difference, finite volume, finite element, or discontinuous
Galerkin methods [10, 23, 13, 31, 11, 24, 12, 21, 22, 14, 7, 15, 28, 25, 35,
36, 8]. Such methods are explicit or semi-implicit in time, so the per time
computation is cheap. But it is often challenging to ensure both unconditional
positivity and discrete energy decay simultaneously. The nonlinearity also
complicates the way to obtain solutions when applying implicit or semi-
implicit solvers.

For the gradient flow

$\partial_{t}\rho=-\nabla^{W}_{\rho}E(\rho),$

with the gradient with respect to the quadratic Wasserstein metric
$W_{2}(\cdot,\cdot)$, the minimizing movement approximation (see [1] and the
references therein)

$\rho^{n+1}=\text{argmin}_{\rho\in
K}\left\\{\frac{1}{2\tau}W^{2}_{2}(\rho^{n},\rho)+E(\rho)\right\\},\quad\rho^{0}=\rho^{in}(x),$
(1.3)

also named Jordan-Kinderlehrer-Otto (JKO) scheme (Jordan et al. [19]), defines
a sequence $\\{\rho^{n}\\}$ in the probability space $K$ to approximate the
solution $\rho(x,n\tau)$, where $\tau>0$ is the time step. Since $\rho^{n}$ is
in the probability space, thus method (1.3) is positivity and mass preserving.
The fact $W_{2}(\rho^{n},\rho)\geq 0$ ensures the energy dissipation for any
$\tau>0$. We refer the interested reader to [26, 5, 30] for some JKO-type
schemes for gradient flows in the probability space.

The authors in [20] constructed the JKO scheme for a two species PNP system
with constant coefficient $D_{i}(x)=1$ and $\epsilon(x)=1$. The existence of
the unique minimizer to the JKO scheme and convergence of the minimizer to the
weak solution of the PNP system have been established in [20]. However, the
fully discrete JKO-type scheme for the multi-species variable coefficient PNP
system on a bounded domain has not been studied yet. This is what we aim to
accomplish in this paper.

###  1.2. Our contributions

In the present work, we will extend Benamou-Brenier’s dynamic formulation [3]
for the Wasserstein distance to the present setting that resolve the
aforementioned issues.

  * •

In a general setting with mixed types of boundary conditions, we identify the
total energy functional which is dissipating along the solution trajectories,
we also establish (Theorem 2.2) lower energy bounds with coercivity for such a
functional.

  * •

We construct a Wasserstein-type distance and formulate a corresponding
variational scheme. The update at each step reduces to solving a constrained
minimization problem, for which we prove unique solvability (Theorem 2.3).
Three solution properties: mass conservation, positivity, and energy
dissipation are shown to be preserved in time (Theorem 2.4).

  * •

We further convert the variational scheme into a dynamic formulation, which
for variable diffusion coefficients extends the classical Bennamou-Breiner
formulation. To reduce computational cost, we use a local approximation for
the artificial time in the constraint transport equation by a one step
difference and the integral in time by a one term quadrature. The resulting
minimization problem is shown to be a first order time consistent scheme for
the PNP system (Theorem 2.5).

  * •

We present a fully discrete scheme, and prove its unique solvability (Theorem
3.1).

  * •

We establish that for any fixed time step and spatial meth size, density
positivity will be propagating over all time steps (Theorem 3.2). This is in
sharp contrast to the work [26], where Fisher information regularization is
added to enforce solution positivity for the aggregation equation.

  * •

The fully-discrete minimization problem reduces to a convex optimization
problem with linear constraints, and can be solved by some efficient
optimization solvers. Our numerical tests are conducted with a simple
projected gradient algorithm.

  * •

Numerical results are provided to demonstrate the superior performance of the
proposed method.

###  1.3. Organization

The paper is organized as follows. In the next section, we provide necessary
background on the dynamical formulation of the PNP system, main solution
properties and its relation with Wasserstein gradient flows. We then derive
the semi-discrete scheme. In Section 3, we introduce a fully discrete scheme
and study the properties of this scheme. Numerical algorithms are given in
Section 4. Numerical results are provided in Section 5, and the paper is
concluded in Section 6.

Notation. We use $[n]$ to denote $\\{1,2,\cdots,n\\}$ for any integer $n$. For
vector $\phi$ (fully discrete case), its amplitude is denoted by $|\phi|$. For
function $\phi$, $\|\phi\|$ is its $L^{2}$ norm.

##  2\. Model background and semi-discretization

In this section we briefly review the model setup and the corresponding
Wasserstein gradient flow.

###  2.1. Boundary conditions

Boundary conditions are a critical component of the PNP model and determine
important qualitative behaviors of the solution. Let $\Omega$ be a bounded
domain with Lipschitz boundary $\partial\Omega$. We use the no-flux boundary
condition for the NP equations, i.e.,

$D_{i}(x)\left(\nabla\rho_{i}+z_{i}\rho_{i}\nabla\phi\right)\cdot\mathbf{n}=0,\quad
x\in\partial\Omega,\quad i=1,\cdots,s.$ (2.1)

Here, $\mathbf{n}$ is the outer unit normal at the boundary point
$x\in\partial\Omega.$

The external electrostatic potential $\phi$ is influenced by applied
potential, which can be modeled by prescribing a boundary condition. Here we
consider a general form of boundary conditions:

$\alpha\phi+\beta\epsilon(x)\frac{\partial\phi}{\partial\mathbf{n}}=\phi^{b},\quad
x\in\partial\Omega.$ (2.2)

Here $\alpha,\beta$ are physical parameters such that $\alpha\cdot\beta\geq
0$, and $\phi^{b}=\phi^{b}(x,t)$ is a given function. With such setup, we are
to solve the following initial-boundary value problem:

$\left\\{\begin{array}[]{lll}&\partial_{t}\rho_{i}=\nabla\cdot\left[D_{i}(x)\left(\nabla\rho_{i}+z_{i}\rho_{i}\nabla\phi\right)\right],&x\in\Omega,\quad
t>0,\ i=1,\cdots,s,\\\
&-\nabla\cdot(\epsilon(x)\nabla\phi)=f(x)+\sum_{i=1}^{s}z_{i}\rho_{i},&x\in\Omega,\quad
t>0,\\\ &\rho_{i}(x,0)=\rho_{i}^{in}(x)\geq 0,&x\in\Omega,\quad
i=1,\cdots,s,\\\
&\left(D_{i}(x)\left(\nabla\rho_{i}+z_{i}\rho_{i}\nabla\phi\right)\right)\cdot\mathbf{n}=0,&x\in\partial\Omega,\quad
t>0,\ i=1,\cdots,s,\\\
&\alpha\phi+\beta\epsilon(x)\frac{\partial\phi}{\partial\mathbf{n}}=\phi^{b},&x\in\partial\Omega,t>0.\end{array}\right.$
(2.3)

######  Remark 2.1.

(2.2) includes three typical forms: (i) the Robin boundary condition
($\alpha=1,\beta>0$) models a capacitor [10], (ii) the Dirichlet boundary
condition ($\alpha=1,\ \beta=0$) models an applied voltage, and (iii) the
Neumann boundary condition ($\alpha=0,\ \beta=1$) models surface changes. The
case of pure Neumann boundary conditions requires the compatibility condition

$\int_{\Omega}\left(f(x)+\sum_{i=1}^{s}z_{i}\rho^{in}_{i}\right)dx+\int_{\partial\Omega}\epsilon(x)\phi^{b}ds=0,$
(2.4)

and an additional constraint such as $\int_{\Omega}\phi(x,t)dx=0$ so that
$\phi$ is uniquely defined.

Any combination of these three types can be applied to $\phi$ on a disjoint
partition of the boundary. In what follows, we set

$\partial\Omega=\Gamma_{D}\cup\Gamma_{N}\cup\Gamma_{R},$

and on each part, one type of boundary condition is imposed, i.e.,

$\alpha=\begin{dcases}1,&\text{ on }\Gamma_{D},\\\ 0,&\text{ on
}\Gamma_{N},\\\ 1,&\text{ on }\Gamma_{R},\\\
\end{dcases}\quad\beta=\begin{dcases}0,&\text{ on }\Gamma_{D},\\\ 1,&\text{ on
}\Gamma_{N},\\\ \beta_{R},&\text{ on }\Gamma_{R},\\\
\end{dcases}\quad\phi^{b}=\begin{dcases}\phi^{b}_{D},&\text{ on
}\Gamma_{D},\\\ \phi^{b}_{N},&\text{ on }\Gamma_{N},\\\ \phi^{b}_{R},&\text{
on }\Gamma_{R}.\\\ \end{dcases}$

The existence and uniqueness of the solution for the nonlinear PNP boundary
value problems with different boundary conditions have been studied in [18,
27, 33] for the 1D case and in [4, 17] for multi-dimensions.

###  2.2. Energy functional: dissipation and coercivity

In the presence of homogeneous boundary conditions on $\phi$, i.e.,
$\phi^{b}=0$, the PNP system is energetically closed in the sense that the
free energy functional associated to (1.1) is of form

$E_{0}=\int_{\Omega}\bigg{(}\sum_{i=1}^{s}\rho_{i}\log\rho_{i}+\frac{1}{2}(f+\sum_{i=1}^{s}z_{i}\rho_{i})\phi\bigg{)}dx,$
(2.5)

which along solution trajectories is dissipating in time. For general boundary
conditions with $\phi^{b}\not=0$, we need to modify the energy so that it is
still dissipating along the solution of the PNP system. To this end, we
differentiate (2.5) along the solution of (1.1), with integration by parts
using (2.1), we have

$\displaystyle\frac{d}{dt}E_{0}(\rho,\phi)(t)$
$\displaystyle=-\int_{\Omega}\sum_{i=1}^{s}D_{i}(x)\rho_{i}|\nabla(\log\rho_{i}+z_{i}\phi)|^{2}dx+\frac{1}{2}\int_{\partial\Omega}\epsilon(x)\left[\phi(\partial_{n}\phi)_{t}-(\partial_{n}\phi)\phi_{t}\right]ds.$

Assume that $\phi^{b}$ does not depend on time, then
$\alpha\phi_{t}+\beta\epsilon(x)\partial_{n}\phi_{t}=0$ on $\partial\Omega$,
this allows us to express the last term as

$\frac{1}{2}\frac{d}{dt}\left[\int_{\Gamma_{D}}\epsilon(x)\phi^{b}_{D}\partial_{n}\phi
ds-\int_{\Gamma_{N}}\phi^{b}_{N}\phi
ds-\frac{1}{\beta_{R}}\int_{\Gamma_{R}}\phi^{b}_{R}\phi ds\right].$

Thus the modified total energy functional can be taken as

$\displaystyle
E=E_{0}-\frac{1}{2}\left[\int_{\Gamma_{D}}\epsilon(x)\phi^{b}_{D}\partial_{n}\phi
ds-\int_{\Gamma_{N}}\phi^{b}_{N}\phi
ds-\frac{1}{\beta_{R}}\int_{\Gamma_{R}}\phi^{b}_{R}\phi ds\right].$ (2.6)

Using the Poisson equation, the total energy can be rewritten as

$E(\rho,\phi)=\int_{\Omega}\bigg{(}\sum_{i=1}^{s}\rho_{i}\log\rho_{i}+\frac{1}{2}\epsilon(x)|\nabla\phi|^{2}\bigg{)}dx-\int_{\Gamma_{D}}\epsilon(x)\phi_{D}^{b}\partial_{n}\phi
ds+\frac{1}{2\beta_{R}}\int_{\Gamma_{R}}|\phi|^{2}ds.$ (2.7)

######  Proposition 2.1.

Assume that $\phi^{b}$ does not depend on time, then the extended energy
functional (2.7) satisfies

$\frac{d}{dt}E(\rho,\phi)(t)=-\int_{\Omega}\sum_{i=1}^{s}D_{i}(x)\rho_{i}|\nabla(\log\rho_{i}+z_{i}\phi)|^{2}dx\leq
0,\quad t>0,$ (2.8)

along the solution of (1.1).

Recall that on $\Gamma_{D}$, the usual strategy for analysis is to transform
it to the case with zero boundary value for $\phi$. This way the modified
energy would include an additional term called the external potential energy.
For simplicity, we take $\phi_{D}^{b}=0$, so that we have the following
result.

######  Theorem 2.2.

(Lower bound and coercivity of $E$ ) Let $\Omega$ be an open, bounded
Lipschitz domain, and $\phi^{b}$ be independent of time with $\phi_{D}^{b}=0$,
$\beta_{R}>0$, and $\epsilon(x)\geq a>0$. Then the energy of form

$E(\rho,\phi)=\int_{\Omega}\bigg{(}\sum_{i=1}^{s}\rho_{i}\log\rho_{i}+\frac{1}{2}\epsilon(x)|\nabla\phi|^{2}\bigg{)}dx+\frac{1}{2\beta_{R}}\int_{\Gamma_{R}}\phi^{2}ds$
(2.9)

is bounded from below. Moreover, there exist constants $c_{0},c_{1}>0$ such
that

$E(\rho,\phi)\geq c_{0}\|\phi\|_{H^{1}}^{2}-c_{1}.$ (2.10)

###### Proof.

For $\rho_{i}\geq 0$, we have
$\int_{\Omega}\sum_{i=1}^{s}\rho_{i}\log(\rho_{i})\geq-s|\Omega|/e=:-c_{1}$.
For the $\phi$-dependent part in $E$, we argue for all possible cases. For
$\Gamma_{D}\not=\emptyset$ we have $\phi_{D}^{b}=0$; for purely Neumann’s
condition we have the additional condition $\int_{\Omega}\phi(x)dx=0$, in
either case we can apply the Poincaré inequality or the Poincaré–Wirtinger
inequality to conclude

$\|\phi\|^{2}_{L^{2}}\leq c^{*}\|\nabla\phi\|^{2}_{L^{2}}$

with constant $c^{*}$ depending on the geometry of $\Omega$, hence

$E\geq-c_{1}+\frac{a}{2}\|\nabla\phi\|^{2}_{L^{2}}\geq
c_{0}\|\phi\|^{2}_{H^{1}}-c_{1},\quad
c_{0}=\frac{a}{4}\min\\{1,\frac{1}{c^{*}}\\}.$

For the case $\partial\Omega=\Gamma_{R}\cup\Gamma_{N}$ with
$\Gamma_{R}\not=\emptyset$, we have

$E(\rho,\phi)\geq\frac{1}{2}\min\\{a,\beta_{R}^{-1}\\}\tilde{E}-c_{1}$

with
$\tilde{E}(\phi):=\int_{\Omega}|\nabla\phi|^{2}dx+\int_{\partial\Omega}|\phi|^{2}ds$.
We claim that

$\tilde{E}(\phi)\geq c\|\phi\|_{H^{1}}^{2}\quad\text{for some}\quad c>0,$

which can be proved with a contradiction argument. Since otherwise we can
assume $\tilde{E}(\phi_{n})<\frac{1}{n}\|\phi_{n}\|^{2}_{H^{1}}$. Set
$w_{n}=\phi_{n}/\|\phi_{n}\|_{H^{1}}$, then $w_{n}\in H^{1}(\Omega)$ with

$\|w_{n}\|_{H^{1}}=1\;\text{and}\;\|\nabla w_{n}\|^{2}_{L^{2}}<1/n.$

By the Rellich-Kondrachov theorem, we can extract a subsequence
$\\{w_{n_{k}}\\}$ weakly converging to $w$ in $H^{1}(\Omega)$ with $\nabla
w_{n_{k}}\to 0$ weakly in $L^{2}(\Omega)$. This allows us to conclude $w\in
H^{1}$, and $\nabla w=0$. From $\int_{\Gamma_{R}}|w_{n_{k}}|^{2}ds<1/n_{k}$
and

$\|w\|_{L^{2}(\Gamma_{R})}\leq\|w_{n_{k}}\|_{L^{2}(\Gamma_{R})}+\|w_{n_{k}}-w\|_{L^{2}(\Gamma_{R})}\leq
1/\sqrt{n_{k}}+C\|w-w_{n_{k}}\|_{H^{1}},$

we obtain $w=0$ on $\Gamma_{R}$. Hence $w=0$ a.e., this is a contradiction. We
complete this case by taking $c_{0}=\frac{c}{2}\min\\{a,\beta_{R}^{-1}\\}$. ∎

###  2.3. Wasserstein distance and JKO scheme for multi-density

In order to derive a variational scheme for the PNP system with multi-density,
we need to introduce a Wasserstein-type distance. Motivated by the well-known
characterization of the Wasserstein distance in a one-component fluid obtained
by Benamou-Brenier [3], we consider to minimize a joint functional over the
set

$\displaystyle K:$
$\displaystyle=\\{\rho=(\rho_{1},\cdots,\rho_{s}),u=(u_{1},\cdots,u_{s}):$
(2.11)
$\displaystyle\quad\partial_{t}\rho_{i}+\nabla\cdot(\rho_{i}u_{i})=0,\quad(\rho_{i}u_{i})\cdot\mathbf{n}=0\quad\text{
on }\partial\Omega\times[0,1],$
$\displaystyle\quad\rho_{i}\in\mathcal{P}(\Omega),\quad\rho_{i}(x,0)=\rho_{i}^{0}(x),\quad\rho_{i}(x,1)=\rho_{i}^{1}(x)\\}.$

For the PNP system of two species $s=2$ with $D_{i}(x)=1$ and $\epsilon(x)=1$
considered in [20], the distance inherited from the 2-Wasserstein distance is
defined by

$d^{2}(\rho^{0},\rho^{1})=\sum_{i=1}^{2}W_{2}^{2}(\rho_{i}^{0},\rho_{i}^{1}).$

This is equivalent to the minimization of the joint functional:

$d^{2}(\rho^{0},\rho^{1}):=\min_{(\rho,u)\in
K}\sum_{i=1}^{2}\int_{0}^{1}\int_{\Omega}|u_{i}|^{2}\rho_{i}dxdt.$ (2.12)

Here $t$ is an artificial time and serves to characterize the optimal curve in
the density space. Following [19], the authors in [20] constructed the
following JKO scheme: Given a time step $\tau$, the scheme defines a sequence
$\rho^{n}$ as

$\rho^{0}=\rho^{in},\quad\rho^{n+1}=\arg\min_{\rho\in[\mathcal{P}(\Omega)]^{2}}\left\\{\frac{1}{2\tau}d^{2}(\rho^{n},\rho)+\mathcal{E}(\rho)\right\\}.$
(2.13)

Here $\mathcal{E}$ is the total free energy, $d^{2}$ is the (squared) distance
on the product space as defined in (2.12). One of the challenges in this
program lies in handling the coupling terms, some intrinsic difficulties arise
due to both the specific Poisson kernel and the system setting. Note that in
[20] with $\epsilon(x)=1$, the electrostatic potential $\phi$ in
$E(\rho,\phi)$ is replaced by

$\phi[\rho]=N*(f+\sum_{i=1}^{2}z_{i}\rho_{i}),\quad x\in\Omega,$

so that $\mathcal{E}(\rho)=E(\rho,\phi[\rho])$. Here the kernel $N\sim
C/|x|^{d-2}$ serves as a counterpart of the Green’s function for the Newton
potential in $\mathbb{R}^{d}$. Even with this treatment derivation of the
corresponding Euler–Lagrange equations is quite delicate. We refer to [20] for
further details.

In order to extend the above JKO-type scheme to the present setting, we face
two new difficulties: (i) $D_{i}(x)$ is no longer a constant, the kinetic
energy corresponding to the squared distance cost needs to be modified; (ii)
$\epsilon(x)$ is a general non-negative function, $\phi$ cannot be expressed
explicitly in terms of $\rho$. As for (i), we follow [16] and consider a
modified functional

$d^{2}(\rho^{0},\rho^{1}):=\min_{(\rho,u)\in
K}\sum_{i=1}^{s}\int_{0}^{1}\int_{\Omega}D_{i}^{-1}|u_{i}|^{2}\rho_{i}dxdt.$
(2.14)

As for (ii), the Poisson equation is treated as a constraint in the resulting
minimization problem. For ease of presentation we define

$\mathcal{A}:=\left\\{(\rho,\phi):-\nabla\cdot(\epsilon(x)\nabla\phi)=f(x)+\sum_{i=1}^{s}z_{i}\rho_{i},\quad\alpha\phi+\beta\frac{\partial\phi}{\partial\mathbf{n}}=\phi^{b},x\in\partial\Omega,\rho\in[\mathcal{P}(\Omega)]^{s}\right\\}.$
(2.15)

For fixed $\rho^{*}\in[\mathcal{P}(\Omega)]^{s}$, and time step $\tau>0$ we
set

$G_{\tau}(\rho,\phi)=\frac{1}{2\tau}d^{2}(\rho^{*},\rho)+E(\rho,\phi),\quad(\rho,\phi)\in\mathcal{A}.$
(2.16)

In order to define a discrete sequence of approximate solutions using the
minimizing scheme, we present a result on the existence of minimizers of
$G_{\tau}$. To establish the uniqueness, we now prepare a technical lemma,
with (iii) to be used later in the proof of Theorem 3.1.

######  Lemma 2.1.

Given $X^{0},X^{1}$, let $X(\theta)=\theta X^{0}+(1-\theta)X^{1}$ for any
$\theta\in(0,1)$.
(i) If $X^{0},X^{1}$ are vectors, then

$|X(\theta)|^{2}-\theta|X^{0}|^{2}-(1-\theta)|X^{1}|^{2}=-\theta(1-\theta)|X^{1}-X^{0}|^{2}.$
(2.17)

(ii) If $X^{0}>0,X^{1}>0$ are scalars, then

$X(\theta)\log X(\theta)-\theta X^{0}\log X^{0}-(1-\theta)X^{1}\log
X^{1}=-\theta(1-\theta)(X^{1}-X^{0})^{2}g(X^{0},X^{1};\theta),$ (2.18)

for some positive function $g$ depending on $X^{0},X^{1}$ and $\theta$.
(iii) If $X^{0}>0,X^{1}>0$, $Y^{0},Y^{1}$ are scalars, then

$\frac{(Y(\theta))^{2}}{X(\theta)}-\theta\frac{(Y^{0})^{2}}{X^{0}}-(1-\theta)\frac{(Y^{1})^{2}}{X^{1}}=-\theta(1-\theta)\frac{(X^{1}Y^{0}-X^{0}Y^{1})^{2}}{X^{0}X^{1}X(\theta)}.$
(2.19)

###### Proof.

We only prove (ii); for (i) and (iii) can be verified by a direct calculation.
Note that

$X(\theta)\log X(\theta)=\theta
X^{0}\log(X(\theta))+(1-\theta)X^{1}\log(X(\theta)).$ (2.20)

Taylor’s expansion of $\log(X(\theta)$ at $X^{0}$ and $X^{1}$, respectively,
gives

$\log(\theta
X^{0}+(1-\theta)X^{1})=\log(X^{0})+\frac{1}{X^{0}}(1-\theta)(X^{0}-X^{1})-\frac{(1-\theta)^{2}(X^{1}-X^{0})^{2}}{(\tilde{X}^{0})^{2}},$

where $\tilde{X}^{0}$ in between $X^{0}$ and $X(\theta)$, and

$\log(\theta
X^{0}+(1-\theta)X^{1})=\log(X^{1})+\frac{1}{X^{1}}+\theta(X^{1}-X^{0})-\frac{\theta^{2}(X^{1}-X^{0})^{2}}{(\tilde{X}^{1})^{2}},$

where $\tilde{X}^{1}$ in between $X^{1}$ and $X(\theta)$. Substituting these
into the right hand side of (2.20) leads to

$\displaystyle X(\theta)\log X(\theta)=$ $\displaystyle\theta X^{0}\log
X^{0}-\theta(1-\theta)^{2}\frac{X^{0}}{(\tilde{X}^{0})^{2}}(X^{1}-X^{0})^{2}$
$\displaystyle+(1-\theta)X^{1}\log
X^{1}-\theta^{2}(1-\theta)\frac{X^{1}}{(\tilde{X}^{1})^{2}}(X^{1}-X^{0})^{2},$

this completes the proof of (ii) by defining
$g(X^{0},X^{1},\theta)=\frac{(1-\theta)X^{0}}{(\tilde{X}^{0})^{2}}+\frac{\theta
X^{1}}{(\tilde{X}^{1})^{2}}>0$. ∎

######  Theorem 2.3.

(Existence of minimizers) Fix $\tau>0$, and
$\rho^{*}\in[\mathcal{P}(\Omega)]^{s}$. Then the functional
$G_{\tau}(\rho,\phi)$ admits a unique minimizer on $\mathcal{A}$.

###### Proof.

By Theorem 2.2, $G_{\tau}$ is bounded from below on $\mathcal{A}$, hence there
is a minimizing sequence $(\rho^{k},\phi^{k})$ and $\rho^{k}$ is tight and
uniformly integrable. By the Dunford–Pettis Theorem one may extract a
subsequence such that $\rho^{k}\to\rho$ in $L^{1}(\Omega)$, which together
with $\rho^{k}\in[\mathcal{P}(\Omega)]^{s}$ ensure that
$\rho\in[\mathcal{P}(\Omega)]^{s}$. In addition, $E(\rho,\cdot)$ is also
coercive in $\phi$ because of (2.10), i.e.,

$E(\rho,\phi)\geq c_{0}\|\phi\|^{2}_{1}-c_{1}.$

Hence one may extract a subsequence such that $\phi^{k}\to\phi$ weakly in
$H^{1}(\Omega)$. The weak $L^{1}$ lower semi-continuity (l.s.c.) of the
squared Wasserstein distance can be easily adapted to the present case. The
lower semicontinuity of $E$ with respect to weak convergence can be seen from
the following inequality

$E(\rho^{k},\phi^{k})\geq E(\rho,\phi)+\int_{\Omega}\left[\sum_{i=1}^{s}{\rm
ln}\rho_{i}(\rho_{i}^{k}-\rho_{i})+\epsilon(x)\nabla\phi\cdot(\nabla\phi^{k}-\nabla\phi)\right]dx+\frac{\alpha}{\beta}\int_{\partial\Omega}\phi(\phi^{k}-\phi)ds.$

Putting all these together we claim that the limit is a minimizer.

Finally, the uniqueness comes from the fact that the admissible set
$\mathcal{A}$ is convex w.r.t. linear interpolation and that the total free
energy is jointly strictly convex in $(\rho,\phi)$ on $\mathcal{A}$. More
precisely, we argue as follows. Let $\theta\in(0,1)$, then
$\rho(\theta)=\theta\rho^{0}+(1-\theta)\rho^{1}$ is a convex linear
combination for $\rho^{0}$ and $\rho^{1}$. Let $\phi^{0}$ and $\phi^{1}$ be
obtained from the Possion equation, corresponding to $\rho^{0}$ and
$\rho^{1}$, respectively. Then
$\phi(\theta)=\theta\phi^{0}+(1-\theta)\phi^{1}$ must be the solution to the
Poisson equation corresponding to $\rho(\theta)$. For the energy of form
(2.9), we evaluate $E(\rho(\theta),\phi(\theta))$ term by term to determine
whether it is strictly convex. Using (2.18) for $\rho^{l}_{i}=X^{l}$ and
(2.17) for $X^{l}=\nabla\phi^{l}$ in $\Omega$, and (2.17) for $X^{l}=\phi^{l}$
on $\Gamma_{R}$, respectively, we obtain

$\displaystyle E(\rho(\theta),\phi(\theta))-\theta
E(\rho^{0},\phi^{0})-(1-\theta)E(\rho^{1},\phi^{1})=-\theta(1-\theta)I$

with

$I=\int_{\Omega}\bigg{(}\sum_{i=1}^{s}(\rho_{i}^{1}-\rho_{i}^{0})^{2}g(\rho_{i}^{0},\rho^{1}_{i};\theta)+\frac{1}{2}\epsilon(x)|\nabla(\phi^{1}-\phi^{0})|^{2}\bigg{)}dx+\frac{1}{2\beta_{R}}\int_{\Gamma_{R}}(\phi^{1}-\phi^{0})^{2}ds.$

Convexity of $E$ follows from $I\geq 0$. Actually this inequality is strict,
unless $\rho^{0}=\rho^{1},\phi^{0}=\phi^{1}$, which can be derived from
letting $I=0$. Hence $E(\rho,\phi)$ is strictly convex under two linear
constraints. ∎

We are now ready to present a variational scheme formulation – a JKO-type
scheme for (2.3): given time step $\tau>0$, recursively we define a sequence
$\\{\rho^{n},\phi^{n}\\}$ by

$\displaystyle\rho^{0}=\rho^{\rm in},\quad(\rho^{n+1},\phi^{n+1})=$
$\displaystyle\arg\min_{(\rho,\phi)\in\mathcal{A}}\left\\{\frac{1}{2\tau}d^{2}(\rho^{n},\rho)+E(\rho,\phi)\right\\}.$
(2.21)

######  Theorem 2.4.

(Solution properties of scheme (2.21) )
(i) (Probability-preserving) If $\rho^{n}\in[\mathcal{P}(\Omega)]^{s}$, so is
$\rho^{n+1}$;
(ii) (Unconditionally energy stability) the inequality

$E(\rho^{n+1},\phi^{n+1})+\frac{1}{2\tau}d^{2}(\rho^{n},\rho^{n+1})\leq
E(\rho^{n},\phi^{n})$

holds for any $n\geq 0$. Furthermore,

$\sum_{n=0}^{\infty}d^{2}(\rho^{n},\rho^{n+1})\leq
2\tau(E(\rho^{0},\phi^{0})-\inf_{(\rho,\phi)\in\mathcal{A}}E(\rho,\phi)).$
(2.22)

###### Proof.

(i) The constraint $\mathcal{A}$ ensures that
$\rho^{n}\in[\mathcal{P}(\Omega)]^{s}$ which is inherited from initial data;
namely the method is both positivity and mass preserving.
(ii) From the definition of the minimizer, it follows

$E(\rho^{n+1},\phi^{n+1})+\frac{1}{2\tau}d^{2}(\rho^{n},\rho^{n+1})\leq
E(\rho^{n},\phi^{n}).$

Here we used $d^{2}(\rho,\rho)=0$ for any $\rho\in[\mathcal{P}(\Omega)]^{s}$.
Finally, summation over $n$ yields (2.22). ∎

###  2.4. Semi-discrete JKO scheme

We proceed to obtain a computable formulation. Let $m_{i}=\rho_{i}u_{i}$, the
dynamic formulation of the distance $d^{2}(\cdot,\cdot)$ in (2.21) can be
expressed as: given $\rho^{n}(x)$, we have

$\displaystyle(\rho^{n+1},\phi^{n+1})=$
$\displaystyle\arg\min_{(\rho,\phi)\in\mathcal{A},m}\left\\{\frac{1}{2\tau}\sum_{i=1}^{s}\int_{0}^{1}\int_{\Omega}F(\rho_{i},m_{i})D^{-1}_{i}dxdt+E(\rho(\cdot,1),\phi(\cdot,1))\right\\},$
(2.23) $\displaystyle\text{ s.t.
}\quad\partial_{t}\rho_{i}+\nabla\cdot(m_{i})=0,\quad
m_{i}\cdot\mathbf{n}=0,\quad x\in\partial\Omega,\;\rho(x,0)=\rho^{n}.$

Here $t$ is an artificial time, and

$F(\rho_{i},m_{i})=\left\\{\begin{array}[]{ll}\frac{|m_{i}|^{2}}{\rho_{i}}&\text{if}\;\rho_{i}>0,\\\
0&\text{if}\;(\rho_{i},m_{i})=(0,0),\\\
+\infty&\text{otherwise}.\end{array}\right.$

The use of $m_{i}$ has enhanced the functional convexity in $m_{i}$ and made
the transport constraint linear (see Breiner [3]), yet causing difficulties
for solutions near $\rho_{i}=0$. We shall prove for the fully discrete case
positivity of $\rho_{i}^{n}$ is preserved for all $n$. Another computational
overhead with (2.23) is dealing with the artificial time $t\in[0,1]$ which is
induced by the optimal transport flow. To overcome this issue, we follow [26]
with a local approximation in the artificial time : approximate the derivative
in $t$ in the constraint transport equation by a one step difference and the
integral in time in the objective function by a one term quadrature. We thus
obtain the following scheme:

$\displaystyle(\rho^{n+1},\phi^{n+1})=$
$\displaystyle\arg\min_{(\rho,\phi)\in\mathcal{A},m}\left\\{\frac{1}{2\tau}\sum_{i=1}^{s}\int_{\Omega}F(\rho_{i},m_{i})D^{-1}_{i}dx+E(\rho,\phi)\right\\},$
(2.24) $\displaystyle\text{ s.t.
}\quad\rho_{i}-\rho_{i}^{n}+\nabla\cdot(m_{i})=0,\quad
m_{i}\cdot\mathbf{n}=0,\quad x\in\partial\Omega.$

######  Theorem 2.5.

The positive minimizer of the variational problem (2.24) is a first-order time
consistent scheme for the PNP system.

###### Proof.

Let (2.24) admit a minimizer with $\rho>0$. We can derive optimal conditions
by the Lagrange multiplier method. Define the Lagrangian as

$\displaystyle L(\rho,\phi,m,v,\xi):=$
$\displaystyle\frac{1}{2\tau}\sum_{i=1}^{s}\int_{\Omega}F(\rho_{i},m_{i})D^{-1}_{i}dx+E(\rho,\phi)+\int_{\partial\Omega}\xi(\alpha\phi+\beta\partial_{n}\phi-\phi^{b})ds$
$\displaystyle+\sum_{i=1}^{s}\int_{\Omega}v_{i}(\rho_{i}-\rho_{i}^{n}+\nabla\cdot
m_{i})dx+\int_{\Omega}v_{s+1}(f+\sum_{i=1}^{s}z_{i}\rho_{i}+\nabla\cdot(\epsilon(x)\nabla\phi))dx.$

The optimality conditions for $x\in\Omega$ are

$\displaystyle\frac{\delta L}{\delta\rho_{i}}=0\quad$ implies
$\displaystyle-\frac{1}{2\tau}\frac{||m_{i}||^{2}}{\rho_{i}^{2}}D^{-1}_{i}+\log(\rho_{i})+1+\frac{1}{2}z_{i}\phi+v_{i}+z_{i}v_{s+1}=0,$
$\displaystyle\quad i=1,\cdots,s,$ $\displaystyle\frac{\delta
L}{\delta\phi}=0\quad$ implies
$\displaystyle\frac{1}{2}(f+\sum_{i=1}^{s}z_{i}\rho_{i})+\nabla\cdot(\epsilon(x)\nabla
v_{s+1})=0,$ $\displaystyle\frac{\delta L}{\delta m_{i}}=0\quad$ implies
$\displaystyle\frac{1}{\tau}\frac{m_{i}}{\rho_{i}}D_{i}^{-1}-\nabla\cdot
v_{i}=0,$ $\displaystyle\quad i=1,\cdots,s,$ $\displaystyle\frac{\delta
L}{\delta v_{i}}=0\quad$ implies
$\displaystyle\rho_{i}-\rho_{i}^{n}+\nabla\cdot m_{i}=0,$ $\displaystyle\quad
i=1,\cdots,s,$ $\displaystyle\frac{\delta L}{\delta v_{s+1}}=0\quad$ implies
$\displaystyle
f+\sum_{i=1}^{s}z_{i}\rho_{i}+\nabla\cdot(\epsilon(x)\nabla\phi)=0.$

For $x\in\Omega$, we thus have

$v_{i}=\frac{1}{2\tau}\frac{||m_{i}||^{2}}{\rho_{i}^{2}}D^{-1}_{i}-\log(\rho_{i})-1-\frac{1}{2}z_{i}\phi-
z_{i}v_{s+1},\quad m_{i}=\tau D_{i}\rho_{i}\nabla v_{i}$

and

$\nabla\cdot(\epsilon(x)\nabla
v_{s+1})=\frac{1}{2}\nabla\cdot(\epsilon(x)\nabla\phi).$

On $\partial\Omega$, from integrating by parts in calculating $\delta L$ there
remain the following boundary terms

$\int_{\partial\Omega}\epsilon(x)\delta(\partial_{n}\phi)v_{s+1}ds-\int_{\partial\Omega}\epsilon(x)\delta\phi\partial_{n}v_{s+1}ds+\int_{\partial\Omega}v_{i}\delta
m_{i}\cdot\mathbf{n}ds,$

where the last term vanishes due to the constraint $m_{i}\cdot\mathbf{n}=0.$
In addition, we need also consider terms arising from

$\delta
B+\delta\int_{\partial\Omega}\xi(\alpha\phi+\beta\partial_{n}\phi-\phi^{b})ds.$

Upon careful regrouping, we have two cases to distinguish:
(i) for $\beta\not=0$, the correction term $B$ in the energy (1.2) is given by

$B=\frac{1}{2\beta}\int_{\partial\Omega}\phi^{b}\phi ds.$

We obtain

$\epsilon(x)v_{s+1}+\beta\epsilon(x)\xi=0,\quad-\epsilon(x)\partial_{n}v_{s+1}+\alpha\xi+\frac{1}{2\beta}\phi^{b}=0,\quad\text{on}\;\partial\Omega;$

(ii) For $\beta=0$,
The correction term $B$ in the energy (1.2) is given by

$B=-\frac{1}{2\alpha}\int_{\partial\Omega}\epsilon(x)\phi^{b}\partial_{n}\phi
ds.$

from which we have

$\epsilon(x)v_{s+1}-\frac{1}{2\alpha}\epsilon(x)\phi^{b}+\beta\xi=0,\quad-\epsilon(x)\partial_{n}v_{s+1}+\alpha\xi=0,\quad\text{on}\;\partial\Omega.$

These ensure that we always have

$\alpha
v_{s+1}+\beta\partial_{n}v_{s+1}=\frac{1}{2}\phi^{b}\quad\text{on}\;\partial\Omega.$

Take $\psi=\frac{1}{2}\phi-v_{s+1}$ we have

$\nabla\cdot(\epsilon(x)\nabla\psi)=0,x\in\Omega;\alpha\psi+\beta\partial_{n}\psi=0\;\text{on}\;\partial\Omega.$

By the uniqueness of the Poisson problem we conclude $\psi\equiv 0$ or
$\psi={\rm cost}$ if $\alpha=0$, i.e.,

$v_{s+1}\equiv\frac{1}{2}\phi+cost.$

Combing the above we have the following update

$\rho_{i}=\rho^{n}_{i}+\tau\nabla\cdot\left(D_{i}\rho_{i}\nabla(\log(\rho_{i})+z_{i}\phi)\right)+O(\tau^{2}).$

This says scheme (2.24) is a first order time discretization of the PNP system
(2.3). ∎

######  Remark 2.2.

A natural question arises: is the discrete transport still preserves
positivity of $\rho_{i}$. We shall address this issue for the fully discrete
scheme, for which positivity propagation is rigorously established in Theorem
3.2.

##  3\. Numerical method

In this section, we detail the spatial discretization. The underlying
principle for spatial discretization is to preserve the structure of
Wasserstein metric tensor in the discrete sense.

###  3.1. Spatial discretization

We only consider the discretization in one dimensional setting. Let
$\Omega=[a,\ b]$ be the computational domain partitioned into $N$ cells
$I_{j}=[x_{j-\frac{1}{2}},\ \ x_{j+\frac{1}{2}}]$, with mesh size $h=(b-a)/N$
and cell center at $x_{j}=x_{j-\frac{1}{2}}+\frac{1}{2}h$,
$j\in\\{1,2,\cdots,N\\}.$ Let numerical solution be
$\\{\phi_{j}\\}_{j=1}^{N}$, $\\{\rho_{ij}\\}_{j=1}^{N}$, and
$\\{m_{i,j+1/2}\\}_{j=1}^{N-1}$ on two grids $x_{j}$ and $x_{j+1/2}$,
respectively. We define the difference operator by

$(D_{h}v)_{j+1/2}:=\frac{v_{j+1}-v_{j}}{h},\quad(d_{h}v)_{j}=\frac{v_{j+1/2}-v_{j-1/2}}{h}$

and average operator by

$\hat{v}_{j}=\frac{v_{j+1/2}+v_{j-1/2}}{2}.$

We also use $\epsilon_{j+1/2}=\epsilon(x_{j+1/2})$, $f_{j}=f(x_{j})$, and
$D_{ij}=D_{i}(x_{j})$.

The transport constraint is discretized with central difference in space as
follows:

$\rho_{ij}-\rho^{n}_{ij}+d_{h}(m_{i})_{j}=0,$ (3.1)

and the zero boundary conditions $m_{i,1/2}=m_{i,N+1/2}=0$ are applied.

For the Possion equation, we consider the Robin boundary condition at both
ends, other types of boundary conditions can be handled in same fashion. We
introduce two ghost values $\phi_{0}$ and $\phi_{N+1}$ for conveniently
approximating the boundary condition (2.2) with center differences:

$\frac{\phi_{0}+\phi_{1}}{2}-\beta_{a}\epsilon(a)\frac{\phi_{1}-\phi_{0}}{h}=\phi^{b}(a),\quad\frac{\phi_{N+1}+\phi_{N}}{2}+\beta_{b}\epsilon(b)\frac{\phi_{N+1}-\phi_{N}}{h}=\phi^{b}(b).$
(3.2)

This together with the center difference approximation of the Poisson equation
gives a coupled linear system:

$\displaystyle(h+2\beta_{a}\epsilon(a))\phi_{0}+(h-2\beta_{a}\epsilon(a))\phi_{1}-2h\phi^{b}(a)=0,$
(3.3) $\displaystyle-d_{h}(\epsilon
D_{h}\phi)_{j}-f_{j}-\sum_{i=1}^{s}z_{i}\rho_{ij}=0,\quad j=1,\cdots,N,$
$\displaystyle(h-2\beta_{b}\epsilon(b))\phi_{N}+(h+2\beta_{b}\epsilon(b))\phi_{N+1}-2h\phi^{b}(b)=0.$

We denote such linear constraint by $L_{h}(\phi,\rho)=0$. The objective
function then writes as

$\displaystyle\mathcal{F}_{h}(\rho,m,\phi)=$
$\displaystyle\frac{h}{2\tau}\sum_{j=1}^{N}\sum_{i=1}^{s}\frac{\hat{m}^{2}_{i,j}}{\rho_{i,j}}D^{-1}_{i,j}+h\sum_{j=1}^{N}\left(\sum_{i=1}^{s}\rho_{i,j}\log\rho_{i,j}+\frac{\epsilon_{j}}{8h^{2}}(\phi_{j+1}-\phi_{j-1})^{2}\right)$
(3.4)
$\displaystyle+\frac{1}{8\beta_{a}}(\phi_{0}+\phi_{1})^{2}+\frac{1}{8\beta_{b}}(\phi_{N}+\phi_{N+1})^{2},$

which is a second order spatial approximation of the objective functional in
(2.24).

To formulate an admissible set for the discrete minimization problem, let the
discrete probability distribution set be: for $\delta>0$

$P_{h,\delta}=\left\\{(\rho_{1},\cdots,\rho_{N}):\quad\rho_{j}\geq\delta,\quad
h\sum_{j=1}^{N}\rho_{j}=1\right\\}.$

Then the constraint set for $(\rho,\phi)$ becomes

$\mathcal{A}_{h,\delta}=\\{(\rho,\phi):\quad\rho\in[P_{h,\delta}]^{s},\quad
L_{h}(\phi,\rho)=0\\}.$

Thus the admissible set for all $(\rho,m,\phi)$ collectively can be written as

$V_{h,\delta}^{n}=\\{(\rho,m,\phi):\quad\rho_{ij}-\rho_{ij}^{n}+d_{h}(m_{i})_{j}=0,\quad(\rho,\phi)\in\mathcal{A}_{h,\delta}\\}$

with $m_{i,1/2}=m_{i,N+1/2}=0$. Thus we have

$V_{h,\delta}^{n}\subset\mathbb{R}^{s(2N-1)+N+2}.$

The one time update with the fully discrete scheme is to find

$\displaystyle\rho^{n+1}=\arg$ $\displaystyle\min_{u\in
V_{h,\delta}^{n}}\bigg{\\{}\mathcal{F}_{h}(u)\bigg{\\}},\quad
u:=(\rho,m,\phi).$ (3.5)

######  Theorem 3.1.

(Unique solvability) Fix $\tau>0,h>0$ and $\\{\rho_{i}^{n}\in
P_{h,\delta}\\}_{i=1}^{s}$ for some $\delta>0$. Then the function
$\mathcal{F}_{h}(\rho,m,\phi)$ admits a unique minimizer in
$V_{h,\delta}^{n}\subset\mathbb{R}^{s(2N-1)+N+2}$.

###### Proof.

The proof proceeds in two steps:
Step 1 (Admissible set is non-empty and convex) The conservative form of the
transport constraint ensures that we always have

$h\sum_{j=1}^{N}\rho_{ij}=1\quad i\in[s].$

For fixed $\delta>0$, take $\rho_{ij}\geq\delta$, we can uniquely determine
$m$ by

$m_{i,j+1/2}=\frac{1}{h}\sum_{l=1}^{j}(\rho_{il}-\rho_{il}^{n}),$ (3.6)

for $j=1,\cdots,N-1$. From the linear system $L(\phi\;,\rho)=0$ we obtain a
unique $\phi=(\phi_{0},\cdots,\phi_{N+1})$ in terms of $f_{j}$ and
$\rho_{ij}\geq\delta$, since its coefficient matrix is tridiagonal, and
diagonally dominated. Hence the admissible set $V_{h,\delta}^{n}$ is non-
empty. The fact that both the transport constraint and $L(\phi,\rho)=0$ are
linear implies that the set $V_{h,\delta}^{n}$ is convex in
$\mathbb{R}^{s(2N-1)+N+2}$.
Step 2 (Objective function is strictly convex under constraints)
With $u=(\rho,m,\phi)$, for any $u^{0},u^{1}\in V_{h,\delta}$ and
$\theta\in(0,1)$, $u(\theta)=\theta u^{0}+(1-\theta)u^{1}$ is a convex linear
combination of $u^{0}$ and $u^{1}$. In addition, as argued in the proof of
Theorem 2.3, we have $u(\theta)\in V_{h,\delta}$. We now show the convexity of
$\mathcal{F}_{h}(u)$ by directly calculating

$\mathcal{F}_{h}(u(\theta))-\theta\mathcal{F}_{h}(u^{0}))-(1-\theta)\mathcal{F}_{h}(u^{1})=-\theta(1-\theta)(I_{1}+I_{2}+I_{3}),$

where applying Lemma 2.1 to each term $I_{i}$, we have

$\displaystyle I_{1}=$
$\displaystyle\frac{h}{2\tau}\sum_{i=1}^{s}\sum_{j=1}^{N}\frac{(\rho^{1}_{i,j}\hat{m}^{0}_{i,j}-\rho^{0}_{i,j}\hat{m}^{1}_{i,j})^{2}}{\rho^{0}_{i,j}\rho^{1}_{i,j}\rho(\theta)_{i,j}}\geq
0,\quad\text{by}\;(iii)\;\text{of Lemma}\;\ref{lem2.1}$ $\displaystyle I_{2}=$
$\displaystyle
h\sum_{i=1}^{s}\sum_{j=1}^{N}g_{i,j}(\rho^{0}_{i,j},\rho^{1}_{i,j},\theta)(\rho^{0}_{i,j}-\rho^{1}_{i,j})^{2}\geq
0,\;\;\text{by}\;(ii)\;\text{of Lemma}\;\ref{lem2.1}$ $\displaystyle I_{3}=$
$\displaystyle\frac{1}{8h}\sum_{j=1}^{N}\epsilon_{j}[(\phi^{0}_{j+1}-\phi^{0}_{j-1})-(\phi^{1}_{j+1}-\phi^{1}_{j-1})]^{2}\quad\text{by}\;(i)\;\text{of
Lemma}\;\ref{lem2.1}$
$\displaystyle+\frac{1}{8\beta_{a}}[(\phi_{0}^{0}+\phi_{1}^{0})-(\phi_{0}^{1}+\phi_{1}^{1})]^{2}+\frac{1}{8\beta_{b}}[(\phi^{0}_{N}+\phi^{0}_{N+1})-(\phi^{1}_{N}+\phi^{1}_{N+1})]^{2}\geq
0.$

Convexity of $\mathcal{F}_{h}$ follows from $I_{1}+I_{2}+I_{3}\geq 0$. To
establish strictly convexity we only need to show $I_{1}+I_{2}+I_{3}=0$ must
lead to $u^{0}=u^{1}$. We argue as follows.

Clearly the equality holds only when $I_{1}=I_{2}=I_{3}=0$. From $I_{2}=0$ it
follows $\rho^{0}=\rho^{1}$. This when combined with $I_{1}=0$ implies
$\hat{m}^{0}_{i,j}=\hat{m}^{1}_{i,j}$, which together with
$m^{l}_{i,1/2}=m^{l}_{i,N+1/2}=0$ yields $m^{0}=m^{1}$. Finally we show
$\phi^{0}=\phi^{1}$ must also hold. Set $\xi_{j}=\phi_{j}^{0}-\phi_{j}^{1}$
for $j=0,\cdots,N+1,$ then $I_{3}=0$ corresponds to the system of linear
equations $\xi_{0}+\xi_{1}=0,$ $\xi_{N}+\xi_{N+1}=0$ and
$\xi_{j+1}-\xi_{j-1}=0$, for $j=1,\cdots,N$. This obviously admits non-zero
solutions. From the constraint for $\phi$ near the boundary we have

$\phi_{0}^{0}+\phi_{1}^{0}=2\phi^{b}(a)+\frac{\beta_{a}}{h}\epsilon(a)(\phi^{0}_{1}-\phi^{0}_{0}),\quad\phi_{0}^{1}+\phi_{1}^{1}=2\phi^{b}(a)+\frac{\beta_{a}}{h}\epsilon(a)(\phi^{1}_{1}-\phi^{1}_{0}),$

this implies
$\xi_{0}+\xi_{1}=\frac{\beta_{a}}{h}\epsilon(a)(\xi_{1}-\xi_{0})$. Using also
$\xi_{0}+\xi_{1}=0$, we can conclude $\xi_{0}=\cdots=\xi_{N+1}=0$, therefore
$\phi^{0}=\phi^{1}.$ Hence $\mathcal{F}_{h}(u)$ is strictly convex on
$V_{h,\delta}$. ∎

The last issue is to find a threshold for $\delta$ so to ensure that solution
positivity for the PNP system is propagated at all time steps.

######  Theorem 3.2.

(Positivity propagation) There exists $\delta_{0}>0$ such that the minimizer
does not touch the boundary of $V_{h,\delta}^{n}$ for all
$0<\delta\leq\delta_{0}$. This implies that $\rho^{n}>0$ for all $n>0$ as long
as $\rho^{0}>0$.

###### Proof.

We use a contradiction argument: suppose there exists a minimizer $u^{*}$ to
the optimization problem (3.5) touching the boundary of $V_{h,\delta}^{n}$ at
some grid points $j_{1}<\cdots<j_{k}$ with $1\leq k\leq N-1$ for $\rho_{i}$,
that is

$\rho_{i,j_{1}}^{*}=\cdots=\rho_{i,j_{k}}^{*}=\delta.$

From $h\sum_{j=1}^{N}\rho_{i,j}^{*}=1$, we see that $\delta<\frac{1}{b-a}$.
Since $\mathcal{F}_{h}$ is convex and differentiable, we only need to find
$u\in\mathcal{A}_{h,\delta}$ such that

$\nabla\mathcal{F}_{h}(u^{*})\cdot(u-u^{*})<0.$ (3.7)

Note that both $m$ and $\phi$ can be uniquely determined by $\rho$ from the
constraints, it suffices to first choose $\rho$ and then express all
components of $u$ in terms of $\rho$. Let $\rho_{i,j_{k+1}}^{*}$ be the
maximum component in vector $\rho^{*}_{i}$, using
$h\sum_{j=1}^{N}\rho_{i,j}^{*}=1$ we thus have

$\frac{1}{b-a}<\rho^{*}_{i,j_{k+1}}<\frac{1}{h}=\frac{N}{b-a}.$ (3.8)

Without loss of generality, we assume $j_{k+1}>j_{k}$, and

$\rho_{i,j}^{*}\geq\delta+r_{p}(h),\quad j_{p}<j<j_{p+1},\quad p=1,\cdots,k,$
(3.9)

where $r_{p}(0)=0$ and $r_{p}(h)>0$ for $h>0$ small. This can be justified by
approximation for sufficiently small $h$. Fix $h>0$, we take for
$0<\gamma<\frac{1}{k}(\frac{1}{b-a}-\delta)$,

$\rho_{l,j}=\left\\{\begin{array}[]{ll}\delta+\gamma,&l=i,\;j=j_{1},\cdots,j_{k},\\\
\rho^{*}_{i,j_{k+1}}-\gamma k,&l=i,j=j_{k+1},\\\
\rho^{*}_{l,j},&\text{else}.\end{array}\right.$

Hence $\tilde{u}=u-u^{*}$ can be determined by

$\tilde{\rho}_{l,j}=\rho_{l,j}-\rho^{*}_{lj}=\left\\{\begin{array}[]{ll}\gamma,&l=i,\;j=j_{1},\cdots,j_{k},\\\
-\gamma k,&l=i,j=j_{k+1},\\\ 0,&\text{else}.\end{array}\right.$

Using $\tilde{m}=m-m^{*}$ and formula (3.6) for both $m$ and $m^{*}$, we have

$\tilde{m}_{l,j+1/2}=\frac{1}{h}\sum_{p=1}^{j}\tilde{\rho}_{l,p}=\left\\{\begin{array}[]{ll}\frac{1}{h}b_{j}\gamma,&l=i,\;j_{1}\leq
j\leq j_{k+1}-1,\\\ 0,&\text{else},\end{array}\right.$ (3.10)

for $0\leq b_{j}\leq k$. Hence

$0\leq\hat{\tilde{m}}_{i,j}\leq\frac{k\gamma}{h},\quad j_{1}\leq j\leq
j_{k+1}.$

For $\tilde{\phi}=\phi-\phi^{*}$, using (3.3) for both $\phi$ and $\phi^{*}$,
we obtain $A\tilde{\phi}=[0,z_{i}h^{2}\tilde{\rho}_{i},0]^{\top}$, where the
coefficient matrix $A$ is non-singular, more precisely, $\tilde{\phi}$ solves

$\displaystyle(h+2\beta_{a}\epsilon(a))\tilde{\phi}_{0}+(h-2\beta_{a}\epsilon(a))\tilde{\phi}_{1}=0,$
$\displaystyle-\epsilon_{j-1/2}\tilde{\phi}_{j-1}+2\hat{\epsilon}_{j}\tilde{\phi}_{j}-\epsilon_{j+1/2}\tilde{\phi}_{j+1}=h^{2}z_{i}\tilde{\rho}_{i,j}\quad
j=1,\cdots,N,$
$\displaystyle(h-2\beta_{b}\epsilon(b))\tilde{\phi}_{N}+(h+2\beta_{b}\epsilon(b))\tilde{\phi}_{N+1}=0.$

The solution of this linear system can be expressed as

$\tilde{\phi}_{l}=\gamma h^{2}(c_{l}-kd_{l})z_{i},\quad l=0,1,\cdots,N+1$

for some $c_{l},d_{l}$ depending on the coefficients in the above system. The
above preparation yields

$\displaystyle\nabla\mathcal{F}_{h}(u^{*})\cdot(u-u^{*})=\nabla\mathcal{F}_{h}(u^{*})\cdot\tilde{u}$
$\displaystyle=\sum_{l=1}^{s}\sum_{j=1}^{N}\partial_{\rho_{l,j}}\mathcal{F}_{h}(u^{*})\tilde{\rho}_{l,j}+\sum_{l=1}^{s}\sum_{j=1}^{N-1}\partial_{m_{l,j+1/2}}\mathcal{F}_{h}(u^{*})\tilde{m}_{l,j+1/2}+\sum_{j=0}^{N+1}\partial_{\phi_{j}}\mathcal{F}_{h}(u^{*})\tilde{\phi}_{j}$
$\displaystyle=\gamma\left[\sum_{p=1}^{k}\partial_{\rho_{i,j_{p}}}\mathcal{F}_{h}(u^{*})-k\partial_{\rho_{i,j_{k+1}}}\mathcal{F}_{h}(u^{*})\right]+\sum_{j=j_{1}}^{j_{k+1}-1}\partial_{m_{i,j+1/2}}\mathcal{F}_{h}(u^{*})\tilde{m}_{i,j+1/2}+\sum_{j=0}^{N+1}\partial_{\phi_{j}}\mathcal{F}_{h}(u^{*})\tilde{\phi}_{j}$
$\displaystyle=:I_{1}+I_{2}+I_{3}.$

In order to estimate $I_{1},I_{2},I_{3}$ we also need to bound $u^{*}$ in
terms of $\rho^{*}$. From (3.6) and (3.8) we have

$|m^{*}_{i,j+1/2}|\leq\frac{1}{h}\sum_{l=1}^{j}\rho^{*}_{il}\leq\frac{Nj}{(b-a)h}\Rightarrow|\hat{m}^{*}_{i,j}|\leq\frac{N^{2}}{(b-a)h}=Nh^{-2}.$

For $\phi^{*}$ satisfying
$A\phi^{*}=h[2\phi^{b}(a),h(f+\sum_{i=1}^{s}z_{i}\rho^{*}_{i}),2\phi^{b}(b)]^{\top}$,
we have

$|\phi^{*}|\leq
h|A^{-1}|[2\phi^{b}(a),h(f+\sum_{i=1}^{s}z_{i}\rho^{*}_{i}),2\phi^{b}(b)]^{\top}|=:C^{*}_{\phi}.$

We proceed as follows: from the definition of the objective function

$\displaystyle\mathcal{F}_{h}(u)=$
$\displaystyle\frac{h}{2\tau}\sum_{j=1}^{N}\sum_{i=1}^{s}\frac{\hat{m}^{2}_{i,j}}{\rho_{i,j}}D^{-1}_{i,j}+h\sum_{j=1}^{N}\left(\sum_{i=1}^{s}\rho_{i,j}\log\rho_{i,j}+\frac{\epsilon_{j}}{8h^{2}}(\phi_{j+1}-\phi_{j-1})^{2}\right)$
$\displaystyle+\frac{1}{8\beta_{a}}(\phi_{0}+\phi_{1})^{2}+\frac{1}{8\beta_{b}}(\phi_{N}+\phi_{N+1})^{2},$

given in (3.4) we have

$\displaystyle I_{1}$
$\displaystyle=\gamma\sum_{p=1}^{k}\left[-\frac{h}{2\tau}\cdot\frac{(\hat{m}^{*}_{i,j_{p}})^{2}}{(\rho_{i,j_{p}}^{*})^{2}}D_{i,j_{p}}^{-1}+h(1+\log\rho_{i,j_{p}}^{*})\right]-\gamma
k\left[-\frac{h}{2\tau}\cdot\frac{(\hat{m}^{*}_{i,j_{k+1}})^{2}}{(\rho_{i,j_{k+1}}^{*})^{2}}D_{i,j_{k+1}}^{-1}+h(1+\log\rho_{i,j_{k+1}}^{*})\right]$
$\displaystyle=-\frac{\gamma
h}{2\tau\delta^{2}}\sum_{p=1}^{k}(\hat{m}^{*}_{i,j_{p}})^{2}D_{i,j_{p}}^{-1}+\gamma
hk{\rm log}\delta+\frac{\gamma
kh}{2\tau}\frac{(\hat{m}^{*}_{i,j_{k+1}})^{2}}{(\rho^{*}_{i,j_{k+1}})^{2}}D_{i,j_{k+1}}^{-1}-\gamma
hk{\rm log}\rho^{*}_{i,j_{k+1}}$ $\displaystyle\leq-\frac{\gamma
h}{2\tau\delta^{2}}\sum_{p=1}^{k}(\hat{m}^{*}_{i,j_{p}})^{2}D_{i,j_{p}}^{-1}+\gamma
hk{\rm log}\delta+\frac{\gamma kN^{4}}{2\tau h}D_{i,j_{k+1}}^{-1}+\gamma
hk{\rm log}(b-a)$ $\displaystyle=-\frac{\gamma
h}{2\tau\delta^{2}}\sum_{p=1}^{k}(\hat{m}^{*}_{i,j_{p}})^{2}D_{i,j_{p}}^{-1}+\frac{\gamma
kN^{4}}{2\tau h}D_{i,j_{k+1}}^{-1}+\gamma hk{\rm log}\delta(b-a).$

Next, we estimate $I_{2}$:

$\displaystyle I_{2}$
$\displaystyle=\frac{h}{2\tau}\sum_{j=j_{1}}^{j_{k+1}-1}\left(\frac{\hat{m}^{*}_{i,j}}{\rho^{*}_{i,j}}D_{i,j}^{-1}+\frac{\hat{m}^{*}_{i,j+1}}{\rho^{*}_{i,j+1}}D_{i,j+1}^{-1}\right)\tilde{m}_{i,j+1/2}$
$\displaystyle=\frac{h}{\tau}\sum_{j=j_{1}}^{j_{k+1}-1}\frac{\hat{m}^{*}_{i,j}}{\rho^{*}_{i,j}}D_{i,j}^{-1}\hat{\tilde{m}}_{i,j}+\frac{h}{2\tau}\left(-\frac{\hat{m}^{*}_{i,j_{1}}}{\rho^{*}_{i,j_{1}}}D_{i,j_{1}}^{-1}\tilde{m}_{i,j_{1}-1/2}+\frac{\hat{m}^{*}_{i,j_{k+1}}}{\rho^{*}_{i,j_{k+1}}}D_{i,j_{k+1}}^{-1}\tilde{m}_{i,j_{k+1}-1/2}\right)$
$\displaystyle=\frac{h}{\tau}\sum_{p=1}^{k}\frac{\hat{m}^{*}_{i,j_{p}}}{\delta}D_{i,j_{p}}^{-1}\hat{\tilde{m}}_{i,j_{p}}+\frac{h}{\tau}\sum_{p=1}^{k}\sum_{j=j_{p}+1}^{j_{p+1}-1}\frac{\hat{m}^{*}_{i,j}}{\rho^{*}_{i,j}}D_{i,j}^{-1}\hat{\tilde{m}}_{i,j}+\frac{h}{2\tau}\cdot\frac{\hat{m}^{*}_{i,j_{k+1}}}{\rho^{*}_{i,j_{k+1}}}D_{i,j_{k+1}}^{-1}\tilde{m}_{i,j_{k+1}-1/2}$
$\displaystyle\leq\frac{h}{\tau}\cdot\frac{k\gamma}{h}\left(\sum_{p=1}^{k}\frac{|\hat{m}^{*}_{i,j_{p}}|}{\delta}D_{i,j_{p}}^{-1}+\sum_{p=1}^{k}\sum_{j=j_{p}+1}^{j_{p+1}-1}D_{i,j}^{-1}\frac{Nh^{-2}}{\delta+r_{p}(h)}+\frac{N^{2}}{h}D_{i,j_{k+1}}^{-1}\right)$
$\displaystyle\leq\frac{k\gamma}{\tau}\left(\frac{\eta}{2\delta^{2}}\sum_{p=1}^{k}(\hat{m}^{*}_{i,j_{p}})^{2}D_{i,j_{p}}^{-1}+\frac{1}{2\eta}\sum_{p=1}^{k}D_{i,j_{p}}^{-1}+\frac{N}{h^{2}}\sum_{p=1}^{k}\sum_{j=j_{p}+1}^{j_{p+1}-1}\frac{D_{i,j}^{-1}}{r_{p}(h)}+\frac{N^{2}}{h}D_{i,j_{k+1}}^{-1}\right),\quad\forall\eta>0.$

Take $\eta$ so that $k\eta=h$, we have

$I_{2}\leq\frac{\gamma
h}{2\tau\delta^{2}}\sum_{p=1}^{k}(\hat{m}^{*}_{i,j_{p}})^{2}D_{i,j_{p}}^{-1}+C_{1},\quad
C_{1}:=\frac{k^{2}\gamma}{2\tau
h}\sum_{p=1}^{k}D_{i,j_{p}}^{-1}+\frac{Nk\gamma}{\tau
h^{2}}\sum_{p=1}^{k}\sum_{j=j_{p}+1}^{j_{p+1}-1}\frac{D_{i,j}^{-1}}{r_{p}(h)}+\frac{N^{2}k\gamma}{\tau
h}D_{i,j_{k+1}}^{-1}.$

Note that

$\partial_{\phi_{j}}\mathcal{F}_{h}(u^{*})=\frac{1}{4h}\left[\epsilon_{j-1}(\phi^{*}_{j}-\phi^{*}_{j-2})+\epsilon_{j+1}(\phi^{*}_{j}-\phi^{*}_{j+2})\right]\quad
j=2,\cdots,N-1,$

this together with derivatives involving boundary terms allows us to estimate
$I_{3}$:

$\displaystyle I_{3}$
$\displaystyle\leq|\partial_{\phi}\mathcal{F}_{h}(u^{*})|\cdot|\tilde{\phi}|$
$\displaystyle\leq\frac{1}{h}(|\epsilon|+\beta_{a}^{-1}+\beta_{b}^{-1})C_{\phi}^{*}\cdot\gamma|z|(|c|+k|d|)h^{2}=C_{0}C^{*}_{\phi}h.$

For $\delta<\frac{1}{2(b-a)}$, we can take $\gamma k=\frac{1}{2(b-a)}$ such
that $k\gamma<1/(b-a)-\delta$ still holds. Hence

$\displaystyle I_{1}+I_{2}+I_{3}$ $\displaystyle\leq\gamma hk{\rm
log}\delta(b-a)+\frac{\gamma kN^{4}}{2\tau
h}D_{i,j_{k+1}}^{-1}+C_{1}+C_{0}C^{*}_{\phi}h$ $\displaystyle=\frac{1}{2N}{\rm
log}\delta(b-a)+\frac{N^{3}}{4\tau
h^{2}}D_{i,j_{k+1}}^{-1}+C_{1}+C_{0}C^{*}_{\phi}h<0$

provided $\delta<\delta_{0}$ with

$\delta_{0}:=\frac{1}{b-a}\min\left\\{\exp\left(-\frac{N^{4}}{2\tau
h^{2}}D_{i,j_{k+1}}^{-1}-2NC_{1}-2C_{0}C_{\phi}^{*}Nh\right),\frac{1}{2}\right\\}.$

This gives (3.7) as we intended to show. Such contradiction allows us to
conclude that a minimizer at nth step can only occur in the interior of
$V_{h,\delta_{0}}^{n}$ for some $\delta_{0}>0$. In order to show such solution
positivity can propagate, we start from $\rho^{0}>0$. Based on the above
conclusion we recursively have

$\rho^{n+1}\in V_{h,\delta}^{n}\subset V_{h,\delta_{0}}^{n}.$

This completes the proof. ∎

##  4\. Optimization algorithms

In this section, we discuss numerical techniques for solving the constrained
optimization problem (3.5). Let $u=(\rho,m,\phi)$, (3.5) can be written as

$\min_{u}\mathcal{F}_{h}(u),s.t.\quad Au=b,\quad Su\geq\delta,$ (4.1)

Where $\mathcal{F}_{h}(u)$ is defined in (3.4), $Au=b$ is the linear system
corresponding to the constraints (3.1) and (3.3), and $S$ is the selection
matrix that only selects $\rho$ component in $u$.

A simple method to solve (4.1) is to do the following update:

$\tilde{u}^{n+1}=u^{n}-\eta G\nabla_{u}\mathcal{F}_{h}(u^{n}),$

with the matrix

$G=I-A^{\top}(AA^{\top})^{-1}A,$

so that $A\tilde{u}^{n+1}=b$ if $Au^{n}=b$. One then applies the projection

$u^{n+1}=\Pi(\tilde{u}^{n+1})$

so that $\rho_{ij}^{n+1}\geq\delta$.

Input: $A,b$, $u^{n}$, $Iter_{max}$, and $\epsilon$.

Output: $u^{n+1}$

initialization;

$G=I-A^{T}(AA^{T})^{-1}A$, $u^{(0)}=u^{n}$.

for  _$k=1:Iter_{max}$_ do

  * •

Compute the update direction by

$\Delta u=-G\nabla_{u}\mathcal{F}_{h}(u^{(k-1)})$

  * •

Use backtracking to determine step size $\eta$;

* •

Update to get

$\tilde{u}=u^{(k-1)}+\eta\Delta u$

* •

Projection $u^{(k)}=\Pi(\tilde{u})$;

if  _$||Au^{(k)}-b||+||\eta\Delta u||\leq\epsilon$_ then

Stop the iteration;

end if

end for

$u^{n+1}=u^{(k)}$.

Algorithm 1 PG Algorithm

The positivity propagation property stated in Theorem 3.2 ensures that
$Su\geq\delta$ will be fulfilled by the scheme as long as $\rho^{0}\geq\delta$
for $\delta$ suitably small. Hence in our numerical tests the projection step
is not enforced, where we select

$\delta=\max\\{\min\\{h^{2},\tau\\},\min\\{\rho^{in}_{i}(x_{j})\\}\\}>0.$

In summary, the numerical solutions $\rho^{n}_{i,j}$ and $\phi^{n}_{j}$ are
updated with the algorithm:

Input: $\rho^{in}_{i}(x)$, final time $T$, and discretization parameters $h$,
$\tau$, $\delta>0$

Output: $\rho^{n}_{i,j}$, $\phi^{n}_{j}$ for $n=1,\cdots,T/\tau$.

initialization: $u^{0}=(\rho^{0},m^{0},\phi^{0})$ with

$\rho^{0}_{ij}=\max\\{\rho^{in}_{i}(x_{j}),\delta\\}$.

$m^{0}_{i,j}=0$, and $\phi^{0}_{j}$ is obtained by solving (3.2) with
$\rho^{0}_{i,j}$.

for  _$n=1:T/\tau$_ do

$\rho^{n+1}=\arg\min_{u\in
V_{h,\delta}^{n}}\bigg{\\{}\mathcal{F}_{h}(u)\bigg{\\}}$ with Algorithm 1.

end for

Algorithm 2 Algorithm for the fully discrete scheme

######  Remark 4.1.

One may also apply the Primal-Dual Interior-Point algorithm (PDIP) [32,
Chapter 19]) to solve the minimization problem in Algorithm 2, which helps to
enforce the positive lower bound for densities, but it is much more expensive
than the PG method, see a comparison in Table 3.

##  5\. Numerical tests

In this section, we present a selected set of numerical tests to demonstrate
the convergence and properties of the proposed scheme. In all tests, the
tolerance for PG method is set as $10^{-6}.$

Errors are measured in the following discrete $l_{2}$ norm:

$err=\left(\sum_{1\leq j\leq N}h|u_{j}^{n}-U_{j}^{n}|^{2}\right)^{1/2}.$

Here $u_{j}^{n}$ and $U_{j}^{n}$ denotes the numerical solutions and reference
solutions at $(x_{j},t_{n})$. In what follows we take
$u^{n}_{j}=\rho_{i,j}^{n}$,or $\phi_{j}^{n}$ at time $t=n\tau.$

###  5.1. 1D multiple species

We apply our scheme to solve the 1D two-species PNP system (1.1) and verify
the proven properties.

######  Example 5.1.

(Accuracy test) We consider the following PNP system

$\displaystyle\partial_{t}\rho_{1}=$
$\displaystyle\partial_{x}\left(\partial_{x}\rho_{1}+\rho_{1}\partial_{x}\phi\right),$
(5.1) $\displaystyle\partial_{t}\rho_{2}=$
$\displaystyle\partial_{x}\left(\partial_{x}\rho_{2}-\rho_{2}\partial_{x}\phi\right),$
$\displaystyle-\partial_{x}^{2}\phi=$ $\displaystyle\rho_{1}-\rho_{2},$

in $[-1,\ 1]$ and $t>0$. This is (1.1) with $D_{1}=D_{2}=\epsilon=1$,
$q_{1}=1$, $q_{2}=-1$, and $f(x)=0.$ The initial and boundary conditions are
chosen as

$\displaystyle\rho_{1}^{in}(x)=2-x^{2},\quad\rho_{2}^{in}(x)=2+sin(\pi x),$
(5.2) $\displaystyle\phi(0,t)=-1,\ \phi(1,t)=1.$

In the accuracy test, we consider the numerical solutions obtained by
$h=1/320$ and $\tau=1/10000$ as the reference solution. We use the time step
$\tau=h$ and $\tau=h^{2}$ to compute numerical solutions. The errors and
orders at $t=0.5$ are listed in Table 1 and Table 2.

Table 1. Accuracy for Example 5.1 with $\tau=h$ h |  $\rho_{1}$ error | order |  $\rho_{2}$ error | order |  $\phi$ error | order
---|---|---|---|---|---|---
1/10 | 2.67958E-02 | - | 9.80117E-03 | - | 1.17890E-03 | -
1/20 | 1.27689E-02 | 1.06937 | 4.12484E-03 | 1.24862 | 5.46161E-04 | 1.11004
1/40 | 6.20098E-03 | 1.04207 | 1.91422E-03 | 1.10758 | 3.18396E-04 | 0.77850
1/80 | 3.04165E-03 | 1.02764 | 9.21957E-04 | 1.05399 | 1.74525E-04 | 0.86739
Table 2. Accuracy for Example 5.1 with $\tau=h^{2}$ h |  $\rho_{1}$ error | order |  $\rho_{2}$ error | order |  $\phi$ error | order
---|---|---|---|---|---|---
1/10 | 9.00817E-03 | - | 3.13854E-03 | - | 1.42932E-03 | -
1/20 | 2.22285E-03 | 2.01882 | 7.47122E-04 | 2.07068 | 3.62387E-04 | 1.97973
1/40 | 5.36781E-04 | 2.05001 | 1.78121E-04 | 2.06849 | 9.15204E-05 | 1.98537
1/80 | 1.15944E-04 | 2.21091 | 3.79348E-05 | 2.23126 | 2.35563E-05 | 1.95798

We see from Table 1 and Table 2 that the scheme is first order in time and
second order in space.

######  Example 5.2.

In this test, still with the initial boundary value problem (5.1)-(5.2), we
show the proven solution properties. We take $h=0.05,\tau=0.01$ to compute the
numerical solutions up to $T=2$. Solutions at $T=0,\ 0.05,\ 0.25,\ 1.5,\ 2$
are given in Figure 1. In Figure 2 are total mass of $\rho_{1}$ and
$\rho_{2}$, and free energy profile. We see from Figure 1 and Figure 2 that
the scheme is positivity preserving, mass conservative, and energy
dissipating.

Figure 1. Solution evolutions for $\rho_{1},\rho_{2}$, and $\phi$.

Figure 2. Energy dissipation and mass conservation

######  Example 5.3.

(Positivity propagation) In this test, we consider the PNP system (5.1) with
following initial and boundary conditions

$\displaystyle\rho_{1}^{in}(x)=\frac{10}{3}\chi_{{}_{[-0.5,0.5]}},\quad\rho_{2}^{in}(x)=2+sin(\pi
x),$ (5.3) $\displaystyle\phi(0,t)=-1,\ \phi(1,t)=1.$

We take $h=0.05,\tau=0.01$ to compute the numerical solutions up to $T=2$.
Solutions at $T=0,\ 0.015,\ 0.1,\ 1,\ 2$ are displayed in Figure 3. In Figure
4 are total mass of $\rho_{1}$, $\rho_{2}$, and free energy profile. From
these results we see that the scheme is positivity preserving, mass
conservative, and energy dissipating. We also observe that steady state
solutions are identical to those in Example 5.2; this suggests that steady
state solutions of the PNP systems with Dirichlet boundary condition only
depends on the total mass and the Dirichlet boundary condition, but not
sensitive to the profile of the initial data.

In Table 3 we compare CPU times (in seconds) for the PDIP method and the PG
method. Here we set $T=0.5$ and choose different number of sub-intervals.

Table 3. CPU times comparison for PDIP method and PG method  h | 1/10 | 1/50 | 1/100 | 1/150 | 1/200 | 1/20 | 1/300
---|---|---|---|---|---|---|---
PDIP | 1.22 | 2.37 | 8.38 | 17.73 | 31.32 | 48.97 | 74.71
PG | 0.19 | 0.39 | 1.15 | 2.02 | 3.12 | 4.56 | 6.38
Figure 3. Solution evolutions for $\rho_{1},\rho_{2}$.

Figure 4. Energy dissipation and mass conservation

######  Example 5.4.

2D single species (Neumann boundary condition). We now apply our scheme to
solve the 2D single-species PNP system

$\displaystyle\partial_{t}\rho=$
$\displaystyle\nabla\cdot\left(\nabla\rho+\rho\nabla\phi\right),$
$\displaystyle-\Delta\phi=$ $\displaystyle\rho+f(x,y),$

on domain $\Omega=[0,1]\times[0,1]$. We consider the initial boundary
conditions

$\rho^{in}(x,y)=-4(x^{2}-x)-8(y^{2}-y),\quad\frac{\partial\phi}{\partial
n}|_{\partial\Omega}=-1.$

The permanent charge $f(x,y)$ is

$f(x,y)=\begin{dcases}32,&\frac{5}{8}\leq
x\leq\frac{7}{8},\quad\frac{5}{8}\leq y\leq\frac{7}{8},\\\ 0,&else.\\\
\end{dcases}$ (5.4)

This problem satisfies the compatibility condition (2.4). We take
$h_{x}=h_{y}=0.025,\tau=0.01$ to compute the numerical solutions up to $T=6$.
Color plot of the solutions at $T=0.01,0.5,1,2,4,6$ are given in Figure 5. In
Figure 6 are total mass of $\rho$ and free energy profile

Figure 5. Solution evolutions for $\rho$.

Figure 6. Energy dissipation and mass conservation

##  6\. Concluding remarks

In this paper, a dynamic mass transport method for the PNP system is
established by drawing ideas from both the JKO-type scheme [19, 20] and the
classical Bennamou-Breiner formulation [3]. The energy estimate resembles the
physical energy law that governs the PNP system in the continuous case, where
the JKO type formulation is an essential component for preserving intrinsic
solution properties. Both mass conservation and the energy stability are shown
to hold, irrespective of the size of time steps. To reduce computational cost,
we use a local approximation for the artificial time in the constraint
transport equation by a one step difference and the integral in time by a one
term quadrature.

Furthermore, by imposing a centered finite difference discretization in
spatial variables, we establish the solvability of the constrained
optimization problem. This also leads to a remarkable result: for any fixed
time step and spatial meth size, density positivity will be propagating over
all time steps, which is desired for any discrete version of the PNP system.

In the previous section, some numerical experiments were carried out to
demonstrate the proven properties of a computed solution. The first experiment
numerically verified that the variational scheme yields convergence to the
solution of the nonlinear PDE with desired accuracy. Secondly, with further
examples the computed solutions are also shown to satisfy the energy law for
the PNP system, mass conservation, and positivity propagation. It is a matter
of future work to prove an error estimate for these numerical solutions. This
is not a standard error analysis due to the nonlinearities in these problems,
as well as the reformulation as a constrained optimization problem.

## Acknowledgments

This research was partially supported by the National Science Foundation under
Grant DMS1812666.

## References

  * [1] L. Ambrosio, N. Gigli, and G. Savaré.  Gradient flows in metric spaces and in the space of probability measures.  Lect. Math., ETH Zürich. Birkhäuser Verlag, Basel., 2005.
  * [2] M.Z. Bazant, K. Thornton, and A. Ajdari.  Diffuse-charge dynamics in electrochemical systems.  Phys. Rev. E., 70: 021506, 2004.
  * [3] J.-D. Benamou and Y. Brenier.  A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem.  Numer. Math., 84: 375–393, 2000.
  * [4] M. Burger, B. Schlake, and M.-T. Wolfram.  Nonlinear Poisson-Nernst-Planck equations for ion flux through confined geometries.  Nonlinearity., 25: 961–990, 2012.
  * [5] J. A. Carrillo, K. Craig, L. Wang, and C. Wei.  Primal dual methods for Wasserstein gradient flows.  Found. Comput. Math., 22:389–443, 2022.
  * [6] D. Chen and R. Eisenberg.  Poisson-Nernst-Planck (PNP) theory of open ionic channels.  Biophys. J., 64: A22, 1993.
  * [7] J. Ding, Z. Wang, and 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: 108864, 2019.
  * [8] J. Ding, C. Wang, and S. Zhou.  Convergence analysis of structure-preserving numerical methods based on Slotboom transformation for the Poisson–Nernst–Planck equations.  arXiv:2202.10931, 2022.
  * [9] B. Eisenberg.  Ionic channels in biological membranes - electrostatic analysis of a natural nanotube.  Contemp. Phys., 39 (6): 447–466, 1998.
  * [10] A. Flavell, M. Machen, R. Eisenberg, J. Kabre, C. Liu, and X. Li.  A conservative finite difference scheme for Poisson-Nernst-Planck equations.  J. Comput. Electron., 15: 1–15, 2013.
  * [11] A. Flavell, J. Kabre, and X. Li.  An energy-preserving discretization for the Poisson-Nernst-Planck equations.  J. Comput. Electron., 16: 431–441, 2017.
  * [12] H. Gao and D. He.  Linearized conservative finite element methods for the Nernst-Planck-Poisson equations.  J. Sci. Comput., 72: 1269–1289, 2017.
  * [13] D. He and K. Pan.  An energy preserving finite difference scheme for the Poisson-Nernst-Planck system.  Appl. Math. Comput., 287: 214–223, 2016.
  * [14] D. He, K. Pan, and X. Yue.  A positivity preserving and free energy dissipative difference scheme for the Poisson-Nernst-Planck system.  J. Sci. Comput., 81: 436–458, 2019.
  * [15] J.W. Hu and X.D. Huang.  A fully discrete positivity-preserving and energy-dissipative finite difference scheme for Poisson-Nernst-Planck equations  Numer. Math., 145: 77–115, 2020.
  * [16] X. Huo, H. Liu, A. E. Tzavaras, and S. Wang.  An Energy Stable and Positivity-Preserving Scheme for the Maxwell-Stefan Diffusion System.  SIAM J. Numer. Anal., 59(5): 2321–2345, 2021.
  * [17] J.W. Jerome.  Consistency of semiconductor modeling: an existence/stability analysis for the sationary Van Roosbroeck system.  SIAM J. Appl. Math., 45 (4): 565–590, 1985.
  * [18] S. Ji and W. Liu.  Poisson-Nernst-Planck systems for ion flow with density functional theory for hard-sphere potential: I-V relations and critical potentials. Part I: analysis.  J. Dyn. Diff. Equat., 24: 955–983, 2012.
  * [19] R. Jordan, D. Kinderlehrer, and F. Otto.  The variational formulation of the Fokker- Plank equation.  SIAM. J. Math. Anal., 29(1):1–17, 1998.
  * [20] D. Kinderlehrer, L. Monsaingeon, and X. Xu.  A Wasserstein gradient flow approach to Poisson-Nernst-Planck equations.  ESAIM: Control Optim. Calc. Var., 23(1):137–164, 2017.
  * [21] H. Liu and W. Maimaitiyiming.  Energy stable and unconditional positive schemes for a reduced Poisson-Nernst- Plank system.  Comm. Comput. Phys., 7(5): 1505–1529, 2020.
  * [22] H. Liu and W. Maimaitiyiming. Efficient, positive, and energy stable schemes for multi-D Poisson–Nernst–Planck Systems.  J. of Sci Comput, 87(3):1–36, 2021.
  * [23] H. Liu and Z. Wang.  A free energy satisfying finite difference method for Poisson-Nernst-Planck equations.  J. Comput. Phys., 268: 363–376, 2014.
  * [24] H. Liu and Z. Wang.  A free energy satisfying discontinues Galerkin method for one-dimensional Poisson-Nernst-Planck systems.  J. Comput. Phys., 328: 413–437, 2017.
  * [25] H. Liu, Z. Wang, P. Yin, and H. Yu.  Positivity-preserving third order DG schemes for Poisson–Nernst–Planck equations.  J. Comput. Phys., 452: 110777, 2022.
  * [26] W.-C. Li, J.-F. Lu, and L. Wang.  Fisher information regularization schemes for Wasserstein gradient flows.  J. Comput. Phys., 416:109449, 2020.
  * [27] W. Liu.  Geometric singular perturbation approach to steady state Poisson-Nernst-Planck systems.  SIAM J. Appl. Math., 65 (3): 754–766, 2005.
  * [28] C. Liu, C. Wang, S. M. Wise, X. Yue, and S. Zhou.  A positivity-preserving, energy stable and convergent numerical scheme for the Poisson-Nernst-Planck system.  Math. Comp., 90(331), 2071–2106, 2021.
  * [29] P.A. Markowich, C.A. Ringhofer, and C. Schmeiser.  Semiconductor Equations.  Springer- Verlag, 1990\.
  * [30] D. Matthes and S. Plazotta.  A variational formulation of the BDF2 method for metric gradient flows.  ESAIM: M2AN, 53: 145–172, 2019 .
  * [31] M.S. Metti, J. Xu, and C. Liu.  Energetically stable discretizations for charge transport and electrokinetic models.  J. Comput. Phys., 306: 1–18, 2016 .
  * [32] J. Nocedal, S. J. Wright.  Numerical Optimization.  Springer, 2006\.
  * [33] J.-H. Park and J. Jerome.  Qualitative properties of steady-state Poisson-Nernst-Planck systems: Mathematical study.  SIAM J. Appl. Math., 57 (3): 609–630, 1997.
  * [34] S. Selberherr.  Analysis and Simulation of Semiconductor Devices.  Springer-Verlag/Wien,New York, 1984\.
  * [35] J. Shen, J. Xu.  Unconditionally positivity preserving and energy dissipative schemes for Poisson-Nernst-Planck equations.  Numer. Math., 148(3): 671–697, 2021.
  * [36] R. Shen, Q. Zhang, ang B. Lu.  An inverse averaging finite element method for solving the size-modified Poisson-Nernst-Planck equations in ion channel simulations.  arXiv:2112.01692, 2021.
  * [37] T. Teorell.  Transport processes and electrical phenomena in ionic membranes.  Progress Biophysics, 3: 305, 1953.