Efficient, positive, and energy stable schemes for multi-D Poisson-Nernst-Planck systems

Hailiang Liu and Wumaier Maimaitiyiming Iowa State University, Mathematics Department, Ames, IA 50011 [email protected] [email protected]

Abstract.

In this paper, we design, analyze, and numerically validate positive and energy-dissipating schemes for solving the time-dependent multi-dimensional system of Poisson-Nernst-Planck (PNP) equations, which has found much use in the modeling of biological membrane channels and semiconductor devices. The semi-implicit time discretization based on a reformulation of the system gives a well-posed elliptic system, which is shown to preserve solution positivity for arbitrary time steps. The first order (in time) fully-discrete scheme is shown to preserve solution positivity and mass conservation unconditionally, and energy dissipation with only a mild $O(1)$ time step restriction. The scheme is also shown to preserve the steady-states. For the fully second order (in both time and space) scheme with large time steps, solution positivity is restored by a local scaling limiter, which is shown to maintain the spatial accuracy. These schemes are easy to implement. Several three-dimensional numerical examples verify our theoretical findings and demonstrate the accuracy, efficiency, and robustness of the proposed schemes, as well as the fast approach to steady states.

Key words and phrases:

Poisson-Nernst-Planck equations, semi-implicit discretization, mass conservation, positivity, energy decay, steady-state

1991 Mathematics Subject Classification:

35Q92, 35J57, 65N08, 65N12, 82D37.

1. Introduction

In this paper, we are concerned with efficient and structure-preserving numerical approximations to a multi-dimensional 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 [43, 6, 7], electrochemical systems [33, 1], and semiconductor devices [30, 38]. In the semiconductor modeling, it is often called the Poisson-drift-diffusion system.

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 we consider the following form


		$\displaystyle\partial_{t}\rho_{i}+\nabla\cdot J_{i}=0,\quad x\in\Omega\subset\mathbb{R}^{d},\quad t>0,$		(1.1a)
		$\displaystyle-J_{i}=D_{i}(x)\bigg{[}\nabla\rho_{i}+\frac{1}{k_{B}T}\rho_{i}(q_{i}\nabla\phi+\nabla\mu_{i})\bigg{]},$		(1.1b)
		$\displaystyle-\nabla\cdot(\epsilon(x)\nabla\phi)=4\pi\left(f(x)+\sum_{i=1}^{m}q_{i}\rho_{i}\right),$		(1.1c)

subject to initial data $\rho_{i}(x,0)=\rho_{i}^{in}(x)\geq 0$ ( $i=1,\cdots,m$ ) and appropriate boundary conditions to be specified in section 2.1. Here $m$ is the number of species, $\rho_{i}=\rho_{i}(x,t)$ is the charge carrier density for the $i$ -th species, and $\phi=\phi(x,t)$ the electrostatic potential. The charge carrier flux is $J_{i}$ , with which $D_{i}(x)$ is the diffusion coefficient, $k_{B}$ is the Boltzmann constant, $T$ is the absolute temperature. The coupling parameter $q_{i}=z_{i}e$ , where $z_{i}$ is the valence (with sign), $e$ is the unit charge. In the Poisson equation, $\epsilon(x)$ is the permittivity, $f(x)$ is the permanent (fixed) charge density of the system. The equations are valid in a bounded domain $\Omega$ with boundary $\partial\Omega$ and for time $t\geq 0$ . For more accurate modeling of collective interactions of charged particles, the chemical potential $\mu_{i}$ is often included and can be modeled by other means (see section 2.3 for more details).

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 methods with comprehensive numerical analysis for the PNP system is highly desirable. This is what we plan to do in this work.

1.1. Related work

In the literature, there are different numerical solvers available for solving both steady and time-dependent PNP problems; see, e.g., [40, 41, 17, 12, 11, 31, 47, 29]. Many existing algorithms were introduced to handle specific issues in complex applications, in which one may encounter different numerical obstacles, such as discontinuous coefficients, singular charges, geometric singularities, and nonlinear couplings to accommodate various phenomena exhibited by biological ion channels. We refer the interested reader to [44] for some variational multiscale models on charge transport and related algorithms.

Solutions to the PNP equations are known to satisfy some important physical properties. It is desirable to maintain these properties at the discrete level, preferably without or with only mild constraints on time step relative to spatial meshes. Under natural boundary conditions, three main properties for the PNP equations are known as (i) Conservation of mass, (ii) Density positivity, and (iii) Free energy dissipation law. The first property requires the scheme to be conservative. The second property is point-wise and also important for the third property. In general, it is rather challenging to obtain both unconditional positivity and discrete energy decay simultaneously. This is evidenced by several recent efforts [26, 8, 14, 32, 9, 27, 10], in which these properties have been partially addressed at the discrete level for PNP equations. With explicit time discretization, the finite difference scheme in [26] preserves solution positivity under a CFL condition $\Delta t=O(\Delta x^{2})$ and the energy decay was shown for the semi-discrete scheme (time is continuous). An arbitrary high order DG scheme in [27] was shown to dissipate the free energy, with solution positivity restored with the aid of a scaling limiter. With implicit time discretization, the second order finite difference scheme in [8] preserves positivity under a CFL condition $\Delta t=O(\Delta x^{2})$ and a constraint on spatial meshes. An energy-preserving version was further given in [9] with a proven second order energy decay rate. The finite element method in [32] employs the fully implicit backward Euler scheme to obtain solution positivity and the discrete energy decay. In some cases, electric energy alone can be shown to decay (see [27]). Such decay has been verified for the finite difference scheme in [14] and the finite element scheme in [10], both with semi-implicit time discretization.

More recent attempts have focused on semi-implicit schemes based on a formulation of the nonlogarithmic Landau type. As a result, all schemes obtained in [22, 23, 15, 5, 16] have been shown to feature unconditional positivity ( see further discussion in section 1.2).

Our goal here is to construct and analyze structure-preserving numerical schemes for PNP equations in a more general setting: multi-dimension, multi-species, also subject to other chemical forces.

1.2. Our contributions

A key step is to reformulate (1.1a)-(1.1b) as

\partial_{t}\rho_{i}=\nabla\cdot(D_{i}(x)e^{-\psi_{i}}\nabla(\rho_{i}e^{\psi_{i}})),

(1.2)

where

\psi_{i}(x,t)=\frac{q_{i}}{k_{B}T}\phi(x,t)+\frac{1}{k_{B}T}\mu_{i}.

Such reformulation, called the Slotboom transformation in the semiconductor literature, converts a drift-diffusion operator into a self-adjoint elliptic operator. It can be more efficiently solved, and in particular more suitable for keeping the positivity-preserving property. In the context of Fokker-Planck equations it is termed as the nonlogarithmic Landau formulation (see, e.g., [2, 25]). Using such reformulation in [25] Liu and Yu constructed an implicit scheme for a singular Fokker-Planck equation and proved that all three solution properties hold for arbitrary time steps, for which implicit time-discretization is essential. Inspired by [25, 26], we adopted a semi-implicit discretization of (1.2) in [22] to construct a first order in time and second order in space scheme for a reduced PNP system, and proved all three solution properties for the resulting scheme with only a mild $O(1)$ time step restriction. We further introduced a second order (in time) extension in [23] again for the reduced PNP system, and a fully second order scheme in [24] for a class of nonlinear nonlocal Fokker-Planck type equations. All schemes in [22, 23, 24] feature unconditional positivity and a conditional discrete energy dissipation law simultaneously.

This paper improves upon the existing results in [22, 23, 24] in the study of (1.1). We first present a semi-implicit time discretization of form

\frac{\rho^{n+1}_{i}-\rho^{n}_{i}}{\tau}=\nabla\cdot(D_{i}(x)e^{-\psi^{n}_{i}}\nabla(e^{\psi^{n}_{i}}\rho^{n+1}_{i}))=:R[\rho^{n+1}_{i},\psi^{n}_{i}],

(1.3)

which is shown to be well-posed and positivity-preserving for time steps of arbitrary size and independent of the Poisson solver. We further construct the following second order scheme

\frac{\rho^{*}_{i}-\rho^{n}_{i}}{\tau/2}=R[\rho^{*}_{i},\frac{3}{2}\psi^{n}_{i}-\frac{1}{2}\psi^{n-1}_{i}],\quad\rho^{n+1}_{i}=2\rho^{*}_{i}-\rho^{n}_{i},

(1.4)

for which solution positivity for large time steps is restored by a positivity-preserving local limiter. For the spatial discretization we use the 2nd order central difference approximation.

Before stating the main results, let us mention some viable options in the use of reformulation (1.2), i.e.,

\partial_{t}\rho_{i}=R[\rho_{i},\psi_{i}],

which is linear in $\rho_{i}$ if $\psi_{i}$ is a priori given. With the second order central difference in spatial discretization, there are several ways to define $\psi_{i}$ on cell interfaces (see section 3.3). For the time discretization, solution positivity is readily available if we take

\frac{\rho^{n+1}_{i}-\rho^{n-k+1}_{i}}{k\tau}=R[\rho^{n+1}_{i},\psi^{*}_{i}],

(1.5)

with a consistent choice for $\psi_{i}^{*}$ and integer $k\geq 1$ . Different options are introduced in [15, 5, 16] for obtaining their respective positive schemes.

It is natural and simple to take $k=1$ and $\psi^{*}=\psi^{n}$ in (1.5), that is (1.3) (again with further central difference in space). But it is subtle to establish a discrete energy dissipation law. A fully discrete scheme using (1.3) was studied in [5], where no energy dissipation law was established. Nonetheless, a discrete energy dissipation law can be verified with other options. Indeed, (1.5) with $k=2$ and $\psi^{*}_{i}=\psi^{n}_{i}$ was considered in [15], where the authors proved unconditional energy decay for a modified energy. In [16], (1.5) with $k=1$ and $\psi_{i}^{*}=(\psi^{n+1}_{i}+\psi^{n}_{i})/2$ was considered, and all three properties are shown to hold simultaneously even for general boundary conditions for the Poisson equation. Obviously these options can bring further computational overheads.

In this work, we formulate simple finite volume schemes for (1.1) by integrating the central difference method for spatial discretization with the semi-implicit time discretization of the reformulation (1.2). We have strived to advance these numerical schemes by presenting a series of theoretical results. We summarize the main contributions as follows:

•

We show that the first order time discretization gives a well-posed elliptic system (1.3) at each time step, and features solution positivity independent of the time steps (Theorem 3.1). Upper bound of numerical solutions for some cases is established as well (Theorem 3.2).
•

For the first order (in time) fully-discrete scheme, beyond the unconditional solution positivity (Theorem 3.3), we further establish a discrete energy dissipation law for time steps of size $O(1/M)$ , where $M$ is the upper-bound of the numerical solutions (Theorem 3.4). This result sharpens the previous estimates in [22] for the reduced PNP system. We also prove that the scheme preserves steady-states, and numerical solutions converge to a steady state as $n\to\infty$ (Theorem 3.5).
•

We design a fully second order (both in time and space) scheme, and solution positivity is shown for small time steps (Theorem 4.1). While solution positivity for large time steps is ensured by using a local limiter. We prove that such limiter does not destroy the 2nd order spatial accuracy (Theorem 4.2).
•

Three-dimensional numerical tests are conducted to evaluate the scheme performance and verify our theoretical findings. The computational cost of the second order scheme is comparable to that of the first order semi-implicit schemes (see section 5).

1.3. Organization.

We organize this paper as follows: In Section 2, we present primary problem settings and solution properties, as well as model variations. In Section 3, we formulate a unified finite volume method for the PNP system subject to mixed boundary conditions and establish solution positivity, energy dissipation, mass conservation, and steady-state preserving properties for the case of natural boundary conditions. Extension to a second order scheme is given in Section 4. In Section 5, we numerically verify good performance of the schemes. Finally in Section 6 some concluding remarks are given.

Throughout this paper, we denote $\rho$ as vector $(\rho_{1},\cdots,\rho_{m})$ , $\partial\Omega$ as the boundary of domain $\Omega$ includes both the Dirichlet boundary $\partial\Omega_{D}$ and the Neumann boundary $\partial\Omega_{N}$ . $|K|$ denotes the volume of domain $K$ . We use $g_{\alpha}$ to denote ${g_{\alpha}}=1/|K_{\alpha}|\int_{K_{\alpha}}g(x)dx$ , for an integral average of function $g(x)$ over a cell $K_{\alpha}$ .

2. Models and related work

2.1. Boundary conditions

Boundary conditions are a critical component of the PNP model and determine important qualitative behavior of the solution. Here we consider the simplest form of boundary conditions of Dirichlet and/or Neumann type [3].

Let $\Omega$ be a bounded domain with Lipschitz boundary $\partial\Omega$ . The external electrostatic potential $\phi$ is influenced by applied potential, which can be modeled by prescribing a Dirichlet boundary condition

\phi(x,t)=\phi^{b}(x,t),\quad x\in\partial\Omega_{D}.

(2.1)

For the remaining part of the boundary $\partial\Omega_{N}=\partial\bar{\Omega}\setminus\partial\Omega_{D},$ a no-flux boundary condition is applied:

\epsilon(x)\nabla\phi\cdot\mathbf{n}=0,\quad x\in\partial\Omega_{N}.

(2.2)

This boundary condition models surface charges, where $\mathbf{n}$ is the outward unit normal vector on the boundary $\partial\Omega_{N}.$ Same types of boundary conditions are imposed for $\rho_{i}$ as

		$\displaystyle\rho_{i}(x,t)=\rho_{i}^{b}(x,t)\geq 0,\quad x\in\partial\Omega_{D},$		(2.3)
		$\displaystyle J_{i}\cdot\mathbf{n}=0,\quad x\in\partial\Omega_{N}.$		(2.4)

In this work we present our schemes by restricting to a rectangular computational domain $\Omega=(0,L_{1})\times\cdots\times(0,L_{d})$ , with $\partial\Omega_{D}=\{x\in\bar{\Omega}|\quad x_{1}=0,\ x_{1}=L_{1}\}$ .

We remark that the boundary conditions for the electrostatic potential are not unique and greatly depend on the problem under investigation. For example, one may use a non-homogeneous Neumann boundary condition ( $\nabla\phi\cdot\mathbf{n}=\sigma$ is used in [27]) or Robin boundary conditions [8, 16]. The existence and uniqueness of the solution for the nonlinear PNP boundary value problems have been studied in [19, 21, 34] for the 1D case and in [3, 18] for multi-dimensions.

2.2. Positivity and energy dissipation law

One important solution property is

\rho_{i}(x,t)\geq 0,\quad x\in\Omega,\ t>0.

(2.5)

Integration of each density equation gives

\frac{d}{dt}\int_{\Omega}\rho_{i}(x,t)dx=\int_{\partial\Omega}J_{i}\cdot\mathbf{n}ds,

which with zero flux $J_{i}\cdot\mathbf{n}=0$ on the whole boundary leads to the mass conservation:

\int_{\Omega}\rho_{i}(x,t)dx=\int_{\Omega}\rho_{i}^{in}(x)dx,\quad t>0,\quad i=1,\cdots,m.

(2.6)

We consider the free energy functional $E$ associated to (1.1) with $\mu_{i}=\mu_{i}(x)$ :

E=\int_{\Omega}\bigg{(}\sum_{i=1}^{m}\rho_{i}(\log\rho_{i}-1)+\frac{1}{2k_{B}T}(f+\sum_{i=1}^{m}q_{i}\rho_{i})\phi+\frac{1}{k_{B}T}\sum_{i=1}^{m}\rho_{i}\mu_{i}\bigg{)}dx.

(2.7)

In virtue of the Poisson equation (1.1c), the free energy may be written as

E=\int_{\Omega}\bigg{(}\sum_{i=1}^{m}\rho_{i}(\log\rho_{i}-1)+\frac{\epsilon}{8\pi k_{B}T}|\nabla\phi|^{2}+\frac{1}{k_{B}T}\sum_{i=1}^{m}\rho_{i}\mu_{i}\bigg{)}dx.

Note that the unscaled free energy $F=k_{B}TE$ is also often used, see [28]. A formal calculation gives

	$\displaystyle\frac{dE}{dt}$	$\displaystyle=-\int_{\Omega}\sum_{i=1}^{m}D_{i}(x)\rho_{i}\|\nabla\psi^{}_{i}\|^{2}dx+\int_{\partial\Omega}\sum_{i=1}^{m}\psi_{i}^{}J_{i}\cdot\mathbf{n}ds$
		$\displaystyle\quad+\frac{1}{8\pi k_{B}T}\int_{\partial\Omega}\epsilon(x)\left[\phi(\partial_{n}\phi)_{t}-\partial_{n}\phi\phi_{t}\right]ds,$

where

\psi^{*}_{i}:=\log\rho_{i}+\frac{q_{i}}{k_{B}T}\phi+\frac{1}{k_{B}T}\mu_{i}.

Clearly, with $\partial\Omega_{D}=\emptyset$ , we have the following energy dissipation law:

\frac{dE}{dt}=-\int_{\Omega}\sum_{i=1}^{m}D_{i}(x)\rho_{i}|\nabla\psi^{*}_{i}|^{2}dx\leq 0.

(2.8)

Otherwise, the Dirichlet boundary condition needs to be carefully handled (see, e.g., [28]).

For time dependent chemical potentials $\mu_{i}(x,t)$ , the total free energy and its dissipation law needs to be modified depending on how the chemical potential is determined.

2.3. Chemical potential

In application, the chemical potential $\mu_{i}$ often includes the ideal chemical potential $\mu^{id}_{i}(x,t)$ and the excess chemical potential $\mu^{ex}_{i}(x,t)$ of the charged particles:

\mu_{i}(x,t)=\mu^{id}_{i}(x,t)+\mu^{ex}_{i}(x,t),

with

\mu^{id}_{i}(x,t)=-\log[\gamma_{i}\rho_{i}(x,t)/\rho^{bulk}_{i}],

where the activity coefficient $\gamma_{i}$ described by the extended Debye-Hückel theory depends on $\rho$ in nonlinear manner. Meanwhile,

\mu^{ex}_{i}(x,t)=\frac{\delta F^{ex}(\rho(x,t))}{\delta\rho_{i}(x,t)}

is the $L^{2}$ variational derivative of the excess chemical functional $F^{ex}$ , which may include hard-sphere components, short-range interactions, Coulomb interactions and electrostatic correlations, where the expression of each component can be found in [42, 31].

We remark that the steric interactions between ions of different species are important in the modeling of ion channels [20, 17]. Such effects can be described by choosing

F^{\rm ex}=\frac{1}{2}\int_{\Omega}\omega_{ij}\rho_{i}\rho_{j},

where $\omega_{ij}$ are the second-order virial coefficients for hard spheres, depending on the size of $i$ -th and $j$ -th ion species [49]. With this addition alone, the flux becomes

-J_{i}=D_{i}(x)\left(\nabla\rho_{i}+\frac{1}{k_{B}T}q_{i}\rho_{i}\nabla\phi+\ \rho_{i}\sum_{j=1}^{m}\omega_{ij}\nabla\rho_{j}\right).

The PNP system with this modified flux has been studied numerically first in [39] without cross steric interactions, and then in [5] with cross interactions.

Our schemes will be constructed so that numerical solutions are updated in an explicit-implicit manner while $\mu_{i}$ needs only to be evaluated off-line. For simplicity, we shall present our schemes assuming $\mu_{i}$ is given while keeping in mind that it can be applied to complex chemical potentials without difficulty.

2.4. Steady states

By the free energy dissipation law (2.8), the solution to (1.1) with zero-flux boundary conditions is expected to converge to the steady-states as time becomes large. In such case the steady states formally satisfy (1.1) with $\partial_{t}\rho_{i}=0$ ; i.e.,

\nabla\cdot(D_{i}(x)\rho_{i}\nabla\psi_{i}^{*}))=\nabla\cdot J_{i}=0,\quad J_{i}\cdot\mathbf{n}=0,\quad x\in\partial\Omega.

This yields $\int_{\Omega}J_{i}\cdot\nabla\psi_{i}^{*}dx=0$ , which ensures that $\psi_{i}^{*}$ must be a constant. This gives the well-known Boltzmann distribution

\rho_{i}=c_{i}e^{-\frac{1}{k_{B}T}(q_{i}\phi+\mu_{i})},

(2.9)

where $c_{i}$ is any constant. Such constant can be uniquely determined by the initial data in the PNP system (1.1) if such steady-state is approached by the solution at large times. Indeed, mass conservation simply gives

c_{i}=\frac{\int_{\Omega}\rho^{in}_{i}dx}{\int_{\Omega}e^{-\frac{1}{k_{B}T}(q_{i}\phi+\mu_{i})}dx}.

(2.10)

This allows us to obtain a closed Poisson-Boltzmann equation (PBE) of form

-\nabla\cdot(\epsilon(x)\nabla\phi)=4\pi\left(f(x)+\sum_{i=1}^{m}q_{i}c_{i}e^{-\frac{1}{k_{B}T}(q_{i}\phi+\mu_{i})}\right),\quad\partial_{\mathbf{n}}\phi|_{\partial\Omega}=0.

(2.11)

We should point out that the numerical method presented in this paper may be used as an iterative algorithm to numerically compute the nonlocal PBE (2.11); hence it serves as a simpler alternative to the iterative DG methods recently developed in [45, 46].

In practical applications, one may describe ions of less interest using the Boltzmann distribution and still solve the NP equations for the target ions so to reduce the computational cost, see [48] for further details on related models. Our numerical method thus provides an alternative path to simulate such models.

3. Numerical method

In this section we will construct positive and energy stable schemes.

3.1. Reformulation

By setting

\psi_{i}(x,t)=\frac{1}{k_{B}T}(q_{i}\phi(x,t)+\mu_{i}),

we reformulate the density equation (1.1a)-(1.1b) as:

\partial_{t}\rho_{i}=\nabla\cdot(D_{i}(x)e^{-\psi_{i}}\nabla(e^{\psi_{i}}\rho_{i})).

(3.1)

In spite of the aforementioned advantages of such reformulation, possible large variation of the transformed diffusion coefficients could result in large condition number of the stiffness matrix [29]. This issue has been recently investigated in [36, 5].

3.2. Time discretization

Let $\tau>0$ be a time step, and $t_{n}=\tau n$ , $n=0,1\cdots,$ be the corresponding temporal grids. We initialize by taking $\rho^{0}(x)=\rho^{in}(x)$ , and obtaining $\phi^{0}$ by solving the Poisson equation (1.1c) using $\rho^{0}(x)$ .

Let $\rho^{n}$ and $\phi^{n}$ be numerical approximations of $\rho(x,t_{n})$ and $\phi(x,t_{n})$ , respectively, we first obtain $\rho^{n+1}$ by solving the following elliptic system:


		$\displaystyle\frac{\rho^{n+1}_{i}-\rho^{n}_{i}}{\tau}=\nabla\cdot(D_{i}(x)e^{-\psi^{n}_{i}}\nabla(e^{\psi^{n}_{i}}\rho^{n+1}_{i}))=:R[\rho^{n+1}_{i},\psi^{n}_{i}],$		(3.2a)
		$\displaystyle\rho^{n+1}_{i}=\rho^{b}_{i}(x,t_{n+1}),\quad x\in\partial\Omega_{D},$		(3.2b)
		$\displaystyle\nabla(e^{\psi^{n}_{i}}\rho^{n+1}_{i})\cdot\mathbf{n}=0,\quad x\in\partial\Omega_{N},$		(3.2c)

where

\psi_{i}^{n}=\frac{1}{k_{B}T}(q_{i}\phi^{n}+\mu_{i}).

Using this obtained $\rho^{n+1}$ , we update to obtain $\phi^{n+1}$ from solving


		$\displaystyle-\nabla\cdot(\epsilon(x)\nabla\phi^{n+1})=4\pi\left(f(x)+\sum_{i=1}^{m}q_{i}\rho^{n+1}_{i}\right),$		(3.3a)
		$\displaystyle\phi^{n+1}(x)=\phi^{b}(x,t_{n+1}),\quad x\in\partial\Omega_{D},$		(3.3b)
		$\displaystyle\nabla\phi^{n+1}\cdot\mathbf{n}=0,\quad x\in\partial\Omega_{N}.$		(3.3c)

This scheme is well-defined for any $\tau>0$ with $\rho^{n}\geq 0$ for all $n\in\mathbb{N}$ . More precisely, we have

Theorem 3.1.

Assume $D_{i}(x)\geq D_{0}>0$ and $\epsilon(x)\geq\epsilon_{0}>0$ , and $\mu_{i}(x)\in C(\bar{\Omega})$ . Then for given $(\rho^{n},\phi^{n})\in C(\bar{\Omega})\cap C^{2}(\Omega)$ , there exists a unique solution $(\rho^{n+1},\phi^{n+1})\in C(\bar{\Omega})\cap C^{2}(\Omega)$ . If $\rho^{n}\geq 0$ and $\rho^{b}(x,t)\geq 0$ , $x\in\partial\Omega_{D}$ , then $\rho^{n+1}\geq 0$ for $x\in\Omega$ .

The proof is deferred to the appendix A.

In some cases density for the PNP problem is known to be uniformly bounded for all time. We shall show this bound property also for the semi-discrete scheme (3.2).

Theorem 3.2.

Let $0\leq\rho_{i}^{in}(x)\leq B_{i}$ , $0\leq\rho^{b}_{i}(x,t)\leq B^{b}_{i},$ $D_{i}(x)/\epsilon(x)=\sigma_{i}$ be constants, $\Omega$ be $C^{1}$ convex domain, all $q_{i}$ have the same sign, and $\mu_{i}$ is smooth with $(\nabla\mu_{i})\cdot\mathbf{n}\geq 0$ on $\partial\Omega_{N}$ . If

\tau<\frac{1}{Q_{i,max}},

then $\rho^{n}$ obtained by scheme (3.2) is uniformly bounded, i.e.,

\rho^{n}_{i}(x)\leq\max\left\{B^{b}_{i},\quad B_{i},\quad\frac{Q_{i,max}}{\gamma_{i}}\right\},

(3.4)

where $Q_{i,max}=\max_{x\in\bar{\Omega}}Q_{i}(x)$ with

Q_{i}(x)=\frac{1}{k_{B}T}\left[\nabla\cdot\left(D_{i}(x)\nabla\mu_{i}\right)-4\pi q_{i}\sigma_{i}f(x)\right],\quad\gamma_{i}=\frac{4\pi q_{i}^{2}\sigma_{i}}{k_{B}T}.

Remark 3.1.

In the case of $q_{i}$ with different sign, density $\rho_{i}$ in (1.1) may not be bounded.

Proof.

We rewrite the semi-discrete scheme

\frac{\rho^{n+1}_{i}-\rho^{n}_{i}}{\tau}=\nabla\cdot\left(D_{i}(x)e^{-\psi^{n}_{i}}\nabla\left(\rho^{n+1}_{i}e^{\psi^{n}_{i}}\right)\right)

into

\frac{\rho^{n+1}_{i}-\rho^{n}_{i}}{\tau}=D_{i}(x)\Delta\rho^{n+1}_{i}+b_{i}\cdot\nabla\rho^{n+1}_{i}+c_{i}\rho^{n+1}_{i},

where

b_{i}=\left(\nabla D_{i}(x)+D_{i}(x)\nabla\psi^{n}_{i}\right),\quad c_{i}=\nabla\cdot\left(D_{i}(x)\nabla\psi^{n}_{i}\right).

In virtue of $\psi^{n}_{i}=\frac{q_{i}}{k_{B}T}\phi^{n}+\frac{1}{k_{B}T}\mu_{i}$ and $D_{i}(x)/\epsilon(x)=\sigma_{i}$ , the coefficient $c_{i}$ can be estimated as

	$\displaystyle c_{i}=$	$\displaystyle\frac{1}{k_{B}T}\left[\nabla\cdot\left(q_{i}D_{i}(x)\nabla\phi^{n}\right)+\nabla\cdot\left(D_{i}(x)\nabla\mu_{i})\right)\right]$
	$\displaystyle=$	$\displaystyle\frac{1}{k_{B}T}\left[q_{i}\sigma_{i}\nabla\cdot\left(\epsilon(x)\nabla\phi^{n}\right)+\nabla\cdot\left(D_{i}(x)\nabla\mu_{i})\right)\right]$
		$\displaystyle(\text{using (\ref{Pst1})})$
	$\displaystyle=$	$\displaystyle\frac{1}{k_{B}T}\left[-4\pi q_{i}\sigma_{i}\left(f(x)+\sum_{j=1}^{m}q_{j}\rho^{n}_{j}\right)+\nabla\cdot\left(D_{i}(x)\nabla\mu_{i})\right)\right]$
		$\displaystyle(\text{using $q_{i}q_{j}>0$ and $\rho_{j}^{n}\geq 0$})$
	$\displaystyle\leq$	$\displaystyle\frac{1}{k_{B}T}\left[\nabla\cdot\left(D_{i}(x)\nabla\mu_{i}\right)-4\pi q_{i}\sigma_{i}f(x)-4\pi q_{i}^{2}\sigma_{i}\rho^{n}_{i}\right]$
	$\displaystyle=$	$\displaystyle Q_{i}(x)-\gamma_{i}\rho^{n}_{i}.$

Hence

\frac{\rho^{n+1}_{i}-\rho^{n}_{i}}{\tau}\leq D_{i}(x)\Delta\rho^{n+1}_{i}+b_{i}\cdot\nabla\rho^{n+1}_{i}+\rho^{n+1}_{i}\left(Q_{i,max}-\gamma_{i}\rho^{n}_{i}\right).

(3.5)

We proceed to distinct three cases, by letting $x^{*}=\operatornamewithlimits{argmax}_{x\in\bar{\Omega}}\rho^{n+1}_{i}(x)$ :

(i) If $x^{*}\in\partial\Omega_{D}$ we have

\rho^{n+1}_{i}(x^{*})=\rho^{b}_{i}(x^{*},t_{n+1})\leq B^{b}_{i}.

(ii) If $x^{*}\in\Omega$ , then (3.5) can be reduced to

\frac{\rho^{n+1}_{i}(x^{*})-\rho^{n}_{i}(x^{*})}{\tau}\leq\rho^{n+1}_{i}(x^{*})\left(Q_{i,max}-\gamma_{i}\rho^{n}_{i}(x^{*})\right).

This using notation $\rho^{n}_{i,max}=\max_{x\in\bar{\Omega}}\rho^{n}_{i}$ yields

\rho^{n+1}_{i}(x)\leq\rho^{n+1}_{i}(x^{*})\leq\frac{\rho^{n}_{i,max}}{1-\tau Q_{i,max}+\tau\gamma_{i}\rho^{n}_{i,max}}=:P(\rho^{n}_{i,max}),

(3.6)

where we used the fact that $P(\cdot):\mathbb{R}^{+}\to\mathbb{R}^{+}$ is non-decreasing.

(iii) If $x^{*}\in\partial\Omega_{N}$ , we must have $\rho^{n+1}_{i}(x^{*})\leq P(\rho^{n}_{i,max})$ . Otherwise assume $\rho^{n+1}_{i}(x^{*})>P(\rho^{n}_{i,max}).$ Set

U(x)=\rho^{n+1}_{i}(x)-\rho^{n+1}_{i}(x^{*}),

and introduce the differential operator

L\xi:=\tau D_{i}(x)\Delta\xi+\tau b_{i}\cdot\nabla\xi-(1-\tau Q_{i,max}+\tau\gamma_{i}\rho^{n}_{i})\xi.

From (3.5) we have

L\rho_{i}^{n+1}\geq-\rho_{i}^{n},

and using (3.6) we obtain

	$\displaystyle LU(x)=$	$\displaystyle L\rho^{n+1}_{i}(x)-L\rho^{n+1}_{i}(x^{*})$
	$\displaystyle\geq$	$\displaystyle-\rho^{n}_{i}+(1-\tau Q_{i,max}+\tau\gamma_{i}\rho^{n}_{i})\rho^{n+1}_{i}(x^{*})$
	$\displaystyle\geq$	$\displaystyle-\rho^{n}_{i}+(1-\tau Q_{i,max}+\tau\gamma_{i}\rho^{n}_{i})P(\rho^{n}_{i,max})$
	$\displaystyle\geq$	$\displaystyle 0.$

Note that $U(x)\leq 0$ on $\partial\Omega$ and $U(x^{*})=0$ . Apply the maximum-principle [35, Theorem 8] we have

(\nabla U(x^{*}))\cdot\mathbf{n}=(\nabla\rho^{n+1}_{i}(x^{*}))\cdot\mathbf{n}>0.

On the other hand, from the no-flux boundary condition (3.2c) and using (3.3c), we have

	$\displaystyle 0=$	$\displaystyle\left(\nabla\left(\rho^{n+1}_{i}e^{\psi^{n}_{i}}\right)\right)\cdot\mathbf{n}$
	$\displaystyle=$	$\displaystyle\left(e^{\psi^{n}_{i}}\nabla\rho^{n+1}_{i}+\frac{1}{k_{B}T}e^{\psi^{n}_{i}}\rho^{n+1}_{i}(q_{i}\nabla\phi^{n}+\nabla\mu_{i})\right)\cdot\mathbf{n}$
	$\displaystyle=$	$\displaystyle e^{\psi^{n}_{i}}\left(\nabla\rho^{n+1}_{i}\cdot\mathbf{n}+\frac{1}{k_{B}T}\nabla\mu_{i}\cdot\mathbf{n})\right)$
	$\displaystyle>$	$\displaystyle e^{\psi^{n}_{i}}\frac{1}{k_{B}T}\nabla\mu_{i}\cdot\mathbf{n}\quad x\in\partial\Omega_{N}.$

This is a contradiction to the assumption $(\nabla\mu_{i})\cdot\mathbf{n}\geq 0$ . Hence for $x\in\Omega\cup\partial\Omega_{N}\cup\partial\Omega_{D}=\bar{\Omega}$ , we have

\rho^{n+1}_{i,max}\leq\max\left\{B^{b}_{i},\quad P(\rho^{n}_{i,max})\right\}.

Again by the monotonicity of $P(\cdot)$ , we obtain

\rho^{n+1}_{i,max}\leq\max\left\{B^{b}_{i},\quad\max\{\rho^{n}_{i,max},\quad\frac{Q_{i,max}}{\gamma_{i}}\}\right\}.

The stated result (3.4) thus follows by induction. ∎

A discrete energy dissipation law can be established by precisely quantifying a sufficient bound on the time step. In order to save space, we present a detailed analysis of the energy dissipation property only for the fully discrete scheme in the next section.

3.3. Spatial discretization

For given positive integers $N_{j}$ $(j=1,\cdots,d)$ , let $h_{j}=L_{j}/N_{j}$ be the mesh size in $j$ -th direction, $\alpha\in\mathbb{Z}^{d}$ be the index vector with $\alpha(j)\in\{1,\cdots,N_{j}\}$ , and $e_{j}\in\mathbb{Z}^{d}$ be a vector with $j$ -th entry equal to one and all other entries equal to zero. We partition the domain $\Omega$ into computational cells

K_{\alpha}=[(\alpha(1)-1)h_{1},\alpha(1)h_{1}]\times\cdots\times[(\alpha(d)-1)h_{d},\alpha(d)h_{d}]

with cell size $|K_{\alpha}|=\prod_{j=1}^{d}h_{j}$ such that $\bigcup_{\alpha\in\mathcal{A}}K_{\alpha}=\Omega,$ where $\mathcal{A}$ denotes the set of all indices $\alpha$ .

3.3.1. Density update

A finite volume approximation of (3.2a) over each cell $K_{\alpha}$ with $\alpha\in\mathcal{A}$ gives

	$\frac{\rho^{n+1}_{i,\alpha}-\rho^{n}_{i,\alpha}}{\tau}=\sum_{j=1}^{d}\frac{C_{i,\alpha+e_{j}/2}-C_{i,\alpha-e_{j}/2}}{h_{j}}=:R_{\alpha}[\rho^{n+1}_{i},\psi^{n}_{i}],$		(3.7a)
where $\rho^{0}_{i,\alpha}:=\rho^{in}_{i,\alpha}.$

Numerical fluxes on interfaces are defined by:

(i) on the interior interfaces,

\displaystyle C_{i,\alpha+e_{j}/2}=\frac{D_{i}(x_{\alpha+e_{j}/2})e^{-\psi^{n}_{i,\alpha+e_{j}/2}}}{h_{j}}\left(\rho^{n+1}_{i,\alpha+e_{j}}e^{\psi^{n}_{i,\alpha+e_{j}}}-\rho^{n+1}_{i,\alpha}e^{\psi^{n}_{i,\alpha}}\right),\quad\text{ for }1<\alpha(j)<N_{j};

(3.7b)

(ii) on the boundary $\partial\Omega_{D}$ ,

		$\displaystyle C_{i,\alpha-e_{1}/2}=\frac{2D_{i}(x_{\alpha-e_{1}/2})e^{-\psi^{b}_{i}(x_{\alpha-e_{1}/2},t_{n})}}{h_{1}}\bigg{(}\frac{\rho^{n+1}_{i,\alpha}}{e^{-\psi^{n}_{i,\alpha}}}-\frac{\rho^{b}_{i}(x_{\alpha-e_{1}/2},t_{n+1})}{e^{-\psi^{b}_{i}(x_{\alpha-e_{1}/2},t_{n})}}\bigg{)},\quad\alpha(1)=1,$		(3.7c)
		$\displaystyle C_{i,\alpha+e_{1}/2}=\frac{2D_{i}(x_{\alpha+e_{1}/2})e^{-\psi^{b}_{i}(x_{\alpha+e_{1}/2},t_{n})}}{h_{1}}\bigg{(}\frac{\rho^{b}_{i}(x_{\alpha+e_{1}/2},t_{n+1})}{e^{-\psi^{b}_{i}(x_{\alpha+e_{1}/2},t_{n})}}-\frac{\rho^{n+1}_{i,\alpha}}{e^{-\psi^{n}_{i,\alpha}}}\bigg{)},\quad\alpha(1)=N_{1};$		(3.7c)

(iii) on the boundary $\partial\Omega_{N}$ ,

		$\displaystyle C_{i,\alpha-e_{j}/2}=0,\ \text{ for }\alpha(j)=1,$		(3.7d)
		$\displaystyle C_{i,\alpha+e_{j}/2}=0,\ \text{ for }\alpha(j)=N_{j}.$		(3.7d)

In (3.7b) $e^{-\psi^{n}_{i,\alpha+e_{j}/2}}$ needs to be evaluated using numerical solutions $\phi^{n}_{\alpha}$ . There are three choices, all are second order approximations:

(i) the harmonic mean

e^{-\psi^{n}_{i,\alpha+e_{i}/2}}=\frac{2e^{-\psi^{n}_{i,\alpha+e_{j}}-\psi^{n}_{i,\alpha}}}{e^{-\psi^{n}_{i,\alpha+e_{j}}}+e^{-\psi^{n}_{i,\alpha}}},

(3.8)

(ii) the geometric mean

e^{-\psi^{n}_{i,\alpha+e_{i}/2}}=\sqrt{e^{-\psi^{n}_{i,\alpha+e_{j}}-\psi^{n}_{i,\alpha}}},

(3.9)

(iii) the algebraic mean

e^{-\psi^{n}_{i,\alpha+e_{i}/2}}=\frac{e^{-\psi^{n}_{i,\alpha+e_{j}}}+e^{-\psi^{n}_{i,\alpha}}}{2}.

(3.10)

It is reported in [36] that the harmonic mean results in a linear system with better condition number than that of the geometric mean. We use the harmonic mean in our numerical tests.

3.3.2. Solving Poisson’s equation

In order to complete the scheme, we need to evaluate $\psi_{i,\alpha}^{n}$ by

\psi_{i,\alpha}^{n}=\frac{1}{k_{B}T}(q_{i}\phi_{\alpha}^{n}+\mu_{i,\alpha}),

and $\phi^{n}_{\alpha}$ is determined from $\rho^{n}_{\alpha}$ by using the following discretization of the equation (3.3a):

	$-\sum_{j=1}^{d}\frac{\Phi^{n}_{\alpha+e_{j}/2}-\Phi^{n}_{\alpha-e_{j}/2}}{h_{j}}=4\pi\left(f_{\alpha}+\sum_{i=1}^{m}q_{i}\rho^{n}_{i,\alpha}\right),$		(3.11a)
where numerical fluxes on cell interfaces are defined by:

(i) on the interior interfaces,

\displaystyle\Phi^{n}_{\alpha+e_{j}/2}

\displaystyle=\epsilon(x_{\alpha+e_{j}/2})\frac{\phi^{n}_{\alpha+e_{j}}-\phi^{n}_{\alpha}}{h_{j}},\quad\text{ for }1<\alpha(j)<N_{j},

(3.11b)

(ii) on the boundary $\partial\Omega_{D}$ ,

	$\displaystyle\Phi^{n}_{\alpha-e_{1}/2}$	$\displaystyle=\epsilon(x_{\alpha-e_{1}/2})\frac{2(\phi^{n}_{\alpha}-\phi^{b}(x_{\alpha-e_{1}/2},t_{n}))}{h_{1}},\quad\text{ for }\alpha(1)=1,$		(3.11c)
	$\displaystyle\Phi^{n}_{\alpha+e_{1}/2}$	$\displaystyle=\epsilon(x_{\alpha+e_{1}/2})\frac{2(\phi^{b}(x_{\alpha+e_{1}/2},t_{n})-\phi^{n}_{\alpha})}{h_{1}},\quad\text{ for }\alpha(1)=N_{1},$		(3.11c)

(iii) on the boundary $\partial\Omega_{N}$ ,

		$\displaystyle\Phi^{n}_{\alpha-e_{j}/2}=0,\ \text{ for }\alpha(j)=1,$		(3.11d)
		$\displaystyle\Phi^{n}_{\alpha+e_{j}/2}=0,\ \text{ for }\alpha(j)=N_{j}.$		(3.11d)

Note that in the case of $\partial\Omega_{D}=\emptyset$ , the solution to (3.11) is unique only up to an additive constant, in such case we take $\phi^{n}_{(1,\cdots,1)}=0$ to obtain a unique solution $\phi^{n}_{\alpha}$ .

3.3.3. Positivity

The following theorem states that the scheme (3.7) preserves positivity of numerical solutions without any time step restriction.

Theorem 3.3.

Let $\rho^{n+1}_{\alpha}$ be obtained from (3.7). If $\rho^{n}_{\alpha}\geq 0$ for all $\alpha\in\mathcal{A}$ , and $\rho^{b}(x,t_{n})\geq 0$ , $x\in\partial\Omega_{D}$ , then

\rho^{n+1}_{\alpha}\geq 0\quad\text{for all}\quad\alpha\in\mathcal{A}.

Proof.

This proof mimics that in [25] for the Fokker-Planck equation. Set $\lambda_{j}=\frac{\tau}{h_{j}^{2}}$ , $\bar{D}_{i,\alpha+e_{j}/2}=D_{i}(x_{\alpha+e_{j}/2})e^{-\psi^{n}_{i,\alpha+e_{j}/2}}$ , $g^{n}_{i,\alpha}=e^{\psi^{n}_{i,\alpha}}$ and

G_{i,\alpha}=\rho_{i,\alpha}^{n+1}g_{i,\alpha}^{n},\quad\alpha\in\mathcal{A}.

Let $\beta$ be such that

G_{i,\beta}=\min_{\alpha\in\mathcal{A}}G_{i,\alpha},

it suffices to prove $G_{i,\beta}\geq 0$ . We discuss in cases:

(i) $K_{\beta}$ is an interior cell. On the cell $K_{\beta}$ we have

	$\displaystyle g_{i,\beta}^{n}G_{i,\beta}$	$\displaystyle=\sum_{j=1}^{d}\lambda_{j}[\bar{D}_{i,\beta+e_{j}/2}(G_{i,\beta+e_{j}}-G_{i,\beta})-\bar{D}_{i,\beta-e_{j}/2}(G_{i,\beta}-G_{i,\beta-e_{j}})]+\rho_{i,\beta}^{n}$
		$\displaystyle\geq\rho_{i,\beta}^{n},$

where we used the fact $G_{i,\beta}\leq G_{i,\beta\pm e_{j}}$ and $\bar{D}_{i,\beta\pm e_{j}/2}>0$ . Since $g_{i,\beta}^{n}>0,$ so $G_{i,\beta}\geq 0.$

(ii) $K_{\beta}$ is a boundary cell( $K_{\beta}\cap\partial\Omega_{D}\neq\emptyset$ ). We only deal with the case $\beta(1)=1,$ remaining cases are similar. In such case,

	$\displaystyle g_{i,\beta}^{n}G_{i,\beta}=$	$\displaystyle\sum_{j=2}^{d}\lambda_{j}[\bar{D}_{i,\beta+e_{j}/2}(G_{i,\beta+e_{j}}-G_{i,\beta})-\bar{D}_{i,\beta-e_{j}/2}(G_{i,\beta}-G_{i,\beta-e_{j}})]$
		$\displaystyle+\lambda_{1}\bar{D}_{i,\beta+e_{1}/2}(G_{i,\beta+e_{1}}-G_{i,\beta})$
		$\displaystyle-2\lambda_{1}D_{i}(x_{\beta-e_{1}/2})g_{i}^{b}(x_{\beta-e_{1}/2},t_{n})\bigg{(}G_{i,\beta}-\frac{\rho^{b}_{i}(x_{\beta-e_{1}/2},t_{n+1})}{g_{i}^{b}(x_{\beta-e_{1}/2},t_{n})}\bigg{)}+\rho_{i,\beta}^{n}.$

Due to $G_{i,\beta}\leq G_{i,\beta\pm e_{j}}$ and $\bar{D}_{i,\beta\pm e_{j}/2}\geq 0,$ we have

\bigg{(}g_{i,\beta}^{n}+2\lambda_{1}D_{i}(x_{\beta-e_{1}/2})g_{i}^{b}(x_{\beta-e_{1}/2},t_{n})\bigg{)}G_{i,\beta}\geq 2\lambda_{1}D_{i}(x_{\beta-e_{1}/2})\rho^{b}_{i}(x_{\beta-e_{1}/2},t_{n+1})+\rho_{i,\beta}^{n}\geq 0,

which with $g_{i,\beta}^{n}+2\lambda_{1}D_{i}(x_{\beta-e_{1}/2})g_{i}^{b}(x_{\beta-e_{1}/2},t_{n})>0$ ensures $G_{i,\beta}\geq 0.$

(iii) $K_{\beta}$ is a boundary cell ( $K_{\beta}\cap\partial\Omega_{N}\neq\emptyset$ ). Again we only deal with the case $\beta(l)=1$ . In such case,

	$\displaystyle g_{i,\beta}^{n}G_{i,\beta}=$	$\displaystyle\sum_{j=1,j\neq l}^{d}\lambda_{j}[\bar{D}_{i,\beta+e_{j}/2}(G_{i,\beta+e_{j}}-G_{i,\beta})-\bar{D}_{i,\beta-e_{j}/2}(G_{i,\beta}-G_{i,\beta-e_{j}})]$
		$\displaystyle+\lambda_{l}\bar{D}_{i,\beta+1/2e_{l}}(G_{i,\beta+e_{l}}-G_{i,\beta})+\rho_{i,\beta}^{n}$
	$\displaystyle\geq$	$\displaystyle\rho_{i,\beta}^{n}\geq 0.$

This also gives $G_{i,\beta}\geq 0.$ The proof is thus complete. ∎

3.3.4. Energy dissipation

If $\partial\Omega_{D}=\emptyset$ , then solutions $\rho^{n+1}_{\alpha}$ obtained by (3.7) are conservative and energy dissipating in addition to the non-negativity. Let a discrete version of the free energy (2.7) be defined as

E^{n}_{h}=\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|\left[\sum_{i=1}^{m}\rho^{n}_{i,\alpha}(\log\rho^{n}_{i,\alpha}-1)+\frac{1}{2k_{B}T}\left(f_{\alpha}+\sum_{i=1}^{m}q_{i}\rho^{n}_{i,\alpha}\right)\phi^{n}_{\alpha}+\frac{1}{k_{B}T}\sum_{i=1}^{m}\rho^{n}_{i,\alpha}\mu_{i,\alpha}\right],

(3.12)

we have the following result.

Theorem 3.4.

Let $\rho_{\alpha}^{n}$ be obtained from (3.7) by using either (3.8), (3.9), or (3.10) for $e^{-\psi^{n}_{i,\alpha+e_{i}/2}}$ . Let $\phi^{n}_{\alpha}$ be obtained from (3.11). If $\partial\Omega_{D}=\emptyset$ , then we have:
(i) Mass conservation:

\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|\rho_{i,\alpha}^{n+1}=\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|\rho_{i,\alpha}^{n}\quad\text{ for }n\geq 0,i=1,\cdots,m;

(ii) Energy dissipation: There exists $\tau^{*}>0$ such that if $\tau\in(0,\tau^{*})$ , then

E^{n+1}_{h}-E^{n}_{h}\leq-\frac{\tau}{2}I^{n},

(3.13)

where

I^{n}=\sum_{i=1}^{m}\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}|K_{\alpha}|\frac{C_{i,\alpha+e_{j}/2}}{h_{j}}\left(\log(\rho^{n+1}_{i,\alpha+e_{j}}e^{\psi^{n}_{\alpha+e_{j}}})-\log(\rho^{n+1}_{i,\alpha}e^{\psi^{n}_{\alpha}})\right)\geq 0.

If we let

\epsilon_{min}=\min_{x\in\bar{\Omega}}\epsilon(x),\ \epsilon_{max}=\max_{x\in\bar{\Omega}}\epsilon(x),\ D_{max}=\max_{i,x\in\bar{\Omega}}D_{i}(x),

then $\tau^{*}$ can be quantified by

\tau^{*}=\frac{k_{B}T\epsilon^{2}_{min}}{4\pi\epsilon_{max}D_{max}\max_{i,\alpha,n}\rho^{n}_{i,\alpha}\sum_{i=1}^{m}q_{i}^{2}}e^{-\max_{i,j,\alpha}|\psi^{n}_{i,\alpha+e_{j}}-\psi^{n}_{i,\alpha}|}.

Remark 3.2.

We remark that $\tau^{*}$ is of size $O(1)$ , though it appears to be dependent on numerical solutions. For $h_{j}$ small, the exponential term is only of size $e^{O(h)}$ , therefore bounded. As $n$ increases, the solution $\{\rho^{n}_{\alpha}\}$ is expected to converge to the steady-state and therefore bounded from above, hence we simply use the notation $\max_{i,\alpha,n}\rho^{n}_{i,\alpha}$ . The boundedness of $\rho^{n}$ in $n$ for some cases has been established in Theorem 3.2 for the corresponding semi-discrete scheme.

The proof is deferred to Appendix B.

3.3.5. Preservation of steady-states

With no-flux boundary conditions, scheme (3.7) can be shown to be steady-state preserving. Based on the discussion in section 2.4, we say a discrete function $\rho_{\alpha}$ is at steady-state if

\rho_{i,\alpha}=c_{i}e^{-\frac{1}{k_{B}T}(q_{i}\phi_{\alpha}+\mu_{i,\alpha})},\quad i=1,\cdots,m,\quad\alpha\in\mathcal{A},

(3.14)

where $\phi_{\alpha}$ satisfies (3.11) with $\rho_{i,\alpha}$ replaced by the above relation, which is a nonlinear algebraic equation for $\phi_{\alpha}$ uniquely determined for each $(c_{1},\cdots c_{m})$ . We have the following theorem.

Theorem 3.5.

Let the assumptions in Theorem 3.4 be met, then

(i) If $\rho^{0}_{\alpha}$ is already at steady-state, then $\rho^{n}_{\alpha}=\rho^{0}_{\alpha}$ for $n\geq 1$ .

(ii) If $E^{n+1}_{h}=E^{n}_{h}$ , then $\rho^{n}_{\alpha}$ must be at steady-state.

(iii) If $\rho^{n}_{i,\alpha},\phi^{n}_{\alpha}$ converge as $n\to\infty$ , then their limits are determined by

\rho^{\infty}_{i,\alpha}=c^{\infty}_{i}e^{-\frac{1}{k_{B}T}(q_{i}\phi^{\infty}_{\alpha}+\mu_{i,\alpha})},\quad c^{\infty}_{i}=\frac{\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|\rho^{0}_{i,\alpha}}{\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|e^{-\frac{1}{k_{B}T}(q_{i}\phi^{\infty}_{\alpha}+\mu_{i,\alpha})}},

where $\phi^{\infty}_{\alpha}$ is obtained by solving (3.11) by using $\rho^{\infty}_{i,\alpha}$ .

Proof.

(i) We only need to prove $\rho^{1}_{i,\alpha}=\rho^{0}_{i,\alpha},$ for all $i=1,\cdots,m,\ \alpha\in\mathcal{A}$ . Summing (3.7) with $n=0$ against $|K_{\alpha}|\rho^{1}_{i,\alpha}/\rho^{0}_{i,\alpha}$ , using summation by parts, we obtain

	$\displaystyle\sum_{\alpha\in\mathcal{A}}\|K_{\alpha}\|(\rho^{1}_{i,\alpha}-\rho^{0}_{i,\alpha})\frac{\rho^{1}_{i,\alpha}}{\rho^{0}_{i,\alpha}}=$	$\displaystyle\tau\sum_{j=1}^{d}\sum_{\alpha\in\mathcal{A}}\|K_{\alpha}\|\frac{1}{h_{j}}(C_{i,\alpha+e_{j}/2}-C_{i,\alpha-e_{j}/2})\frac{\rho^{1}_{i,\alpha}}{\rho^{0}_{i,\alpha}}$		(3.15)
	$\displaystyle=$	$\displaystyle-\tau\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\|K_{\alpha}\|\frac{1}{h_{j}}C_{i,\alpha+e_{j}/2}\left(\frac{\rho^{1}_{i,\alpha+e_{j}}}{\rho^{0}_{i,\alpha+e_{j}}}-\frac{\rho^{1}_{i,\alpha}}{\rho^{0}_{i,\alpha}}\right).$		(3.15)

Substituting $\rho^{0}_{i,\alpha}=c_{i}e^{-\psi^{0}_{i,\alpha}}$ into $C_{i,\alpha+e_{j}/2}$ , the right hand side of (3.15) becomes

\displaystyle RHS=

\displaystyle-\tau c_{i}\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}|K_{\alpha}|\frac{D_{i,\alpha+e_{j}/2}e^{-\psi^{0}_{i,\alpha+e_{j}/2}}}{h^{2}_{j}}\left(\frac{\rho^{1}_{i,\alpha+e_{j}}}{\rho^{0}_{i,\alpha+e_{j}}}-\frac{\rho^{1}_{i,\alpha}}{\rho^{0}_{i,\alpha}}\right)^{2}\leq 0.

Adding $\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|(\rho^{0}_{i,\alpha}-\rho^{1}_{i,\alpha})=0$ to the left hand side of (3.15) leads to

	$\displaystyle LHS=$	$\displaystyle\sum_{\alpha\in\mathcal{A}}\|K_{\alpha}\|\left[(\rho^{1}_{i,\alpha}-\rho^{0}_{i,\alpha})\frac{\rho^{1}_{i,\alpha}}{\rho^{0}_{i,\alpha}}+(\rho^{0}_{i,\alpha}-\rho^{1}_{i,\alpha})\right]$
	$\displaystyle=$	$\displaystyle\sum_{\alpha\in\mathcal{A}}\|K_{\alpha}\|\frac{(\rho^{1}_{i,\alpha}-\rho^{0}_{i,\alpha})^{2}}{\rho^{0}_{i,\alpha}}\geq 0.$

Hence $LHS=RHS\equiv 0$ , we must have

\rho^{1}_{i,\alpha}=\rho^{0}_{i,\alpha},\quad i=1,\cdots,m,\quad\alpha\in\mathcal{A}.

(ii) The inequality (3.13) when combined with $E^{n+1}_{h}=E^{n}_{h}$ leads to $I^{n}=0$ . From the proof of Theorem 3.4 in Appendix B it follows

\rho^{n+1}_{i,\alpha}=\rho^{n}_{i,\alpha}.

(iii) Since $E^{n}_{h}$ is non-increasing in $n$ , and we can verify that $E_{h}^{n}$ is bounded from below, hence

\lim_{n\to\infty}E^{n}_{h}=\inf\{E^{n}_{h}\}.

Taking the limit in (3.13), we have $\lim_{n\to\infty}I^{n}=0$ , which implies

\rho^{\infty}_{i,\alpha}=c^{\infty}_{i}e^{-\psi^{\infty}_{i,\alpha}}.

Conservation of mass gives

c^{\infty}_{i}=\frac{\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|\rho^{0}_{i,\alpha}}{\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|e^{-\psi^{\infty}_{i,\alpha}}}.\quad i=1,\cdots,m,\quad\alpha\in\mathcal{A},

where $\phi^{\infty}_{\alpha}$ in $\psi^{\infty}_{i,\alpha}=\frac{1}{k_{B}T}(q_{i}\phi^{\infty}_{\alpha}+\mu_{i,\alpha})$ is obtained by solving (3.11) using $\rho^{\infty}_{i,\alpha}$ . ∎

4. Second order in time discretization

The semi-discrete scheme (3.2a) is first order accurate, one can design higher order in time schemes based on (3.1).

The following is a second order time discretization,

\frac{\rho^{n+1}_{i}-\rho^{n}_{i}}{\tau}=R[(\rho^{n+1}_{i}+\rho^{n}_{i})/2,\frac{3}{2}\psi^{n}_{i}-\frac{1}{2}\psi^{n-1}_{i}].

This can be expressed as a prediction-correction method,

\frac{\rho^{*}_{i}-\rho^{n}_{i}}{\tau/2}=R[\rho^{*}_{i},\frac{3}{2}\psi^{n}_{i}-\frac{1}{2}\psi^{n-1}_{i}],\quad\rho^{n+1}_{i}=2\rho^{*}_{i}-\rho^{n}_{i}.

(4.1)

As argued for the first order scheme, this scheme is well-defined.

4.1. Second order fully-discrete scheme

With central spatial difference, our fully discrete second order (in both space and time) scheme reads


		$\displaystyle\frac{\rho^{}_{i,\alpha}-\rho^{n}_{i,\alpha}}{\tau/2}=R_{\alpha}[\rho^{}_{i},\frac{3}{2}\psi^{n}_{i}-\frac{1}{2}\psi^{n-1}_{i}],$		(4.2a)
		$\displaystyle\rho^{n+1}_{i,\alpha}=2\rho^{*}_{i,\alpha}-\rho^{n}_{i,\alpha}.$		(4.2b)

Positivity of $\rho^{n+1}_{\alpha}$ can be ensured if time steps are sufficient small.

Theorem 4.1.

Let $\rho^{n+1}_{\alpha}$ be obtained from (4.2). If $\rho^{n}_{\alpha}\geq 0$ for all $\alpha\in\mathcal{A}$ , and $\rho^{b}(x,t)\geq 0$ for $x\in\partial\Omega_{D}$ , then

\rho^{n+1}_{\alpha}\geq 0,\quad\alpha\in\mathcal{A}

provided $\tau$ is sufficiently small.

Proof.

Inserting (4.2b) into (4.2a) leads to the following compact form of the scheme (4.2):

\rho^{n+1}_{i,\alpha}-\frac{\tau}{2}R_{\alpha}[\rho^{n+1}_{i},\frac{3}{2}\psi^{n}_{i}-\frac{1}{2}\psi^{n-1}_{i}]=\rho^{n}_{i,\alpha}+\frac{\tau}{2}R_{\alpha}[\rho^{n}_{i},\frac{3}{2}\psi^{n}_{i}-\frac{1}{2}\psi^{n-1}_{i}],

(4.3)

where we have used the linearity of $R_{\alpha}[\cdot,\cdot]$ on the first entry.

Set

g^{*}_{i,\alpha}=e^{\frac{3}{2}\psi^{n}_{i,\alpha}-\frac{1}{2}\psi^{n-1}_{i,\alpha}},\ \bar{D}^{*}_{i,\alpha+e_{j}/2}=D_{i,\alpha+e_{j}/2}e^{-\frac{3}{2}\psi^{n}_{i,\alpha+e_{j}/2}+\frac{1}{2}\psi^{n-1}_{i,\alpha+e_{j}/2}},\ G^{n}_{i,\alpha}=\rho^{n}_{i,\alpha}g^{*}_{i,\alpha},

then the scheme (4.3) can be rewritten as

		$\displaystyle g^{}_{i,\alpha}G^{n+1}_{i,\alpha}-\sum_{j=1}^{d}\frac{\tau}{h_{j}^{2}}[\bar{D}^{}_{i,\alpha+e_{j}/2}(G^{n+1}_{i,\alpha+e_{j}}-G^{n+1}_{i,\alpha})-\bar{D}^{*}_{i,\alpha-e_{j}/2}(G^{n+1}_{i,\alpha}-G^{n+1}_{i,\alpha-e_{j}})]$		(4.4)
		$\displaystyle=g^{}_{i,\alpha}G^{n}_{i,\alpha}+\sum_{j=1}^{d}\frac{\tau}{h_{j}^{2}}[\bar{D}^{}_{i,\alpha+e_{j}/2}(G^{n}_{i,\alpha+e_{j}}-G^{n}_{i,\alpha})-\bar{D}^{*}_{i,\alpha-e_{j}/2}(G^{n}_{i,\alpha}-G^{n}_{i,\alpha-e_{j}})].$		(4.4)

Let $\beta$ be such that

G^{n+1}_{i,\beta}=\min_{\alpha\in\mathcal{A}}G^{n+1}_{i,\alpha},

it suffices to prove $G^{n+1}_{i,\beta}\geq 0$ . We prove the result when $K_{\beta}$ is an interior cell, the result for boundary cells can be proved similarly.

Since $G^{n+1}_{i,\beta}\leq G^{n+1}_{i,\beta\pm e_{j}}$ and $G^{n}_{i,\beta\pm j}\geq 0$ , thus equation (4.4) on cell $K_{\beta}$ reduces to the inequality:

g^{*}_{i,\beta}G^{n+1}_{i,\beta}\geq\left(g^{*}_{i,\beta}-\tau\sum_{j=1}^{d}\frac{1}{h_{j}^{2}}(\bar{D}^{*}_{i,\beta+e_{j}/2}+\bar{D}^{*}_{i,\beta-e_{j}/2})\right)G^{n}_{i,\beta},

we see that $G^{n+1}_{i,\beta}\geq 0$ is insured if

\tau\leq\min_{\alpha}\left\{\frac{g^{*}_{i,\alpha}}{\sum_{j=1}^{d}\frac{1}{h_{j}^{2}}(\bar{D}^{*}_{i,\alpha+e_{j}/2}+\bar{D}^{*}_{i,\alpha-e_{j}/2})}\right\}.

The stated result thus follows. ∎

We should point out that numerical density $\{\rho^{n}_{\alpha}\}$ obtained by the second order scheme (4.2) may not be non-negative for large time step $\tau$ , though $\{\rho^{*}_{\alpha}\}$ stays positive. We shall restore solution positivity by using a local limiter, which was first introduced in [23] for one-dimensional case.

4.2. Positivity-preserving limiter

We present a local limiter to restore positivity of $\rho$ if

\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|\rho_{\alpha}>0,

but $\rho_{\beta}<0$ for some $\beta\in\mathcal{A}$ . The idea is to find a neighboring index set $S_{\beta}$ such that the local average

\bar{\rho}_{\beta}=\frac{1}{|S_{\beta}|}\sum_{\gamma\in S_{\beta}}|K_{\gamma}|\rho_{\gamma}>0,

where $|S_{\beta}|$ denotes the minimum number of indices for which $\rho_{\gamma}\neq 0$ and $\bar{\rho}_{\beta}>0$ , then use this local average as a reference to define the following scaling limiter

\tilde{\rho}_{\alpha}=\theta\rho_{\alpha}+(1-\theta)\bar{\rho}_{\beta}/|K_{\alpha}|,\quad\alpha\in S_{\beta},

(4.5)

where

\theta=\min\left\{1,\frac{\bar{\rho_{\beta}}}{\bar{\rho}_{\beta}-\rho_{\min}}\right\},\quad\rho_{\min}=\min_{\gamma\in S_{\beta}}|K_{\gamma}|\rho_{\gamma}.

Recall the result stated in Lemma 5.1 in [24], such limiter restores solution positivity and respects the local mass conservation. In addition, for any sequence $g_{\alpha}$ with $g_{\alpha}\geq 0$ , we have

\displaystyle|\tilde{\rho}_{\alpha}-g_{\alpha}|\leq(1+|S_{\beta}|\Lambda)\max_{\gamma\in S_{\beta}}|\rho_{\gamma}-g_{\gamma}|,\quad\alpha\in S_{\beta},

(4.6)

where $\Lambda$ is the upper bound of mesh ratio $|K_{\gamma}|/|K_{\alpha}|$ . Let $\rho_{\alpha}$ be the approximation of $\rho(x)\geq 0$ , we let $g_{\alpha}=\rho(x_{\alpha})$ or the average of $\rho$ on $K_{\alpha}$ , so we can assert that the accuracy is not destroyed by the limiter as long as $|S_{\beta}|\Lambda$ is uniformly bounded. Boundedness of $|S_{\beta}|$ for shape-regular meshes was rigorously proved in [24] for the one-dimensional case. We restate such result in the present setting in the following.

Theorem 4.2.

Let $\{\rho_{\alpha}\}$ be an approximation of $\rho(x)\geq 0$ over shape regular meshes, and $\rho\in C^{k}(\Omega)$ ( $k\geq 2$ ). If $\rho_{\beta}<0$ (or only finite number of neighboring values are negative), then there exists $C^{*}>0$ finite such that

|\tilde{\rho}_{\alpha}-\rho(x_{\alpha})|\leq C^{*}\max_{\alpha\in S_{\beta}}|\rho_{\alpha}-\rho(x_{\alpha})|,\quad\forall\alpha\in S_{\beta},

where $C^{*}$ may depend on the local meshes associated with $S_{\beta}$ .

Proof.

For simplicity, we prove only for the case of uniform meshes (e.g. uniform in each dimension). Let $h=\min_{1\leq j\leq d}h_{j}\leq 1$ and $h_{j}\leq\Lambda h$ for some $\Lambda>0$ . From (4.6) we see that it suffices to show there exists $A^{*}>0$ finite such that $|S_{\beta}|\leq A^{*},$ with which we will have $C^{*}=1+A^{*}\Lambda$ . Under the smoothness assumption of $\rho$ we may assume $|\rho_{\alpha}-\rho(x_{\alpha})|\leq Ch^{k}$ . Under the assumption $\rho_{\beta}<0$ , $\rho$ must touch zero near $x_{\beta}$ . We discuss the case where $\rho(x^{*})=0$ and $\nabla\rho(x^{*})=\vec{0}$ with $\rho(x)>0$ for $x(j)\geq x^{*}(j)$ , $j=1,\cdots,d,$ locally with $x^{*}\in K_{\beta}.$ To be concrete, we consider $\beta=(1,\cdots,1)$ and $\int_{K_{\beta}}\rho(x)dx>0$ . From the limiter construction we have $S_{\beta}$ such that

\sum_{\alpha\in S_{\beta}}|K_{\alpha}|\rho_{\alpha}>0.

(4.7)

The rest of the proof is devoted to bounding $|S_{\beta}|$ . The assumed error bound gives

\rho_{\alpha}\geq\rho(x_{\alpha})-Ch^{k}.

(4.8)

From $\rho\in C^{k}(\Omega)(k\geq 2)$ , we have

\rho(x_{\alpha})\geq\bar{\rho}_{\alpha}-\lambda\Lambda^{2}h^{2},

(4.9)

with $\lambda=\frac{d}{24}\max_{j=1,\cdots,d}|\partial_{x_{j}x_{j}}\rho|$ and the cell average $\bar{\rho}_{\alpha}=\frac{1}{|K_{\alpha}|}\int_{K_{\alpha}}\rho(x)dx$ . From (4.8) and (4.9), we see that the left hand side of (4.7) is bounded from below by

	$\displaystyle\sum_{\alpha\in S_{\beta}}\|K_{\alpha}\|\rho_{\alpha}\geq$	$\displaystyle\sum_{\alpha\in S_{\beta}}\|K_{\alpha}\|\left(\bar{\rho}_{\alpha}-(C+\lambda\Lambda^{2})h^{2}\right)$		(4.10)
	$\displaystyle=$	$\displaystyle\int_{\cup_{\alpha\in S_{\beta}}K_{\alpha}}\rho(x)dx-(C+\lambda\Lambda^{2})h^{2}\sum_{\alpha\in S_{\beta}}\|K_{\alpha}\|.$		(4.10)

Without loss of generality we assume $\cup_{\alpha\in S_{\beta}}K_{\alpha}$ is a rectangle in $\mathbb{R}^{d}$ ; otherwise we could add more cells to complete the rectangle. Let

|\cup_{\alpha\in S_{\beta}}K_{\alpha}|=\Pi_{j=1}^{d}\eta_{j},\quad h_{j}\leq\eta_{j}\leq L_{j},

and $\vec{\eta}=(\eta_{1},\cdots,\eta_{d}),\ \vec{h}=(h_{1},\cdots,h_{d}).$ Rewriting integral in (4.10) we have

\sum_{\alpha\in S_{\beta}}|K_{\alpha}|\rho_{\alpha}\geq\left[g(\eta)-(C+\lambda\Lambda^{2})h^{2}\right]\sum_{\alpha\in S_{\beta}}|K_{\alpha}|,

where

g(\eta):=\int_{0}^{1}\!\cdots\!\int_{0}^{1}\!\!\rho\left({\rm diag}(\vec{\theta})\vec{\eta}+x_{\beta}-\frac{1}{2}\vec{h}\right)d\theta_{1}\cdots d\theta_{d}.

From the fact $h^{d}\leq\frac{\eta_{1}\cdots\eta_{d}}{|S_{\beta}|},$ we can see that the term in the bracket is bounded from below by

g(\eta)-(C+\lambda\Lambda^{2})\left(\frac{\eta_{1}\cdots\eta_{d}}{|S_{\beta}|}\right)^{2/d},

which is positive if

|S_{\beta}|>(C+\lambda\Lambda^{2})^{d/2}{g(\eta)}^{-d/2}\eta_{1}\cdots\eta_{d}.

This can be insured if we take

|S_{\beta}|=\lfloor A\rfloor+1,

where

A=(C+\lambda\Lambda^{2})^{d/2}\max_{\eta_{j}\in[h_{j},L_{j}],j=1,\cdots,d}{g(\eta)}^{-d/2}\eta_{1}\cdots\eta_{d}.

This is bounded and may depend on the local mesh of $K_{\beta}$ . ∎

Note that our numerical solutions feature the following property: if $\rho^{n}_{i,\alpha}=0$ , then

\rho^{n+1}_{i,\alpha}=2\rho^{*}_{i,\alpha}-\rho^{n}_{i,\alpha}\geq 0

due to the fact $\rho^{*}_{i,\alpha}\geq 0$ . This means that if $\rho^{\rm in}(x)=0$ on an interval, then $\rho^{1}_{i,\alpha}$ cannot be negative in most of nearby cells. Thus negative values appear only where the exact solution turns from zero to a positive value, and the number of these values are finitely many. Our result in Theorem 4.2 is thus applicable.

4.3. Algorithm

The following algorithm is only for the second order scheme with limiter.

(1)

Initialization: From the initial data $\rho^{in}_{i}(x)$ , obtain

$\rho^{0}_{i,\alpha}=\frac{1}{|K_{\alpha}|}\int_{K_{\alpha}}\rho_{i}^{in}(x)dx,\ i=1,\cdots,m,\ \alpha\in\mathcal{A},$

by using central point quadrature.
(2)

Update to get $\{\rho^{1}_{i,\alpha}\}$ : Compute $\{\phi^{0}_{\alpha}\}$ from (3.11), then obtain $\{\rho_{i,\alpha}^{1}\}$ by the first order scheme (3.7).
(3)

Update from $\{\rho^{n}_{i,\alpha}\}$ : For $n\geq 1$ , compute $\{\phi_{\alpha}^{n}\}$ from scheme (3.11) then get $\{\rho_{i,\alpha}^{n+1}\}$ from (4.2).
(4)

Reconstruction: if necessary, locally replace $\rho_{i,\alpha}^{n+1}$ by $\tilde{\rho}_{i,\alpha}^{n+1}$ using the limiter defined in (4.5).

The following algorithm can be called to find an admissible set $S_{\alpha}$ used in (4.5).

(i)

Start with $S_{\beta}=\{\beta\}$ , $p=1$ .
(ii)

For $l_{j}=\max\{1,\alpha(j)-p\}:\min\{\alpha(j)+p,N_{j}\}$ with $j=1,\cdots,d$ .
If $\alpha:=(l_{1},\cdots,l_{d})\notin S_{\beta}$ and $\rho^{n+1}_{i,\alpha}\neq 0$ , then set $S_{\beta}=S_{\beta}\cup\left\{\alpha\right\}$ .
If $\bar{\rho}_{\beta}>0$ , then stop, else go to (iii).
(iii)

Set $p=p+1$ and go to (ii).

Remark 4.1.

The coefficient matrices of the linear systems obtained by (3.7), (3.11), and (4.2a) are sparse, diagonally dominant, and symmetric, hence more efficient linear system solvers, such as the ILU preconditioner + FGMRES (see e.g., [37]), ILU preconditioner + Bicgstab (see e.g., [4]), can be used.

5. Numerical tests

In this section, we implement the fully discrete schemes (3.7) and (4.2) to demonstrate their orders of convergence and capacity to preserve solution properties. In both schemes the numerical solution $\phi^{n}_{\alpha}$ is computed by the scheme (3.11). Errors in the accuracy tests are measured in the following discrete $l^{1}$ norm:

\text{error}=\sum_{\alpha\in\mathcal{A}}|K_{\alpha}||\tilde{g}_{\alpha}-g_{\alpha}|.

Here $g_{\alpha}$ denotes the numerical solution, say $g_{\alpha}=\rho_{i,\alpha}^{n}$ or $\phi_{\alpha}^{n}$ at time $t=n\tau$ , and $\tilde{g}_{\alpha}$ indicates the cell average of the corresponding exact solutions.

In our numerical tests, the sparse linear systems obtained by (3.7), (3.11), and (4.2a) are solved by ILU preconditioned FGMRES [37] algorithm using compressed row format of the coefficient matrices. In the three-dimensional case, the coefficient matrices of the linear systems are 7-diagonal matrices. It is worth to mention that the compressed row format allows us to store a $l\times l$ 7-diagonal matrix by using at most $15l$ storage locations with $l=N_{x}\times N_{y}\times N_{z}$ . With $30\times 30\times 30$ cells, we can save $99\%$ of the storage space needed for storing the resulting coefficient matrices.

In our three examples below we consider the computational domain

\Omega=(0,1)\times(0,1)\times(0,1).

Example 5.1.

(Accuracy test) In this test we numerically verify the accuracy and order of schemes (3.7) and (4.2) by using manufactured solutions. Consider

\left\{\begin{array}[]{rl}\hfill\rho_{1}(\mathbf{x},t)&=4(x^{2}(1-x)^{2}+y(1-y))e^{-t},\\ \hfill\rho_{2}(\mathbf{x},t)&=(y(1-y)+z^{2}(1-z)^{2})e^{-t},\\ \hfill\phi(\mathbf{x},t)&=(x^{2}(1-x)^{2}+y(1-y)+z^{2}(1-z)^{2})e^{-t}\end{array}\right.

(5.1)

and

\partial\Omega_{D}=\{\mathbf{x}\in\bar{\Omega}:y=0,1\},\quad\partial\Omega_{N}=\partial\bar{\Omega}\setminus\partial\Omega_{D},

then they are exact solutions to the following problem

\left\{\begin{array}[]{rl}\hfill\partial_{t}\rho_{1}=&\nabla\cdot(\nabla\rho_{1}+\rho_{1}\nabla\phi)+f_{1}(\mathbf{x},t),\hfill\ \ \ \mathbf{x}\in\Omega,\ t>0,\\ \hfill\partial_{t}\rho_{2}=&\nabla\cdot(\nabla\rho_{2}-\rho_{2}\nabla\phi)+f_{2}(\mathbf{x},t),\hfill\ \ \ \mathbf{x}\in\Omega,\ t>0,\\ \hfill-\Delta\psi=&\rho_{1}-\rho_{2}+f_{3}(\mathbf{x},t),\hfill\ \ \ \mathbf{x}\in\Omega,\ t>0,\end{array}\right.

(5.2)

where source terms $f_{1}(\mathbf{x},t),f_{2}(\mathbf{x},t)$ and $f_{3}(\mathbf{x},t),$ and the initial and boundary conditions are determined by the exact solutions.

We first test the accuracy of the semi-implicit scheme (3.7) by using various spatial step size $h$ , errors and orders at $t=1$ are listed in Table 1 (with $\tau=h$ ) and in Table 2 (with $\tau=h^{2}$ ), respectively. We observe the first order accuracy in time and the second order accuracy in space. We then test the accuracy of the scheme (4.2) with time step size $\tau=h$ . From Table 3, we see the second order accuracy in both time and space.

Table 1. Scheme (3.7) with

\tau=h

$N_{x}\times N_{y}\times N_{z}$	$\rho_{1}$ error	order	$\rho_{2}$ error	order	$\phi$ error	order
$8\times 8\times 8$	4.7508E-02	-	1.3904E-02	-	5.7213E-03	-
$16\times 16\times 16$	2.1283E-02	1.1585	5.8701E-03	1.2440	2.0987E-03	1.4468
$32\times 32\times 32$	1.0060E-02	1.0811	2.6956E-03	1.1228	8.6460E-04	1.2794
$64\times 64\times 64$	4.8890E-03	1.0410	1.2915E-03	1.0616	3.8667E-04	1.1609

Table 2. Scheme (3.7) with

\tau=h^{2}

$N_{x}\times N_{y}\times N_{z}$	$\rho_{1}$ error	order	$\rho_{2}$ error	order	$\phi$ error	order
$8\times 8\times 8$	1.1252E-02	-	4.0301E-03	-	3.1194E-03	-
$16\times 16\times 16$	2.7824E-03	2.0158	9.8548E-04	2.0319	7.7117E-04	2.0161
$32\times 32\times 32$	6.9369E-04	2.0040	2.4502E-04	2.0079	1.9225E-04	2.0041
$64\times 64\times 64$	1.7330E-04	2.0010	6.1170E-05	2.0020	4.8028E-05	2.0010

Table 3. Scheme (4.2) with

\tau=h

$N_{x}\times N_{y}\times N_{z}$	$\rho_{1}$ error	order	$\rho_{2}$ error	order	$\phi$ error	order
$8\times 8\times 8$	5.5476E-03	-	2.3247E-03	-	2.7378E-03	-
$16\times 16\times 16$	1.5073E-03	1.8799	6.0465E-04	1.9429	6.7758E-04	2.0146
$32\times 32\times 32$	3.9635E-04	1.9271	1.5851E-04	1.9315	1.6895E-04	2.0038
$64\times 64\times 64$	1.0182E-04	1.9608	4.0875E-05	1.9553	4.2206E-05	2.0011

Example 5.2.

(Solution positivity) We consider the two-species PNP system with initial data of form

\left\{\begin{array}[]{rl}\hfill\partial_{t}\rho_{1}=&\nabla\cdot(\nabla\rho_{1}+\rho_{1}\nabla\phi),\hfill\ \ \ \mathbf{x}\in\Omega,\ t>0,\\ \hfill\partial_{t}\rho_{2}=&\nabla\cdot(\nabla\rho_{2}-\rho_{2}\nabla\phi),\hfill\ \ \ \mathbf{x}\in\Omega,\ t>0,\\ \hfill-\Delta\psi=&\rho_{1}-\rho_{2}+10\chi_{{}_{[0.2,0.4]\times[0.2,0.4]\times[0.2,0.4]}},\hfill\ \ \ \mathbf{x}\in\Omega,\ t>0,\\ \hfill\rho_{1}^{in}(\mathbf{x})=&\chi_{{}_{[0,0.25]\times[0,0.25]\times[0,0.25]}},\\ \hfill\rho_{2}^{in}(\mathbf{x})=&2\chi_{{}_{[0,0.25]\times[0,0.25]\times[0,0.25]}}.\end{array}\right.

(5.3)

This corresponds to (1.1) with $D_{1}=D_{2}=1,$ $q_{1}=-q_{2}=1,$ $k_{B}T=1,$ $\epsilon(\mathbf{x})=4\pi$ , $\mu_{i}=0$ , and $f(\mathbf{x})=10\chi_{{}_{[0.2,0.4]\times[0.2,0.4]\times[0.2,0.4]}}.$

With $\partial\Omega_{D}=\{\mathbf{x}\in\bar{\Omega}:y=0,1\},$ and $\partial\Omega_{N}=\partial\bar{\Omega}\setminus\partial\Omega_{D}$ , we solve the problem subject to mixed boundary conditions

\left\{\begin{array}[]{rl}\hfill(\nabla\phi)\cdot\mathbf{n}=&0,\ (\nabla\rho_{1}+\rho_{1}\nabla\phi)\cdot\mathbf{n}=0,\ (\nabla\rho_{2}-\rho_{2}\nabla\phi)\cdot\mathbf{n}=0,\hfill\ \ \ \mathbf{x}\in\partial\Omega_{N},\\ \hfill\phi^{b}(\mathbf{x},t)=&(x^{2}(1-x)^{2}+z^{2}(1-z)^{2})e^{-t},\hfill\ \ \ \mathbf{x}\in\partial\Omega_{D},\\ \hfill\rho^{b}_{1}(\mathbf{x},t)=&4x^{2}(1-x)^{2}e^{-t},\hfill\ \ \ \mathbf{x}\in\partial\Omega_{D},\\ \hfill\rho^{b}_{2}(\mathbf{x},t)=&z^{2}(1-z)^{2}e^{-t},\hfill\ \ \ \mathbf{x}\in\partial\Omega_{D}.\end{array}\right.

(5.4)

We use $30\times 30\times 30$ cells with $\tau=0.5h$ to compute numerical solutions up to $t=2$ . Given in Fig.1 are the time evolution of numerical solutions (top three rows) and the minimum of $\rho_{1},\rho_{2}$ (bottom row) obtained by the scheme (3.7), showing non-negative approximations for both $\rho_{1}$ and $\rho_{2}$ . Results obtained by the scheme (4.2) are given in Fig.2. Note that the positivity preserving limiter keeps being invoked when we use the scheme (4.2). The CPU time (average of 10 simulations) needed for running the schemes (3.7) and (4.2) are $207.27$ seconds and $203.15$ seconds, respectively, from which we see that the second-order scheme is as efficient as the first order scheme.

Refer to caption — Figure 1. Example 5.2: $\rho_{1},\rho_{2},\phi$ computed by scheme (3.7)

Example 5.3.

(Mass conservation and energy dissipation) In this numerical example we test both mass conservation and energy dissipation properties of our schemes.

We consider system (5.3) with zero flux boundary conditions:

(\nabla\phi)\cdot\mathbf{n}=0,\ (\nabla\rho_{1}+\rho_{1}\nabla\phi)\cdot\mathbf{n}=0,\ (\nabla\rho_{2}-\rho_{2}\nabla\phi)\cdot\mathbf{n}=0,\hfill\ \ \ x\in\partial\Omega.

Numerical approximations of $\rho_{1}$ and $\rho_{2}$ at $t=2$ obtained by the scheme (3.7) are given in Fig.3. We can see by comparing Fig.3 and Fig.1 that boundary conditions have strong effects on the solution profiles. In Fig.4 (left) are the time evolution of the total mass and free energy obtained by the scheme (3.7), the results verify our theoretical findings in Theorem 3.4. In Fig.4 (right) are plots of the free energy and total mass obtained by (4.2). In this test the second order scheme looks also energy dissipative and mass conservative.

6. Concluding remarks

In this paper, we have developed unconditional structure-preserving schemes for PNP equations in more general settings. These schemes are shown to preserve several important physical laws at the fully discrete level including: mass conservation, solution positivity, and free energy dissipation. The non-logarithmic Landau reformulation of the model is important, enabling us to construct a simple, easy-to-implement fully discrete scheme (first order in time, second order in space), which proved to satisfy all three desired properties of the continuous model with only $O(1)$ time step restriction. We further designed a second order (in both time and space) scheme, which has the same computational complexity as the first order scheme. For such second order scheme, we employed a local scaling limiter to restore solution positivity where necessary. Moreover, we rigorously proved that the limiter does not destroy the desired accuracy. Three-dimensional numerical tests are conducted to evaluate the scheme performance and verify our theoretical findings. Our schemes presented with $\mu_{i}$ given can be applied to complex chemical potentials without difficulty.

Acknowledgments

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

Appendix A Proof of Theorem 3.1

Proof.

The elliptic problem (3.2) can be rewritten in $w=\rho^{n+1}_{i}e^{\psi^{n}_{i}}$ as


		$\displaystyle e^{-\psi^{n}_{i}}w-\tau\nabla\cdot(D_{i}(x)e^{-\psi^{n}_{i}}\nabla w)=\rho^{n}_{i},$		(A.1a)
		$\displaystyle w=\rho^{b}_{i}(x,t_{n+1})e^{\psi^{b}_{i}(x,t_{n})},\quad x\in\partial\Omega_{D},$		(A.1b)
		$\displaystyle(\nabla w)\cdot\mathbf{n}=0,\quad x\in\partial\Omega_{N}.$		(A.1c)

Let $\gamma_{0}$ be the trace operator on $\partial\Omega_{D}$ . The above problem admits a variational formulation of form

\displaystyle B[u,v]=Lv,\quad u,v\in H,

(A.2)

where for a Dirichlet lift $G\in H^{2}(\Omega)$ with trace $\gamma_{0}(G)=\rho^{b}_{i}(x,t_{n+1})e^{\psi^{b}_{i}(x,t_{n})}$ , we find

w=u+G.

Here


		$\displaystyle H=\{v\in H^{1}(\Omega):\gamma_{0}(v)=0\text{ on }\partial\Omega_{D}\},$		(A.3a)
		$\displaystyle B[u,v]=\int_{\Omega}(\tau D_{i}(x)e^{-\psi^{n}_{i}}\nabla u\cdot\nabla v+e^{-\psi^{n}_{i}}uv)dx,$		(A.3b)
		$\displaystyle Lv=\int_{\Omega}(\rho^{n}_{i}-e^{-\psi^{n}_{i}}G)v-\tau D_{i}(x)e^{-\psi^{n}_{i}}\nabla G\cdot\nabla vdx.$		(A.3c)

Under the assumptions, the celebrated Lax-Milgram theorem ([13] Theorem 5.8) ensures that the variational problem (A.2) admits a unique solution $u\in H$ . We thus obtain

\rho^{n+1}_{i}=e^{-\psi^{n}_{i}}(u+G).

Regularity for $\rho_{i}^{n+1}$ follows from the classical elliptic regularity for $u$ .

Similarly, the variational problem for (3.3) can also be written as (A.2) with


		$\displaystyle B[u,v]=\int_{\Omega}\epsilon(x)\nabla u\cdot\nabla vdx,$		(A.4a)
		$\displaystyle Lv=\int_{\Omega}4\pi\left(f(x)+\sum_{i=1}^{m}q_{i}\rho^{n+1}_{i}\right)v-\epsilon(x)\nabla G\cdot\nabla vdx,$		(A.4b)

where the Dirichlet lift $G\in H^{2}(\Omega)$ with $\gamma_{0}(G)=\phi^{b}(x,t_{n+1})\text{ on }x\in\partial\Omega_{D}.$ Here one can use the Poincaré-Friedrichs’ inequality of form $\|u\|_{L^{2}}\leq C_{F}\|\nabla u\|_{L^{2}}$ , which holds if $u=0$ on a set of $\partial\Omega$ with non-vanishing measure, to regain coercivity of $B$ on $H$ . The variational problem is thus well-posed, and we obtain

\phi^{n+1}=u+G.

Regularity for $\phi^{n+1}$ follows from the classical elliptic regularity for $u$ and regularity for $\rho^{n+1}$ .

If $\partial\Omega_{D}=\emptyset$ , then $B[u,1]=0$ requires the compatibility condition for the source

\int_{\Omega}\left(f(x)+\sum_{i=1}^{m}q_{i}\rho^{n+1}_{i}\right)dx=0.

Due to conservation of mass, this can be ensured by

\int_{\Omega}\left(f(x)+\sum_{i=1}^{m}q_{i}\rho^{\rm in}_{i}\right)dx=0.

With such compatibility condition the solution of this variational formulation exists but is not unique. In such case one can replace $H$ by

H_{*}=\left\{v\in H^{1},\quad\int_{\Omega}vdx=0\right\},

then by the Poincaré-Wirtinger inequality, $B$ is actually $H_{*}$ - coercive. The new variational problem hence admits a unique solution and is well-posed.

Finally we prove positivity of $\rho^{n+1}$ if $\rho^{n}\geq 0$ . Since $w=\rho_{i}^{n+1}e^{\psi^{n}_{i}}\in C(\bar{\Omega})\cap C^{2}(\Omega)$ , we let $x^{*}=\operatornamewithlimits{argmin}_{x\in\bar{\Omega}}w(x)$ , and distinct three cases:

(i) If $x^{*}\in\partial\Omega_{D}$ , then

w(x)\geq w(x^{*})=\rho^{b}_{i}(x^{*},t_{n+1})e^{\psi^{b}_{i}(x^{*},t_{n})}\geq 0,\quad x\in\bar{\Omega}.

(ii) If $x^{*}\in\Omega$ , then we can show that

w(x)\geq w(x^{*})\geq\rho^{n}_{i}(x^{*})e^{\psi^{n}_{i}(x^{*})}\geq 0,\quad x\in\bar{\Omega}.

In fact, from (A.1a) it follows

	$\displaystyle\rho^{n}_{i}(x)=$	$\displaystyle e^{-\psi^{n}_{i}(x)}w(x)-\tau\nabla\cdot(D_{i}(x)e^{-\psi^{n}_{i}(x)}\nabla w(x))$
	$\displaystyle=$	$\displaystyle e^{-\psi^{n}_{i}(x)}w(x)-\tau\nabla(D_{i}(x)e^{-\psi^{n}_{i}(x)})\cdot\nabla w(x)-\tau D_{i}(x)e^{-\psi^{n}_{i}(x)}\Delta w(x).$

This when evaluated at $x^{*}$ , using $\nabla w(x^{*})=0$ and $\Delta w(x^{*})\geq 0$ , gives

\rho^{n}_{i}(x^{*})\leq e^{-\psi^{n}_{i}(x^{*})}w(x^{*}).

(iii) For $x^{*}\in\partial\Omega_{N}$ . If $w(x^{*})\geq 0$ , the proof is complete. We proceed with the case that

w(x^{*})<0,\quad x^{*}\in\partial\Omega_{N}.

This is possible by the Hopf strong minimum principle.

Define the differential operator

L\xi:=\tau D_{i}(x)e^{-\psi^{n}_{i}(x)}\Delta\xi+\tau\nabla(D_{i}(x)e^{-\psi^{n}_{i}(x)})\cdot\nabla\xi-e^{-\psi^{n}_{i}(x)}\xi.

We then have $Lw=-\rho_{i}^{n}(x)\leq 0$ , and $w(x)\geq w(x^{*})$ for all $x\in\Omega$ . These together with $w(x^{*})<0$ allow us to apply Theorem 8 in [35] to conclude $(\nabla w(x^{*}))\cdot\mathbf{n}<0$ . This is a contradiction.
Collecting all three cases, we have $w(x)\geq 0$ for all $x\in\bar{\Omega}$ . ∎

Appendix B Proof of Theorem 3.4.

Proof.

(i) For fixed $i$ we sum (3.7) over all cells to get

\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|(\rho^{n+1}_{i,\alpha}-\rho^{n}_{i,\alpha})=\tau\sum_{j=1}^{d}\sum_{\alpha\in\mathcal{A}}\frac{|K_{\alpha}|}{h_{j}}(C_{i,\alpha+e_{j}/2}-C_{i,\alpha-e_{j}/2})=0,

where we used summation by parts and $C_{i,\alpha+e_{j}/2}=0$ for $x_{\alpha+e_{j}/2}\in\partial\Omega$ .
(ii) Set

S^{n}_{\alpha}=f_{\alpha}+\sum_{i=1}^{m}q_{i}\rho_{i,\alpha}^{n}

and

\psi^{*}_{i,\alpha}=\log\rho^{n+1}_{i,\alpha}+\frac{1}{k_{B}T}q_{i}\phi^{n}_{\alpha}+\frac{1}{k_{B}T}\mu_{i,\alpha}.

Using (3.12) we find that

	$\displaystyle E^{n+1}_{h}-E^{n}_{h}=$	$\displaystyle\sum_{\alpha\in\mathcal{A}}\sum_{i=1}^{m}\|K_{\alpha}\|\left((\rho^{n+1}_{i,\alpha}-\rho^{n}_{i,\alpha})\psi^{*}_{i,\alpha}+\rho^{n}_{i,\alpha}\log\frac{\rho^{n+1}_{i,\alpha}}{\rho^{n}_{\alpha}}\right)$		(B.1)
		$\displaystyle+\frac{1}{k_{B}T}\sum_{\alpha\in\mathcal{A}}\|K_{\alpha}\|\left(\frac{1}{2}S^{n+1}_{\alpha}\phi^{n+1}_{\alpha}+\frac{1}{2}S^{n}_{\alpha}\phi^{n}_{\alpha}-S^{n+1}_{\alpha}\phi^{n}_{\alpha}\right).$		(B.1)

Using $\log X\leq X-1$ for $X>0$ and the mass conservation, we have

\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|\rho^{n}_{i,\alpha}\log\frac{\rho^{n+1}_{i,\alpha}}{\rho^{n}_{\alpha}}\leq\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|(\rho^{n+1}_{i,\alpha}-\rho^{n}_{\alpha})=0.

Also one can verify that

\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|S^{n+1}_{\alpha}\phi^{n}_{\alpha}=\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|S^{n}_{\alpha}\phi^{n+1}_{\alpha},

with which we obtain

\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|\left(\frac{1}{2}S^{n+1}_{\alpha}\phi^{n+1}_{\alpha}+\frac{1}{2}S^{n}_{\alpha}\phi^{n}_{\alpha}-S^{n+1}_{\alpha}\phi^{n}_{\alpha}\right)=\frac{1}{2}\sum_{\alpha\in\mathcal{A}}|K_{\alpha}|(S^{n+1}_{\alpha}-S^{n}_{\alpha})(\phi^{n+1}_{\alpha}-\phi^{n}_{\alpha}).

Insertion of these into (B.1) gives

E^{n+1}_{h}-E^{n}_{h}\leq-\tau I^{n}+\tau^{2}I\!I^{n},

(B.2)

where

	$\displaystyle I^{n}$	$\displaystyle=-\sum_{\alpha\in\mathcal{A}}\sum_{i=1}^{m}\|K_{\alpha}\|\left(\frac{\rho^{n+1}_{i,\alpha}-\rho^{n}_{i,\alpha}}{\tau}\right)\psi^{*}_{i,\alpha},$
	$\displaystyle I\!I^{n}$	$\displaystyle=\frac{1}{2k_{B}T}\sum_{\alpha\in\mathcal{A}}\|K_{\alpha}\|\left(\frac{S^{n+1}_{\alpha}-S^{n}_{\alpha}}{\tau}\right)\left(\frac{\phi^{n+1}_{\alpha}-\phi^{n}_{\alpha}}{\tau}\right).$

By using (3.7) and summation by parts, we have

	$\displaystyle I^{n}=$	$\displaystyle-\sum_{i=1}^{m}\sum_{\alpha\in\mathcal{A}}\sum_{j=1}^{d}\|K_{\alpha}\|\left(\frac{C_{i,\alpha+e_{j}/2}-C_{i,\alpha-e_{j}/2}}{h_{j}}\right)\psi^{*}_{i,\alpha}$		(B.3)
	$\displaystyle=$	$\displaystyle\sum_{i=1}^{m}\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\frac{\|K_{\alpha}\|}{h_{j}}\left(\psi^{}_{i,\alpha+e_{j}}-\psi^{}_{i,\alpha}\right)C_{i,\alpha+e_{j}/2}.$		(B.3)

Note that

C_{i,\alpha+e_{j}/2}=\frac{D_{i}(x_{\alpha+e_{j}/2})e^{-\psi^{n}_{i,\alpha+e_{j}/2}}}{h_{j}}\left(e^{\psi^{*}_{i,\alpha+e_{j}}}-e^{\psi^{*}_{i,\alpha}}\right),

hence $I^{n}\geq 0$ .

We pause to discuss the special case with $I^{n}=0$ . In such case we must have $\psi^{*}_{i,\alpha+e_{j}}=\psi^{*}_{i,\alpha}$ for each $i,j$ and $\alpha\in\mathcal{A}$ , which implies $C_{i,\alpha+e_{j}/2}=0$ for each $i$ , $j$ and $\alpha\in\mathcal{A}$ . Thus, we have $\rho^{n+1}_{i,\alpha}=\rho^{n}_{i,\alpha}$ , hence

S^{n}_{\alpha}=f_{\alpha}+\sum_{i=1}^{m}q_{i}\rho_{i,\alpha}^{n}=f_{\alpha}+\sum_{i=1}^{m}q_{i}\rho_{i,\alpha}^{n+1}=S^{n+1}_{\alpha},\quad\forall\alpha\in\mathcal{A},

therefore $I\!I^{n}=0$ and $E^{n+1}_{h}-E^{n}_{h}\leq 0$ , this is (3.13) with $I^{n}=0$ .

From now on we only consider the case $I^{n}>0$ . We proceed to estimate $I\!I^{n}$ ,

$\displaystyle I\!I^{n}=$	$\displaystyle\frac{1}{2k_{B}T}\sum_{\alpha\in\mathcal{A}}\|K_{\alpha}\|\left(\frac{S^{n+1}_{\alpha}-S^{n}_{\alpha}}{\tau}\right)\left(\frac{\phi^{n+1}_{\alpha}-\phi^{n}_{\alpha}}{\tau}\right)$	(B.4)
$\displaystyle=$	$\displaystyle-\frac{1}{8\pi k_{B}T}\sum_{\alpha\in\mathcal{A}}\sum_{j=1}^{d}\frac{\|K_{\alpha}\|}{\tau^{2}h_{j}}(\Phi_{i,\alpha+e_{j}/2}^{n+1}-\Phi_{i,\alpha-e_{j}/2}^{n+1}-\Phi_{i,\alpha+e_{j}/2}^{n}+\Phi_{i,\alpha-e_{j}/2}^{n})(\phi^{n+1}_{\alpha}-\phi^{n}_{\alpha})$
$\displaystyle=$	$\displaystyle\frac{1}{8\pi k_{B}T}\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\frac{\|K_{\alpha}\|}{\tau^{2}h_{j}}(\Phi_{i,\alpha+e_{j}/2}^{n+1}-\Phi_{i,\alpha+e_{j}/2}^{n})(\phi^{n+1}_{\alpha+e_{j}}-\phi^{n}_{\alpha+e_{j}}-\phi^{n+1}_{\alpha}+\phi^{n}_{\alpha})$
$\displaystyle=$	$\displaystyle\frac{1}{8\pi k_{B}T}\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\|K_{\alpha}\|\frac{\epsilon_{\alpha+e_{j}/2}}{\tau^{2}h^{2}_{j}}(\phi^{n+1}_{\alpha+e_{j}}-\phi^{n}_{\alpha+e_{j}}-\phi^{n+1}_{\alpha}+\phi^{n}_{\alpha})^{2}\geq 0.$

Here the second equality is obtained by using the equation (3.11), the last equality is obtained by using the definition (3.11b) of $\Phi^{n}_{i,\alpha+e_{j}/2}$ .

From (B.3) and (B.4), we see that the energy dissipation inequality (3.13) is satisfied if

\tau\leq\tau^{*}\leq\frac{I^{n}}{2I\!I^{n}}.

(B.5)

In the remaining of the proof we will quantify $\tau^{*}$ from estimating the lower bound of $\frac{I^{n}}{2I\!I^{n}}.$

Subtracting (3.11) at time level $t=t_{n+1}$ and $t=t_{n}$ , one has

-\sum_{j=1}^{d}\frac{\Phi^{n+1}_{\alpha+e_{j}/2}-\Phi^{n+1}_{\alpha-e_{j}/2}-\Phi^{n}_{\alpha+e_{j}/2}+\Phi^{n}_{\alpha-e_{j}/2}}{h_{j}}=4\pi\sum_{i=1}^{m}q_{i}(\rho^{n+1}_{i,\alpha}-\rho^{n}_{i,\alpha}),

(B.6)

multiplying by $|K_{\alpha}|(\phi^{n+1}_{\alpha}-\phi^{n}_{\alpha})$ and summing over $\alpha\in\mathcal{A}$ leads to

		$\displaystyle-\sum_{j=1}^{d}\sum_{\alpha\in\mathcal{A}}\frac{K_{\alpha}}{h_{j}}(\Phi^{n+1}_{\alpha+e_{j}/2}-\Phi^{n+1}_{\alpha-e_{j}/2}-\Phi^{n}_{\alpha+e_{j}/2}+\Phi^{n}_{\alpha-e_{j}/2})(\phi^{n+1}_{\alpha}-\phi^{n}_{\alpha})$		(B.7)
		$\displaystyle=4\pi\sum_{i=1}^{m}\sum_{\alpha\in\mathcal{A}}q_{i}\|K_{\alpha}\|(\rho^{n+1}_{i,\alpha}-\rho^{n}_{i,\alpha})(\phi^{n+1}_{\alpha}-\phi^{n}_{\alpha}).$		(B.7)

Similar to (B.4), the left hand side of (B.7) reduces to

LHS=\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}|K_{\alpha}|\frac{\epsilon_{\alpha+e_{j}/2}}{h^{2}_{j}}(\phi^{n+1}_{\alpha+e_{j}}-\phi^{n}_{\alpha+e_{j}}-\phi^{n+1}_{\alpha}+\phi^{n}_{\alpha})^{2}.

(B.8)

We estimate the right hand side of (B.7) by using the equation (3.7):

$\displaystyle RHS=$	$\displaystyle 4\pi\sum_{i=1}^{m}\sum_{\alpha\in\mathcal{A}}q_{i}\|K_{\alpha}\|(\rho^{n+1}_{i,\alpha}-\rho^{n}_{i,\alpha})(\phi^{n+1}_{\alpha}-\phi^{n}_{\alpha})$	(B.9)
$\displaystyle=$	$\displaystyle 4\pi\tau\sum_{i=1}^{m}\sum_{\alpha\in\mathcal{A}}\sum_{j=1}^{d}q_{i}\|K_{\alpha}\|\frac{1}{h_{j}}(C_{i,\alpha+e_{j}/2}-C_{i,\alpha-e_{j}/2})(\phi^{n+1}_{\alpha}-\phi^{n}_{\alpha})$
$\displaystyle=$	$\displaystyle-4\pi\tau\sum_{i=1}^{m}\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}q_{i}\|K_{\alpha}\|\frac{1}{h_{j}}C_{i,\alpha+e_{j}/2}(\phi^{n+1}_{\alpha+e_{j}}-\phi^{n}_{\alpha+e_{j}}-\phi^{n+1}_{\alpha}+\phi^{n}_{\alpha}).$

Note that

LHS\geq\epsilon_{min}\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}|K_{\alpha}|\left(\frac{\phi^{n+1}_{\alpha+e_{j}}-\phi^{n}_{\alpha+e_{j}}-\phi^{n+1}_{\alpha}+\phi^{n}_{\alpha}}{h_{j}}\right)^{2}.

Using the Cauchy-Schwarz inequality we see that

$\displaystyle RHS$	$\displaystyle\leq 4\pi\tau\sum_{i=1}^{m}\|q_{i}\|\left(\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\|K_{\alpha}\|C_{i,\alpha+e_{j}/2}\left(\frac{\phi^{n+1}_{\alpha+e_{j}}-\phi^{n}_{\alpha+e_{j}}-\phi^{n+1}_{\alpha}+\phi^{n}_{\alpha}}{h_{j}}\right)\right)$	(B.10)
	$\displaystyle\leq 4\pi\tau\sum_{i=1}^{m}\|q_{i}\|\left[\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\|K_{\alpha}\|C^{2}_{i,\alpha+e_{j}/2}\right]^{1/2}$
	$\displaystyle\qquad\times\left[\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\|K_{\alpha}\|\left(\frac{\phi^{n+1}_{\alpha+e_{j}}-\phi^{n}_{\alpha+e_{j}}-\phi^{n+1}_{\alpha}+\phi^{n}_{\alpha}}{h_{j}}\right)^{2}\right]^{1/2}.$

We thus obtain

		$\displaystyle\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\|K_{\alpha}\|\left(\frac{\phi^{n+1}_{\alpha+e_{j}}-\phi^{n}_{\alpha+e_{j}}-\phi^{n+1}_{\alpha}+\phi^{n}_{\alpha}}{h_{j}}\right)^{2}$		(B.11)
		$\displaystyle\frac{\leq 16\pi^{2}\tau^{2}}{\epsilon^{2}_{min}}\left[\sum_{i=1}^{m}\|q_{i}\|\left(\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\|K_{\alpha}\|C^{2}_{i,\alpha+e_{j}/2}\right)^{1/2}\right]^{2}$
		$\displaystyle\leq\frac{16\pi^{2}\tau^{2}}{\epsilon^{2}_{min}}\left(\sum_{i=1}^{m}q^{2}_{i}\right)\sum_{i=1}^{m}\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\|K_{\alpha}\|C^{2}_{i,\alpha+e_{j}/2}.$

Upon insertion into (B.4)

I\!I^{n}\leq C\sum_{i=1}^{m}\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}|K_{\alpha}|C^{2}_{i,\alpha+e_{j}/2},

(B.12)

where $C=\frac{2\epsilon_{max}\pi\sum_{i=1}^{m}q_{i}^{2}}{\epsilon^{2}_{min}k_{B}T}$ . We use (B.3) and (B.12) to obtain:

$\displaystyle\frac{I^{n}}{2I\!I^{n}}\geq$	$\displaystyle\frac{\sum_{i=1}^{m}\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\frac{\|K_{\alpha}\|}{h_{j}}C_{i,\alpha+e_{j}/2}(\psi^{}_{i,\alpha+e_{j}}-\psi^{}_{i,\alpha})}{2C\sum_{i=1}^{m}\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\|K_{\alpha}\|C^{2}_{i,\alpha+e_{j}/2}}$	(B.13)
$\displaystyle\geq$	$\displaystyle\frac{1}{2C}\min_{i,j,\alpha}\left\{\frac{\psi^{}_{i,\alpha+e_{j}}-\psi^{}_{i,\alpha}}{h_{j}C_{i,\alpha+e_{j}/2}}\right\}$
$\displaystyle=$	$\displaystyle\frac{1}{2C}\min_{i,j,\alpha}\left\{\frac{\psi^{}_{i,\alpha+e_{j}}-\psi^{}_{i,\alpha}}{D_{i,\alpha+e_{j}/2}e^{-\psi^{n}_{i,\alpha+e_{j}/2}}(e^{\psi^{}_{i,\alpha+e_{j}}}-e^{\psi^{}_{i,\alpha}})}\right\}\quad\text{ by the mean-value theorem}$
$\displaystyle=$	$\displaystyle\frac{1}{2C}\min_{i,j,\alpha}\left\{\frac{1}{D_{i,\alpha+e_{j}/2}e^{-\psi^{n}_{i,\alpha+e_{j}/2}}e^{(\theta\psi^{}_{i,\alpha+e_{j}}+(1-\theta)\psi^{}_{i,\alpha})}}\right\},$

where $\theta\in(0,1).$ By using the harmonic mean for $e^{-\psi^{n}_{i,\alpha+e_{j}/2}}$ , we have

	$\displaystyle\frac{1}{e^{-\psi^{n}_{i,\alpha+e_{j}/2}}e^{(\theta\psi^{}_{i,\alpha+e_{j}}+(1-\theta)\psi^{}_{i,\alpha})}}=$	$\displaystyle\frac{e^{((\theta-1)\psi^{n}_{i,\alpha}-\theta\psi^{n}_{i,\alpha+e_{j}})}}{(\rho^{n+1}_{i,\alpha+e_{j}})^{\theta}(\rho^{n+1}_{i,\alpha})^{1-\theta}}\cdot\frac{2e^{\psi^{n}_{i,\alpha+e_{j}}+\psi^{n}_{i,\alpha}}}{e^{\psi^{n}_{i,\alpha+e_{j}}}+e^{\psi^{n}_{i,\alpha}}}$
	$\displaystyle=$	$\displaystyle\frac{1}{(\rho^{n+1}_{i,\alpha+e_{j}})^{\theta}(\rho^{n+1}_{i,\alpha})^{1-\theta}}\cdot\frac{2e^{(1-\theta)\psi^{n}_{i,\alpha+e_{j}}+\theta\psi^{n}_{i,\alpha}}}{e^{\psi^{n}_{i,\alpha+e_{j}}}+e^{\psi^{n}_{i,\alpha}}}$
	$\displaystyle\geq$	$\displaystyle\frac{2e^{\min\left\{\psi^{n}_{i,\alpha+e_{j}},\psi^{n}_{i,\alpha}\right\}}}{2Me^{\max\left\{\psi^{n}_{i,\alpha+e_{j}},{\psi^{n}_{i,\alpha}}\right\}}}$
	$\displaystyle=$	$\displaystyle\frac{e^{-\|\psi^{n}_{i,\alpha+e_{j}}-\psi^{n}_{i,\alpha}\|}}{M},$

where $M=\max_{i,\alpha,n}\rho^{n}_{i,\alpha}$ , thus

\frac{I^{n}}{2I\!I^{n}}\geq\frac{1}{2CD_{max}M}e^{-\max_{i,j,\alpha}|\psi^{n}_{i,\alpha+e_{j}}-\psi^{n}_{i,\alpha}|}.

(B.14)

For geometric mean or algebraic mean when used for the evaluation of $e^{-\psi^{n}_{i,\alpha+e_{j}/2}}$ we can verify either the same or bigger bound than the right hand side of in (B.14). ∎

References

[1] M.Z. Bazant, K. Thornton, and A. Ajdari. Diffuse-charge dynamics in electrochemical systems. Phys. Rev. E., 70: 021506, 2004.
[2] C. Buet and S. Dellacheris. On the Chang and Cooper scheme applied to a linear Fokker–Planck equation. Commun. Math. Sci., 8: 1079–1090, 2010.
[3] M. Burger, B. Schlake, and M.-T. Wolfram. Nonlinear Poisson-Nernst-Planck equations for ion flux through confined geometries. Nonlinearity., 25: 961–990, 2012.
[4] J. Chen, L.C. McInnes, and H. Zhang. Analysis and practical use of flexible BiCGStab. J. Sci. Comput., 68(2): 803–825, 2016.
[5] 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.
[6] D. Chen and R. Eisenberg. Poisson-Nernst-Planck (PNP) theory of open ionic channels. Biophys. J., 64: A22, 1993.
[7] B. Eisenberg. Ionic channels in biological membranes - electrostatic analysis of a natural nanotube. Contemp. Phys., 39 (6): 447–466, 1998.
[8] 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.
[9] A. Flavell, J. Kabre, and X. Li. An energy-preserving discretization for the Poisson-Nernst-Planck equations. J. Comput. Electron., 16: 431–441, 2017.
[10] H. Gao and D. He. Linearized conservative finite element methods for the Nernst-Planck-Poisson equations. J. Sci. Comput., 72: 1269–1289, 2017.
[11] C. Gardner, W. Nonner, and R. S. Eisenberg. Electrodiffusion model simulation of ionic channels: 1d simulation. J. Comput. Electron., 3: 25–31, 2004.
[12] D. Gillespie, W. Nonner, and R. S. Eisenberg. Density functional theory of charged hard-sphere fluids. Phys. Rev. E., 68: 0313503, 2003.
[13] D. Gilbarg and N. S. Trudinger. Elliptic Partial Differential Equations of Second Order. Grundlagen der mathematischen Wissenschaften, Vol. 224, Springer, Berlin, 1997.
[14] D. He and K. Pan. An energy preserving finite difference scheme for the Poisson-Nernst-Planck system. Appl. Math. Comput., 287: 214–223, 2016.
[15] 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.
[16] J.W. Hu and X.D. Huang. A fully discrete positivity-preserving and energy-dissipative finite difference scheme for Poisson-Nernst-Planck equations Preprint (2019).
[17] Y. Hyon, B. Eisenberg, and C. Liu. A mathematical model for the hard sphere repulsion in ionic solutions. Commun. Math. Sci., 9 (2): 459–475, 2011.
[18] 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.
[19] 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.
[20] B. Li. Continuum electrostatics for ionic solutions with nonuniform ionic sizes. Nonlinearity., 22: 811–833, 2009.
[21] W. Liu. Geometric singular perturbation approach to steady state Poisson-Nernst-Planck systems. SIAM J. Appl. Math., 65 (3): 754–766, 2005.
[22] H. Liu and W. Maimaitiyiming. Energy stable and unconditional positive schemes for a reduced Poisson-Nernst- Plank system. Comm. Comput. Phys (in press)., arXiv:1909.13161., 2019.
[23] H. Liu and W. Maimaitiyiming. A second order positive scheme for the reduced Poisson-Nernst-Planck system. J. Comp. Appl. Math., (2019), submitted
[24] H. Liu and W. Maimaitiyiming. Positive and free energy satisfying schemes for diffusion with interaction potentials. arXiv:1910.00151., 2019.
[25] H. Liu and H. Yu. An entropy satisfying conservative method for the Fokker-Planck equation of the finitely extensible nonlinear elastic dumbbell model. SIAM J. Numer. Anal., 50 (3): 1207–1239, 2012.
[26] H. Liu and Z. Wang. A free energy satisfying finite difference method for Poisson-Nernst-Planck equations. J. Comput. Phys., 268: 363–376, 2014.
[27] 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.
[28] X. Liu, Y. Qiao, and B. Lu. Analysis of the mean field free energy functional of electrolyte solution with non-homogenous boundary conditions and the generalized PB/PNP equations with inhomogeneous dielectric permittivity. SIAM J. Appl. Math., 78 (2): 1131–1154, 2018.
[29] B. Lu, M.J. Holst, J.A. McCammon, and Y. Zhou. Poisson-Nernst-Planck equations for simulating biomolecular diffusion-reaction processes I: finite element solutions. Journal of Computational Physics., 229: 6979–6994, 2010.
[30] P.A. Markowich, C.A. Ringhofer, and C. Schmeiser. Semiconductor Equations. Springer- Verlag, 1990.
[31] D. Meng, B. Zheng, G. Lin, and M. L. Sushko. Numerical Solution of 3D Poisson-Nernst-Planck Equations Coupled with Classical Density Functional Theory for Modeling Ion and Electron Transport in a Confined Environment. Communications in Computational Physics., 16: 1298–1322, 2014.
[32] M.S. Metti, J. Xu, and C. Liu. Energetically stable discretizations for charge transport and electrokinetic models. J. Comput. Phys., 306: 1–18, 2016 .
[33] J. Newman. Electrochemical Systems. Prentice Hall, 1991.
[34] 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.
[35] M.H. Protter and H. F. Weinberger. Maximum Principles in Differential Equations. Prentice-Hall, New York., 1967.
[36] Y. Qian, Z. Wang, and S. Zhou. A conservative, free energy dissipating, and positivity preserving finite difference scheme for multi-dimensional nonlocal Fokker-Planck equation. Journal of Computational Physics., 386: 22–36, 2019.
[37] Y. Saad. A flexible inner-outer preconditioned GMRES algorithm. SIAM J. Sci. Comput., 14(2): 461–469, 1993.
[38] S. Selberherr. Analysis and Simulation of Semiconductor Devices. Springer-Verlag/Wien,New York, 1984.
[39] F. Siddiqua, Z. Wang, and S. Zhou. A modified Poisson-Nernst-Planck model with excluded volume effect: theory and numerical implementation. Commun. Math. Sci., 16 (1): 251–271, 2018.
[40] T. Sokalski and A. Lewenstam. Application of Nernst-Planck and Poisson equations for interpretation of liquid-junction and membrane potentials in realtime and space domains. Electrochem. Commun., 3(3): 107–112, 2001.
[41] T. Sokalski, P. Lingenfelter and A. Lewenstam. Numerical solution of the coupled Nernst-Planck and Poisson equations for liquid junction and ion selective membrane potentials. J. Phys. Chem. B., 107(11): 2443–2452, 2003.
[42] M. L. Sushko, K. M. Rosso, J.-G. J. Zhang, and J. Liu. Multiscale simulations of Li ion conductivity in solid electrolyte. Chem. Phys. Lett., 2: 2352–2356, 2011.
[43] T. Teorell. Transport processes and electrical phenomena in ionic membranes. Progress Biophysics, 3: 305, 1953.
[44] G.W. Wei, Q. Zheng, Z. Chen, and K. Xia. Variational multiscale models for charge transport. SIAM Rev., 54: 699–754, 2012.
[45] P. Yin, Y. Huang, and H. Liu. An iterative discontinuous Galerkin method for solving the nonlinear Poisson-Boltzmann equation. Commun. Comput. Phys., 16: 491–515, 2014.
[46] P. Yin, Y. Huang, and H. Liu. Error estimates for the iterative discontinuous Galerkin method to the nonlinear Poisson-Boltzmann equation. Commun. Comput. Phys., 23: 168–197, 2018.
[47] Q. Zheng, D. Chen, and G.-W. Wei. Second-order Poisson-Nernst-Planck solver for ion transport. Journal of Computational Physics., 230: 5239–5262, 2011.
[48] Q. Zheng and G.W. Wei. Poisson-Boltzmann-Nernst-Planck model. Journal of Chemical Physics., 134: 194101, 2011.
[49] J. Zhou, Y. Jiang, and M. Doi. Cross interaction drives stratification in drying film of binary colloidal mixtures. Phys. Rev. Lett., 10: 108002, 2017.

$\displaystyle RHS$	$\displaystyle\leq 4\pi\tau\sum_{i=1}^{m}\|q_{i}\|\left(\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\|K_{\alpha}\|C_{i,\alpha+e_{j}/2}\left(\frac{\phi^{n+1}_{\alpha+e_{j}}-\phi^{n}_{\alpha+e_{j}}-\phi^{n+1}_{\alpha}+\phi^{n}_{\alpha}}{h_{j}}\right)\right)$	(B.10)
	$\displaystyle\leq 4\pi\tau\sum_{i=1}^{m}\|q_{i}\|\left[\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\|K_{\alpha}\|C^{2}_{i,\alpha+e_{j}/2}\right]^{1/2}$
	$\displaystyle\qquad\times\left[\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\|K_{\alpha}\|\left(\frac{\phi^{n+1}_{\alpha+e_{j}}-\phi^{n}_{\alpha+e_{j}}-\phi^{n+1}_{\alpha}+\phi^{n}_{\alpha}}{h_{j}}\right)^{2}\right]^{1/2}.$

		$\displaystyle\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\|K_{\alpha}\|\left(\frac{\phi^{n+1}_{\alpha+e_{j}}-\phi^{n}_{\alpha+e_{j}}-\phi^{n+1}_{\alpha}+\phi^{n}_{\alpha}}{h_{j}}\right)^{2}$		(B.11)
		$\displaystyle\frac{\leq 16\pi^{2}\tau^{2}}{\epsilon^{2}_{min}}\left[\sum_{i=1}^{m}\|q_{i}\|\left(\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\|K_{\alpha}\|C^{2}_{i,\alpha+e_{j}/2}\right)^{1/2}\right]^{2}$
		$\displaystyle\leq\frac{16\pi^{2}\tau^{2}}{\epsilon^{2}_{min}}\left(\sum_{i=1}^{m}q^{2}_{i}\right)\sum_{i=1}^{m}\sum_{j=1}^{d}\sum_{\alpha(j)\neq N_{j}}\|K_{\alpha}\|C^{2}_{i,\alpha+e_{j}/2}.$