A Second Order Accurate, Positivity Preserving Numerical Method for the Poisson–Nernst–Planck System and Its Convergence Analysis

Abstract

A second order accurate (in time) numerical scheme is proposed and analyzed for the Poisson–Nernst–Planck equation (PNP) system, reformulated as a non-constant mobility $H^{-1}$ gradient flow in the Energetic Variational Approach (EnVarA). The centered finite difference is taken as the spatial discretization. Meanwhile, the highly nonlinear and singular nature of the logarithmic energy potentials has always been the essential difficulty to design a second order accurate scheme in time, while preserving the variational energetic structures. The mobility function is updated with a second order accurate extrapolation formula, for the sake of unique solvability. A modified Crank–Nicolson scheme is used to approximate the logarithmic term, so that its inner product with the discrete temporal derivative exactly gives the corresponding nonlinear energy difference; henceforth the energy stability is ensured for the logarithmic part. In addition, nonlinear artificial regularization terms are added in the numerical scheme, so that the positivity-preserving property could be theoretically proved, with the help of the singularity associated with the logarithmic function. Furthermore, an optimal rate convergence analysis is provided in this paper, in which the higher order asymptotic expansion for the numerical solution, the rough error estimate and refined error estimate techniques have to be included to accomplish such an analysis. This work combines the following theoretical properties for a second order accurate numerical scheme for the PNP system: (i) second order accuracy in both time and space, (ii) unique solvability and positivity, (iii) energy stability, and (iv) optimal rate convergence. A few numerical results are also presented.

Appendices

Appendix A: Proof of Proposition 4.1

In terms of the temporal discretization, the following local truncation error can be derived by a Taylor expansion in time, combined with the projection estimate (4.3):

$$\begin{aligned} \frac{{{\textsf{N}}}_N^{m+1} - {{\textsf{N}}}_N^m}{{\Delta t}}= & {} \nabla \cdot \big ( \breve{{\textsf{N}}}_N^{m + \nicefrac {1}{2}} \nabla \big ( G_{{{\textsf{N}}}_N^m}^1 \big ( {{\textsf{N}}}_N^{m+1} \big ) + (-\Delta )^{-1} \big ( {{\textsf{N}}}_N^{m+\nicefrac {1}{2}} - {{\textsf{P}}}_N^{m+\nicefrac {1}{2}} \big ) \nonumber \\{} & {} + {\Delta t}\big ( \ln {{\textsf{N}}}_N^{m+1} - \ln {{\textsf{N}}}_N^m \big ) \big ) \big ) + {\Delta t}^2 \big ( G_n^{(0)} \big )^{m+\nicefrac {1}{2}} + O\big ({\Delta t}^3\big ) + O \big (h^{m_0}\big ), \nonumber \\ \end{aligned}$$

(A.1)

$$\begin{aligned} \frac{ {{\textsf{P}}}_N^{m+1} - {{\textsf{P}}}_N^m}{{\Delta t}}= & {} \nabla \cdot \big ( D \breve{\mathsf P}_N^{m + \nicefrac {1}{2}} \nabla \big ( G_{{{\textsf{P}}}_N^m}^1 \big ( {{\textsf{N}}}_P^{m+1} \big ) + (-\Delta )^{-1} \big ( {{\textsf{P}}}_N^{m+\nicefrac {1}{2}} - {{\textsf{N}}}_N^{m+\nicefrac {1}{2}} \big )\nonumber \\{} & {} + {\Delta t}\big ( \ln {{\textsf{P}}}_N^{m+1} - \ln {{\textsf{P}}}_N^m \big ) \big ) \big ) + {\Delta t}^2 \big ( G_p^{(0)} \big )^{m+\nicefrac {1}{2}} + O\big ({\Delta t}^3\big ) + O \big (h^{m_0}\big ), \nonumber \\ \breve{{\textsf{N}}}_N^{m + \nicefrac {1}{2}}= & {} \frac{3}{2} {{\textsf{N}}}_N^m - \frac{1}{2} {{\textsf{N}}}_N^{m-1}, \quad \breve{{\textsf{P}}}_N^{m + \nicefrac {1}{2}} = \frac{3}{2} {{\textsf{P}}}_N^m - \frac{1}{2} {{\textsf{P}}}_N^{m-1}, \nonumber \\ {{\textsf{N}}}_N^{m+\nicefrac {1}{2}}= & {} \frac{1}{2} \big ( {{\textsf{N}}}_N^{m+1} + {\mathsf N}_N^m \big ), \quad {{\textsf{P}}}_N^{m+\nicefrac {1}{2}} = \frac{1}{2} \big ( {\mathsf P}_N^{m+1} + {{\textsf{P}}}_N^m \big ), \end{aligned}$$

(A.2)

with the projection accuracy order $m_0 \ge 4$. In fact, the spatial functions $G_n^{(0)}$ $G_p^{(0)}$ are smooth enough in the sense that their derivatives are bounded.

Subsequently, the leading order temporal correction function $({{\textsf{N}}}_{{\Delta t}, 1}, {{\textsf{P}}}_{{\Delta t}, 1})$ turns out to be the solution of the following linear equations:

$$\begin{aligned} \partial _t {{\textsf{N}}}_{{\Delta t}, 1}= & {} \nabla \cdot \Big ( {{\textsf{N}}}_{{\Delta t}, 1} \nabla \big ( \ln {{\textsf{N}}}_N + (-\Delta )^{-1} \big ( {{\textsf{N}}}_N - {{\textsf{P}}}_N \big ) \big ) \nonumber \\{} & {} \quad + {{\textsf{N}}}_N \nabla \big ( \frac{1}{{{\textsf{N}}}_N} {{\textsf{N}}}_{{\Delta t}, 1} + (-\Delta )^{-1} \big ( {{\textsf{N}}}_{{\Delta t}, 1} - {{\textsf{P}}}_{{\Delta t}, 1} \big ) \big ) \Big ) - G_n^{(0)}, \end{aligned}$$

(A.3)

$$\begin{aligned} \partial _t {{\textsf{P}}}_{{\Delta t}, 1}= & {} \nabla \cdot \Big ( D {{\textsf{P}}}_{{\Delta t}, 1} \nabla \big ( \ln {{\textsf{P}}}_N + (-\Delta )^{-1} \big ( {{\textsf{P}}}_N - {{\textsf{N}}}_N \big ) \big ) \nonumber \\{} & {} \quad + D {{\textsf{P}}}_N \nabla \big ( \frac{1}{{{\textsf{P}}}_N} {{\textsf{P}}}_{{\Delta t}, 1} + (-\Delta )^{-1} \big ( {{\textsf{P}}}_{{\Delta t}, 1} - {{\textsf{N}}}_{{\Delta t}, 1} \big ) \big ) \Big ) - G_p^{(0)}. \end{aligned}$$

(A.4)

In fact, existence of a solution of the above linear and parabolic PDE system is straightforward. It depends only on the projection solution $({{\textsf{N}}}_N, {{\textsf{P}}}_N)$. And also, the derivatives of $({{\textsf{N}}}_{{\Delta t}, 1}, {{\textsf{P}}}_{{\Delta t}, 1})$ in various orders are bounded. In turn, an application of the semi-implicit discretization to (A.3)–(A.4) gives

$$\begin{aligned} \frac{ {{\textsf{N}}}_{{\Delta t}, 1}^{m+1} - {{\textsf{N}}}_{{\Delta t}, 1}^m}{{\Delta t}}= & {} \nabla \cdot \Big ( \breve{{\textsf{N}}}_{{\Delta t}, 1}^{m+\nicefrac {1}{2}} \nabla \big ( G_{{{\textsf{N}}}_N^m}^1 \big ( {{\textsf{N}}}_N^{m+1} \big ) + \big (-\Delta )^{-1} \big ( {{\textsf{N}}}_N^{m+\nicefrac {1}{2}} - {{\textsf{P}}}_N^{m+\nicefrac {1}{2}} \big ) \big ) \nonumber \\{} & {} \, \, + \breve{{\textsf{N}}}_N^{m+\nicefrac {1}{2}} \nabla ( \frac{1}{{\mathsf N}_N^{m+\nicefrac {1}{2}} } {{\textsf{N}}}_{{\Delta t}, 1}^{m+\nicefrac {1}{2}} + (-\Delta )^{-1} \big ( {{\textsf{N}}}_{{\Delta t}, 1}^{m+\nicefrac {1}{2}} - {{\textsf{P}}}_{{\Delta t}, 1}^{m+\nicefrac {1}{2}} \big ) \big ) \Big ) \nonumber \\{} & {} \quad - \big ( G_n^{(0)} \big )^{m+\nicefrac {1}{2}} + {\Delta t}^2 {\varvec{h}}_1^{m+\nicefrac {1}{2}} + O \big ({\Delta t}^3\big ), \end{aligned}$$

(A.5)

$$\begin{aligned} \frac{ {{\textsf{P}}}_{{\Delta t}, 1}^{m+1} - {{\textsf{P}}}_{{\Delta t}, 1}^m}{{\Delta t}}= & {} \nabla \cdot \Big ( D \breve{{\textsf{P}}}_{{\Delta t}, 1}^{m+\nicefrac {1}{2}} \nabla \big ( G_{{{\textsf{P}}}_N^m}^1 \big ( {{\textsf{P}}}_N^{m+1} \big ) + \big (-\Delta )^{-1} \big ( {{\textsf{P}}}_N^{m+\nicefrac {1}{2}} - {{\textsf{N}}}_N^{m+\nicefrac {1}{2}} \big ) \big ) \nonumber \\{} & {} \, \, + D \breve{{\textsf{P}}}_N^{m+\nicefrac {1}{2}} \nabla ( \frac{1}{{{\textsf{P}}}_N^{m+\nicefrac {1}{2}} } {{\textsf{P}}}_{{\Delta t}, 1}^{m+\nicefrac {1}{2}} + \big (-\Delta )^{-1} \big ( {{\textsf{P}}}_{{\Delta t}, 1}^{m+\nicefrac {1}{2}} - {{\textsf{N}}}_{{\Delta t}, 1}^{m+\nicefrac {1}{2}} \big ) \big ) \Big ) \nonumber \\{} & {} \quad - \big ( G_p^{(0)} \big )^{m+\nicefrac {1}{2}} + {\Delta t}^2 {\varvec{h}}_2^{m+\nicefrac {1}{2}} + O ({\Delta t}^3), \nonumber \\ \breve{{\textsf{N}}}_{{\Delta t}, 1}^{m+\nicefrac {1}{2}}= & {} \frac{3}{2} {{\textsf{N}}}_{{\Delta t}, 1}^m - \frac{1}{2} {\mathsf N}_{{\Delta t}, 1}^{m-1}, \quad \breve{{\textsf{P}}}_{{\Delta t}, 1}^{m+\nicefrac {1}{2}} = \frac{3}{2} {{\textsf{P}}}_{{\Delta t}, 1}^m - \frac{1}{2} {{\textsf{P}}}_{{\Delta t}, 1}^{m-1}, \nonumber \\ {{\textsf{N}}}_{{\Delta t}, 1}^{m+\nicefrac {1}{2}}= & {} \frac{1}{2} \big ( {{\textsf{N}}}_{{\Delta t}, 1}^{m+1} + {{\textsf{N}}}_{{\Delta t}, 1}^m \big ), \quad {{\textsf{P}}}_{{\Delta t}, 1}^{m+\nicefrac {1}{2}} = \frac{1}{2} \big ( {{\textsf{P}}}_{{\Delta t}, 1}^{m+1} + {\mathsf P}_{{\Delta t}, 1}^m \big ). \end{aligned}$$

(A.6)

A combination of (A.1)–(A.2) and (A.5)–(A.6) results in the third order temporal truncation error for $\hat{{\textsf{N}}}_1:= {\mathsf N}_N + {\Delta t}^2 {{{\mathcal {P}}}}_N {{\textsf{N}}}_{{\Delta t}, 1}$, $\hat{{\textsf{P}}}_1:= {{\textsf{P}}}_N + {\Delta t}^2 {{{\mathcal {P}}}}_N {{\textsf{P}}}_{{\Delta t}, 1}$:

$$\begin{aligned} \frac{\hat{{\textsf{N}}}_1^{m+1} - \hat{{\textsf{N}}}_1^m}{{\Delta t}}= & {} \nabla \cdot \Big ( \big ( \frac{3}{2} \hat{{\textsf{N}}}_1^m - \frac{1}{2} \hat{\mathsf N}_1^{m-1} \big ) \nabla \big ( G_{\hat{{\textsf{N}}}_1^m}^1 \big ( \hat{\mathsf N}_1^{m+1} \big ) + (-\Delta )^{-1} \big ( \hat{{\textsf{N}}}_1^{m+\nicefrac {1}{2}} - \hat{{\textsf{P}}}_1^{m+\nicefrac {1}{2}} \big ) \nonumber \\{} & {} + {\Delta t}\big ( \ln \hat{{\textsf{N}}}_1^{m+1} - \ln \hat{{\textsf{N}}}_1^m \big ) \big ) \Big ) + {\Delta t}^3 \big ( G_n^{(1)} \big )^{m+\nicefrac {1}{2}} + O({\Delta t}^4) + O \big (h^{m_0}\big ), \nonumber \\ \end{aligned}$$

(A.7)

$$\begin{aligned} \frac{\hat{{\textsf{P}}}_1^{m+1} - \hat{{\textsf{P}}}_1^m}{{\Delta t}}= & {} \nabla \cdot \Big ( D \big ( \frac{3}{2} \hat{{\textsf{P}}}_1^m - \frac{1}{2} \hat{{\textsf{P}}}_1^{m-1} \big ) \nabla \big ( G_{\hat{{\textsf{P}}}_1^m}^1 ( \hat{{\textsf{P}}}_1^{m+1} \big ) + (-\Delta )^{-1} \big ( \hat{{\textsf{P}}}_1^{m+\nicefrac {1}{2}} - \hat{\mathsf N}_1^{m+\nicefrac {1}{2}} \big ) \nonumber \\{} & {} + {\Delta t}\big ( \ln \hat{{\textsf{P}}}_1^{m+1} - \ln \hat{{\textsf{P}}}_1^m \big ) \big ) \Big ) + {\Delta t}^3 \big ( G_p^{(1)} \big )^{m+\nicefrac {1}{2}} + O({\Delta t}^4) + O \big (h^{m_0}\big ), \nonumber \\ \hat{{\textsf{N}}}_1^{m+\nicefrac {1}{2}}= & {} \frac{1}{2} \big ( \hat{{\textsf{N}}}_1^{m+1} + \hat{{\textsf{N}}}_1^m \big ), \quad \check{{\textsf{P}}}_1^{m+\nicefrac {1}{2}} = \frac{1}{2} \big ( \hat{{\textsf{P}}}_1^{m+1} + \hat{{\textsf{P}}}_1^m \big ). \end{aligned}$$

(A.8)

In the derivation of (A.7)–(A.8), the following linearized expansions have been applied:

$$\begin{aligned} \begin{aligned} G_{\hat{{\textsf{N}}}_1^m}^1 \big ( \hat{{\textsf{N}}}_1^{m+1} \big ) =&\, G_{{{\textsf{N}}}_N^m}^1 \big ( {{\textsf{N}}}_N^{m+1} \big ) + \frac{1}{ 2 {\mathsf N}_N^m} \big ( \hat{{\textsf{N}}}_1^m - {{\textsf{N}}}_N^m \big ) + \frac{1}{ 2 {{\textsf{N}}}_N^{m+1} } \big ( \hat{{\textsf{N}}}_1^{m+1} - {\mathsf N}_N^{m+1} \big ) + O ({\Delta t}^3)\\ =&\,G_{{{\textsf{N}}}_N^m}^1 \big ( {{\textsf{N}}}_N^{m+1} \big ) + \frac{1}{2} {\Delta t}^2 \Big ( \frac{1}{ {{\textsf{N}}}_N^m} \cdot {{\textsf{N}}}_{{\Delta t}, 1}^m + \frac{1}{ {{\textsf{N}}}_N^{m+1} } \cdot {{\textsf{N}}}_{{\Delta t}, 1}^{m+1} \Big ) + O ({\Delta t}^3),\\ \frac{1}{{{\textsf{N}}}_N^{m+\nicefrac {1}{2}} } {{\textsf{N}}}_{{\Delta t}, 1}^{m+\nicefrac {1}{2}} =&\, \frac{1}{2} \Big ( \frac{1}{ {{\textsf{N}}}_N^m} \cdot {{\textsf{N}}}_{{\Delta t}, 1}^m + \frac{1}{ {{\textsf{N}}}_N^{m+1} } \cdot {{\textsf{N}}}_{{\Delta t}, 1}^{m+1} \Big ) + O \big ({\Delta t}^2 \big ), \end{aligned} \end{aligned}$$

(A.9)

in which property (3) of $G_a^1 (x)$ (as stated in Lemma 3.1) is recalled. The corresponding expansions for $G_{\hat{{\textsf{P}}}_1^m}^1 ( \hat{{\textsf{P}}}_1^{m+1} )$ could be similarly derived, and the technical details are skipped for the sake of brevity.

Similarly, the next order temporal correction function $\big ({\mathsf N}_{{\Delta t}, 2}, {{\textsf{P}}}_{{\Delta t}, 2}\big )$ turns out to be the solution of following linear equations:

$$\begin{aligned} \partial _t {{\textsf{N}}}_{{\Delta t}, 2}= & {} \nabla \cdot \big ( {{\textsf{N}}}_{{\Delta t}, 2} \nabla \big ( \ln \check{{\textsf{N}}}_1 + (-\Delta )^{-1} \big ( \check{{\textsf{N}}}_1 - \check{{\textsf{P}}}_1 \big ) \big ) \nonumber \\{} & {} \quad + \check{{\textsf{N}}}_1 \nabla \big ( \frac{1}{\check{\mathsf N}_1} {{\textsf{N}}}_{{\Delta t}, 2} + (-\Delta )^{-1} \big ( {{\textsf{N}}}_{{\Delta t}, 2} - {{\textsf{P}}}_{{\Delta t}, 2} \big ) \big ) \big ) - G_n^{(1)}, \end{aligned}$$

(A.10)

$$\begin{aligned} \partial _t {{\textsf{P}}}_{{\Delta t}, 2}= & {} \nabla \cdot \big ( D {{\textsf{P}}}_{{\Delta t}, 2} \nabla \big ( \ln \check{{\textsf{P}}}_1 + (-\Delta )^{-1} \big ( \check{{\textsf{P}}}_1 - \check{{\textsf{N}}}_1 \big ) \big ) \nonumber \\{} & {} \quad + D \check{{\textsf{P}}}_1 \nabla \big ( \frac{1}{\check{\mathsf P}_1} {{\textsf{P}}}_{{\Delta t}, 2} + (-\Delta )^{-1} \big ( {{\textsf{P}}}_{{\Delta t}, 2} - {{\textsf{N}}}_{{\Delta t}, 2} \big ) \big ) \big ) - G_p^{(1)}. \end{aligned}$$

(A.11)

Again, the solution depends only on the exact solution $({\mathsf N}, {{\textsf{P}}})$, with derivatives of various orders stay bounded. Of course, an application of the semi-implicit discretization to (A.10)–(A.11) gives

$$\begin{aligned} \frac{{{\textsf{N}}}_{{\Delta t}, 2}^{m+1} - {{\textsf{N}}}_{{\Delta t}, 2}^m}{{\Delta t}}= & {} \nabla \cdot \big ( \breve{{\textsf{N}}}_{{\Delta t}, 2}^{m+\nicefrac {1}{2}} \nabla \big ( G_{\hat{{\textsf{N}}}_1^m}^1 \big ( \hat{{\textsf{N}}}_1^{m+1} \big ) + (-\Delta )^{-1} \big ( \hat{{\textsf{N}}}_1^{m+\nicefrac {1}{2}} - \hat{\mathsf P}_1^{m+\nicefrac {1}{2}} \big ) \big ) \nonumber \\{} & {} + \bigg ( \frac{3}{2} \hat{{\textsf{N}}}_1^m - \frac{1}{2} \hat{\mathsf N}_1^{m-1} \bigg ) \nabla \big ( \frac{1}{\hat{{\textsf{N}}}_1^{m+\nicefrac {1}{2}} } {\mathsf N}_{{\Delta t}, 2}^{m+\nicefrac {1}{2}} + (-\Delta )^{-1} \big ( {{\textsf{N}}}_{{\Delta t}, 2}^{m+\nicefrac {1}{2}} - {{\textsf{P}}}_{{\Delta t}, 2}^{m+\nicefrac {1}{2}} \big ) \big ) \big ) \nonumber \\{} & {} - \big ( G_n^{(1)} \big )^{m+\nicefrac {1}{2}} + O \big ({\Delta t}^2\big ), \end{aligned}$$

(A.12)

$$\begin{aligned} \frac{{{\textsf{P}}}_{{\Delta t}, 2}^{m+1} - {{\textsf{P}}}_{{\Delta t}, 2}^m}{{\Delta t}}= & {} \nabla \cdot \Big ( D \breve{{\textsf{P}}}_{{\Delta t}, 2}^{m+\nicefrac {1}{2}} \nabla ( G_{\hat{{\textsf{P}}}_1^m}^1 \big ( \hat{\mathsf P}_1^{m+1} \big ) + (-\Delta )^{-1} \big ( \hat{{\textsf{P}}}_1^{m+\nicefrac {1}{2}} - \hat{{\textsf{N}}}_1^{m+\nicefrac {1}{2}} \big ) \big ) \nonumber \\{} & {} + D \big ( \frac{3}{2} \hat{{\textsf{P}}}_1^m - \frac{1}{2} \hat{\mathsf P}_1^{m-1} \big ) \nabla ( \frac{1}{\hat{{\textsf{P}}}_1^{m+\nicefrac {1}{2}} } {\mathsf P}_{{\Delta t}, 2}^{m+\nicefrac {1}{2}} + (-\Delta )^{-1} \big ( {{\textsf{P}}}_{{\Delta t}, 2}^{m+\nicefrac {1}{2}} - {{\textsf{N}}}_{{\Delta t}, 2}^{m+\nicefrac {1}{2}} \big ) \big ) \big ) \nonumber \\{} & {} - \big ( G_p^{(1)} \big )^{m+\nicefrac {1}{2}} + O \big ({\Delta t}^2\big ), \nonumber \\ \breve{{\textsf{N}}}_{{\Delta t}, 2}^{m+\nicefrac {1}{2}}= & {} \frac{3}{2} {\mathsf N}_{{\Delta t}, 2}^{m} - \frac{1}{2} {{\textsf{N}}}_{{\Delta t}, 2}^{m-1}, \quad \breve{{\textsf{P}}}_{{\Delta t}, 2}^{m+\nicefrac {1}{2}} = \frac{3}{2} {{\textsf{P}}}_{{\Delta t}, 2}^{m} - \frac{1}{2} {{\textsf{P}}}_{{\Delta t}, 2}^{m-1}. \end{aligned}$$

(A.13)

As a result, a combination of (A.10)–(A.11) and (A.12)–(A.13) yields the fourth order temporal truncation error for $\hat{{\textsf{N}}}_2:= \hat{{\textsf{N}}}_1 + {\Delta t}^3 {{{\mathcal {P}}}}_N {{\textsf{N}}}_{{\Delta t}, 2}$, $\hat{{\textsf{P}}}_2:= \hat{{\textsf{P}}}_1 + {\Delta t}^3 {{{\mathcal {P}}}}_N {\mathsf P}_{{\Delta t}, 2}$:

$$\begin{aligned} \frac{\hat{{\textsf{N}}}_2^{m+1} - \hat{{\textsf{N}}}_2^m}{{\Delta t}}= & {} \nabla \cdot \big ( \big ( \frac{3}{2} \hat{{\textsf{N}}}_2^m - \frac{1}{2} \hat{\mathsf N}_2^{m-1} \big ) \nabla \big ( G_{\hat{{\textsf{N}}}_2^m}^1 \big ( \hat{\mathsf N}_2^{m+1} ) + (-\Delta )^{-1} \big ( \hat{{\textsf{N}}}_2^{m+\nicefrac {1}{2}} - \hat{{\textsf{P}}}_2^{m+\nicefrac {1}{2}} \big ) \big ) \nonumber \\{} & {} + {\Delta t}\big ( \ln \hat{{\textsf{N}}}_2^{m+1} - \ln \hat{{\textsf{N}}}_2^m \big ) \big ) + {\Delta t}^4 \big ( G_n^{(2)} \big )^{m+\nicefrac {1}{2}} + O\big ({\Delta t}^5\big ) + O \big (h^{m_0}\big ), \nonumber \\ \end{aligned}$$

(A.14)

$$\begin{aligned} \frac{\hat{{\textsf{P}}}_2^{m+1} - \hat{{\textsf{P}}}_2^m}{{\Delta t}}= & {} \nabla \cdot \Big ( D \big ( \frac{3}{2} \hat{{\textsf{P}}}_2^m - \frac{1}{2} \hat{{\textsf{P}}}_2^{m-1} \big ) \nabla \big ( G_{\hat{{\textsf{P}}}_2^m}^1 \big ( \hat{{\textsf{P}}}_2^{m+1} \big ) + (-\Delta )^{-1} \big ( \hat{{\textsf{P}}}_2^{m+\nicefrac {1}{2}} - \hat{\mathsf N}_2^{m+\nicefrac {1}{2}} \big ) \big ) \nonumber \\{} & {} + {\Delta t}\big ( \ln \hat{{\textsf{P}}}_2^{m+1} - \ln \hat{{\textsf{P}}}_2^m \big ) \Big ) + {\Delta t}^4 \big ( G_p^{(2)} \big )^{m+\nicefrac {1}{2}} + O\big ({\Delta t}^5\big ) + O \big (h^{m_0}\big ), \nonumber \\ \end{aligned}$$

(A.15)

in which similar linearized expansions (as in (A.9)) have been used in the derivation.

In terms of spatial discretization, we construct the spatial correction term $({{\textsf{N}}}_{h, 1}, {{\textsf{P}}}_{h, 1})$ to improve the spatial accuracy order. The following truncation error estimate for the spatial discretization is available, by using a straightforward Taylor expansion for the constructed profile $(\hat{{\textsf{N}}}_2, \hat{{\textsf{P}}}_2)$:

$$\begin{aligned} \frac{\hat{{\textsf{N}}}_2^{m+1} - \hat{{\textsf{N}}}_2^m}{{\Delta t}}= & {} \nabla _h \cdot \big ( \big ( \frac{3}{2} \hat{{\textsf{N}}}_2^m - \frac{1}{2} \hat{{\textsf{N}}}_2^{m-1} \big ) \nabla _h \big ( G_{\hat{{\textsf{N}}}_2^m}^1 \big ( \hat{{\textsf{N}}}_2^{m+1} \big ) + (-\Delta _h)^{-1} \big ( \hat{\mathsf N}_2^{m+\nicefrac {1}{2}} - \hat{{\textsf{P}}}_2^{m+\nicefrac {1}{2}} \big ) \big ) \nonumber \\{} & {} + {\Delta t}\big ( \ln \hat{{\textsf{N}}}_2^{m+1} - \ln \hat{{\textsf{N}}}_2^m \big ) \Big ) + h^2 \big ( H_n^{(0)} \big )^{m+\nicefrac {1}{2}} + O \big ({\Delta t}^4 + h^4\big ), \end{aligned}$$

(A.16)

$$\begin{aligned} \frac{\hat{{\textsf{P}}}_2^{m+1} - \hat{{\textsf{P}}}_2^m}{{\Delta t}}= & {} \nabla _h \cdot \big ( D \big ( \frac{3}{2} \hat{{\textsf{P}}}_2^m - \frac{1}{2} \hat{{\textsf{P}}}_2^{m-1} \big ) \nabla _h \big ( G_{\hat{{\textsf{P}}}_2^m}^1 \big ( \hat{{\textsf{P}}}_2^{m+1} \big ) + (-\Delta _h)^{-1} \big ( \hat{{\textsf{P}}}_2^{m+\nicefrac {1}{2}} - \hat{\mathsf N}_2^{m+\nicefrac {1}{2}} \big ) \big ) \nonumber \\{} & {} + {\Delta t}\big ( \ln \hat{{\textsf{P}}}_2^{m+1} - \ln \hat{{\textsf{P}}}_2^m \big ) \big ) + h^2 \big ( H_p^{(0)} \big )^{m+\nicefrac {1}{2}} + O ({\Delta t}^4 + h^4), \end{aligned}$$

(A.17)

in which the average operator is taken in a similar form as (2.9). Similarly, the spatially discrete functions $H_n^{(0)}$, $H_p^{(0)}$ are smooth enough in the sense that their discrete derivatives are bounded. Because of the symmetry in the centered finite difference approximation, there is no $O (h^3)$ truncation error term. In turn, the spatial correction function $({{\textsf{N}}}_{h, 1}, {{\textsf{P}}}_{h, 1})$ is determined by solving the solution of the following linear PDE system:

$$\begin{aligned} \partial _t {{\textsf{N}}}_{h, 1}= & {} \nabla \cdot \Big ( {{\textsf{N}}}_{h, 1} \nabla \big ( \ln \check{{\textsf{N}}}_2 + (-\Delta )^{-1} \big ( \check{{\textsf{N}}}_2 - \check{{\textsf{P}}}_2 \big ) \big ) \nonumber \\{} & {} \quad + \check{{\textsf{N}}}_2 \nabla \big ( \frac{1}{\check{\mathsf N}_2} {{\textsf{N}}}_{h, 1} + (-\Delta )^{-1} \big ( {{\textsf{N}}}_{h, 1} - {{\textsf{P}}}_{h, 1} \big ) \big ) \big ) - H_n^{(0)}, \end{aligned}$$

(A.18)

$$\begin{aligned} \partial _t {{\textsf{P}}}_{h, 1}= & {} \nabla \cdot \big ( D {{\textsf{P}}}_{h, 1} \nabla \big ( \ln \check{\mathsf P}_2 + (-\Delta )^{-1} ( \check{{\textsf{P}}}_2 - \check{{\textsf{N}}}_2 \big ) ) \nonumber \\{} & {} \quad + D \check{{\textsf{P}}}_1 \nabla \big ( \frac{1}{\check{\mathsf P}_2} {{\textsf{P}}}_{h, 1} + (-\Delta )^{-1} \big ( {{\textsf{P}}}_{h, 1} - {{\textsf{N}}}_{h, 1} \big ) \big ) \big ) - H_p^{(0)}. \end{aligned}$$

(A.19)

Again, the solution depends only on the exact solution $({\mathsf N}, {{\textsf{P}}})$, with the divided differences of various orders stay bounded. An application of a full discretization to (A.18)–(A.19) leads to

$$\begin{aligned} \frac{{{\textsf{N}}}_{h, 1}^{m+1} - {{\textsf{N}}}_{h, 1}^m}{{\Delta t}}= & {} \nabla _h \cdot \Big ( \Big ( \breve{{\textsf{N}}}_{h, 1}^{m+\nicefrac {1}{2}} \Big ) \nabla _h \Big ( G_{\hat{{\textsf{N}}}_2^m}^1 \Big ( \hat{{\textsf{N}}}_2^{m+1} \Big ) + (-\Delta _h)^{-1} \Big ( \hat{{\textsf{N}}}_2^{m+\nicefrac {1}{2}} - \hat{\mathsf P}_2^{m+\nicefrac {1}{2}} \Big ) \Big ) \nonumber \\{} & {} \, \, + \Big ( \frac{3}{2} \hat{{\textsf{N}}}_2^m - \frac{1}{2} \hat{\mathsf N}_2^{m-1} \Big ) \nabla _h \Big ( \frac{1}{\hat{{\textsf{N}}}_2^{m+\nicefrac {1}{2}} } {{\textsf{N}}}_{h, 1}^{m+\nicefrac {1}{2}} + (-\Delta _h)^{-1} \Big ( {{\textsf{N}}}_{h, 1}^{m+\nicefrac {1}{2}} - {{\textsf{P}}}_{h, 1}^{m+\nicefrac {1}{2}} \Big ) \Big ) \Big ) \nonumber \\{} & {} \quad - \big ( H_n^{(0)} \big )^{m+\nicefrac {1}{2}} + O \big ({\Delta t}^2 + h^2\big ), \end{aligned}$$

(A.20)

$$\begin{aligned} \frac{{{\textsf{P}}}_{h, 1}^{m+1} - {{\textsf{P}}}_{h, 1}^m}{{\Delta t}}= & {} \nabla _h \cdot \big ( D \big ( \breve{{\textsf{P}}}_{h, 1}^{m+\nicefrac {1}{2}} \big ) \nabla _h \big ( G_{\hat{\mathsf P}_2^m}^1 \big ( \hat{{\textsf{P}}}_2^{m+1} \big ) + (-\Delta _h)^{-1} \big ( \hat{{\textsf{P}}}_2^{m+\nicefrac {1}{2}} - \hat{{\textsf{N}}}_2^{m+\nicefrac {1}{2}} \big ) \big )\nonumber \\{} & {} \, \, + D \big ( \frac{3}{2} \hat{{\textsf{P}}}_2^m - \frac{1}{2} \hat{{\textsf{P}}}_2^{m-1} \big ) \nabla _h \big ( \frac{1}{\hat{\mathsf P}_2^{m+\nicefrac {1}{2}} } {{\textsf{P}}}_{h, 1}^{m+\nicefrac {1}{2}} + (-\Delta _h)^{-1} \big ( {{\textsf{P}}}_{h, 1}^{m+\nicefrac {1}{2}} - {{\textsf{P}}}_{h, 1}^{m+\nicefrac {1}{2}} \big ) \big ) \big )\nonumber \\{} & {} \quad - \big ( H_p^{(0)} \big )^{m+\nicefrac {1}{2}} + O \big ({\Delta t}^2 + h^2 \big ), \nonumber \\ \breve{{\textsf{N}}}_{h, 1}^{m+\nicefrac {1}{2}}= & {} \frac{3}{2} {{\textsf{N}}}_{h, 1}^m - \frac{1}{2} {{\textsf{N}}}_{h, 1}^{m-1}, \quad \breve{{\textsf{P}}}_{h, 1}^{m+\nicefrac {1}{2}} = \frac{3}{2} {{\textsf{P}}}_{h, 1}^m - \frac{1}{2} {{\textsf{P}}}_{h, 1}^{m-1}. \end{aligned}$$

(A.21)

Finally, a combination of (A.18)–(A.19) and (A.20)–(A.21) yields the higher order truncation error for $(\hat{{\textsf{N}}}, \hat{\mathsf P})$, as given by (4.10)–(4.11). Of course, the linear expansions have been extensively utilized.

Moreover, we see that trivial initial data ${{\textsf{N}}}_{{\Delta t}, j} (\, \cdot \,, t=0), {{\textsf{P}}}_{{\Delta t}, j} (\, \cdot \,, t=0) \equiv 0$ could be taken ($j=1, 2$) as in (A.3)–(A.4), (A.10)–(A.11), respectively, as well as $({{\textsf{N}}}_{h, 1}, {{\textsf{P}}}_{h, 1})$ in (A.18)–(A.19). Consequently, using similar arguments as in (4.4)–(4.5), we arrive at the mass conservative identities (4.12), (4.13). Notice that the first step of (4.13) is based on the fact that $\hat{{\textsf{N}}} \in {{{\mathcal {B}}}}^K$, and the second step comes from the mass conservative property of $\hat{{\textsf{N}}}$ at the continuous level. And also, the mass conservative property of $(\hat{{\textsf{N}}}, \hat{{\textsf{P}}})$ is stated in (4.13), and we conclude that the local truncation error $\tau _n$, $\tau _p$ has a similar property, so that (4.14) is proved.

Based on the fact that the temporal and spatial correction functions $({{\textsf{N}}}_{{\Delta t}, j}, {{\textsf{P}}}_{{\Delta t}, j})$, $({{\textsf{N}}}_{h, 1}, {{\textsf{P}}}_{h, 1})$ are bounded, we recall the separation property (4.2) for the exact solution. In turn, a similar property (4.15) becomes available for the constructed profile $(\hat{{\textsf{N}}}, \hat{{\textsf{P}}})$, in which the projection estimate (4.3) has been recalled. Of course, ${\Delta t}$ and h could be taken sufficiently small so that (4.15) is valid for a modified value $\epsilon _0^\star$, such as $\epsilon _0^\star = \frac{1}{4} \epsilon _0$.

Furthermore, we recall the fact that the correction functions stay bounded, in terms of both the spatial and temporal derivatives, since they only depend on $({{\textsf{N}}}_N, {{\textsf{P}}}_N)$ and the exact solution. Therefore, a discrete $W^{1,\infty }$ bound for $(\hat{{\textsf{N}}}, \hat{{\textsf{P}}})$ could be derived as in (4.16), as well as a bound (4.17) for its discrete temporal derivative. This completes the proof of Proposition 4.1.

Appendix B: Proof of Lemma 4.1

For the nonlinear inner product associated with the artificial regularization, the following fact is observed:

$$\begin{aligned} \begin{aligned}&\ln \hat{{\textsf{N}}}^{m+1} - \ln n^{m+1} = \frac{1}{\zeta _n^{(m+1)}} {\tilde{n}}^{m+1}, \quad \zeta _n^{(m+1)} \text{ between } \hat{{\textsf{N}}}^{m+1} \text{ and } n^{m+1}, \quad \text{(Taylor } \text{ expansion) },\\&\bigg \langle \ln \hat{{\textsf{N}}}^{m+1} - \ln n^{m+1}, {\tilde{n}}^{m+1} \bigg \rangle = \bigg \langle \frac{1}{\zeta _n^{(m+1)}} {\tilde{n}}^{m+1}, {\tilde{n}}^{m+1} \bigg \rangle \ge 0, \end{aligned} \end{aligned}$$

(B.1)

in which the positivity-preserving property of $n^{m+1}$ and $\check{{\textsf{N}}}^{m+1}$ has been applied. For the additional term in the artificial regularization part, we have the following estimates:

$$\begin{aligned}&\ln \hat{{\textsf{N}}}^m - \ln n^m = \frac{1}{\zeta _n^{(m)}} {\tilde{n}}^m , \quad \zeta _n^{(m)} \text{ between } \hat{\mathsf N}^m \text{ and } n^m , \quad \text{(Taylor } \text{ expansion) } , \end{aligned}$$

(B.2)

$$\begin{aligned}&| \frac{1}{\zeta _n^{(m)}} | \le \max \Big ( \frac{1}{\hat{\mathsf N}^m} , \frac{1}{n^m} \Big ) \le \frac{2}{\epsilon _0^*} , \quad (\text{ by }~(4.15), (4.30)) , \end{aligned}$$

(B.3)

$$\begin{aligned}&| \ln \hat{{\textsf{N}}}^m - \ln n^m | = | \frac{1}{\zeta _n^{(m)}} | \cdot | {\tilde{n}}^m | \le \frac{2}{\epsilon _0^*} | {\tilde{n}}^m | , \quad | \langle \ln \hat{{\textsf{N}}}^m - \ln n^m , {\tilde{n}}^{m+1} \rangle | \le \frac{2}{\epsilon _0^*} | \langle {\tilde{n}}^m , {\tilde{n}}^{m+1} \rangle | . \end{aligned}$$

(B.4)

For the error term $\langle {\tilde{n}}^{m+1}, G_{\hat{\mathsf N}^m}^1 ( \hat{{\textsf{N}}}^{m+1} ) - G_{n^m}^1 ( n^{m+1} ) \rangle$, we begin with the following decomposition:

$$\begin{aligned} G_{\hat{{\textsf{N}}}^m}^1 \big ( \hat{{\textsf{N}}}^{m+1}\big ) - G_{n^m}^1 \big ( n^{m+1} \big ){} & {} = \mathcal{NLE}_{1}^m + \mathcal{NLE}_{2}^m, \quad \mathcal{NLE}_{1}^m = G_{n^m}^1 \big ( \hat{{\textsf{N}}}^{m+1} \big ) - G_{n^m}^1 \big ( n^{m+1} \big ), \nonumber \\ \mathcal{NLE}_{2}^m{} & {} = G_{\hat{{\textsf{N}}}^{m+1} }^1 \big ( \hat{\mathsf N}^m \big ) - G_{ \hat{{\textsf{N}}}^{m+1} }^1 \big ( n^m \big ), \nonumber \\ \end{aligned}$$

(B.5)

in which $G_a^1 (x)$ has been introduced in (3.1). For the $\mathcal{NLE}_{2}^m$ error term, we apply the intermediate value theorem and obtain

$$\begin{aligned} \begin{aligned} \mathcal{NLE}_{2}^m&= G_{\hat{{\textsf{N}}}^{m+1} }^1 \big ( \hat{{\textsf{N}}}^m \big ) - G_{ \hat{{\textsf{N}}}^{m+1} }^1 ( n^m ) \\&= \big ( G_{\hat{\mathsf N}^{m+1} }^1 \big )' \big (\eta _n^{(m)} \big ) \big ( \hat{{\textsf{N}}}^m - n^m \big ), \, \, \eta _n^{(m)} \text{ between } \hat{{\textsf{N}}}^m \text{ and } n^m. \end{aligned} \end{aligned}$$

(B.6)

Meanwhile, by property (3) in Lemma 3.1, we conclude that

$$\begin{aligned} \big ( G_{\hat{{\textsf{N}}}^{m+1} }^1 \big )' \big (\eta _n^{(m)} \big ) = \frac{1}{2 \xi _n^{(m)} }, \quad \xi _n^{(m)} \text{ between } \hat{\mathsf N}^{m+1} \text{ and } \eta _n^{(m)}. \end{aligned}$$

(B.7)

By the combined fact that, $\xi _n^{(m)}$ is between $\hat{\mathsf N}^{m+1}$ and $\eta _n^{(m)}$, $\eta$ is between $\hat{{\textsf{N}}}^m$ and $n^m$, the following bound is available

$$\begin{aligned} \xi _n^{(m)} \text{ is } \text{ between } \text{ the } \text{ values } \text{ of } \hat{\mathsf N}^{m+1}, \hat{{\textsf{N}}}^m \text{ and } n^m, \, \, \, \Big | \frac{1}{\xi _n^{(m)}} \Big | \le \max \Big ( \frac{1}{\hat{\mathsf N}^{m+1} }, \frac{1}{\hat{{\textsf{N}}}^m}, \frac{1}{n^m} \Big ) \le \frac{2}{\epsilon _0^*}, \nonumber \\ \end{aligned}$$

(B.8)

in which the phase separation properties (4.15), (4.30) have been recalled. Then we arrive at

$$\begin{aligned} \big | \mathcal{NLE}_{2}^m \big | = \big | \big ( G_{\hat{{\textsf{N}}}^{m+1} }^1 \big )' \big (\eta _n^{(m)} \big ) \big | \cdot \big | {\tilde{n}}^m \big | = \Big | \frac{1}{2 \xi _n^{(m)} } \Big | \cdot \big | {\tilde{n}}^m \big | \le \frac{1}{\epsilon _0^*} \big | {\tilde{n}}^m \big |. \end{aligned}$$

(B.9)

Based on this point-wise bound, the following inequality is available

$$\begin{aligned} \langle \mathcal{NLE}_{2}^m, {\tilde{n}}^{m+1} \rangle \ge - \frac{1}{\epsilon _0^*} | \langle {\tilde{n}}^m, {\tilde{n}}^{m+1} \rangle | \ge - \frac{h^3}{\epsilon _0^*} \sum _{i,j,k} | {\tilde{n}}^m_{i,j,k} | \cdot | {\tilde{n}}^{m+1}_{i,j,k} |. \end{aligned}$$

(B.10)

For the $\mathcal{NLE}_{1}^m$ error term, a similar nonlinear analysis could be performed:

$$\begin{aligned} \begin{aligned}&\mathcal{NLE}_{1}^m = G_{n^m}^1 \big ( \hat{{\textsf{N}}}^{m+1} \big ) - G_{n^m}^1 \big ( n^{m+1} \big ) = \big ( G_{\hat{n^m} }\big )' \big (\eta _n^{(m+1)} \big ) \big ( \hat{\mathsf N}^{m+1} - n^{m+1} \big ),\\&\big ( G_{n^m }^1 \big )' \big (\eta _n^{(m+1)} \big ) = \frac{1}{2 \xi _n^{(m+1)} }, \quad \xi _n^{(m+1)} \text{ between } n^m \text{ and } \eta _n^{(m+1)},\\&\xi _n^{(m+1)} \text{ is } \text{ between } \text{ the } \text{ values } \text{ of } n^m, \hat{{\textsf{N}}}^{m+1} \text{ and } n^{m+1}, \, \, \, \frac{1}{\xi _n^{(m+1)}} \ge \min \Big ( \frac{1}{\hat{\mathsf N}^{m+1} }, \frac{1}{n^{m+1}}, \frac{1}{n^m} \Big ). \end{aligned}\nonumber \\ \end{aligned}$$

(B.11)

This in turn leads to

$$\begin{aligned} \Bigg \langle \mathcal{NLE}_{1}^m, {\tilde{n}}^{m+1} \Bigg \rangle = \Bigg \langle \frac{1}{2 \xi _n^{(m+1)} }, \big ( {\tilde{n}}^{m+1} \big )^2 \Bigg \rangle = \frac{1}{2} h^3 \sum _{i,j,k} \frac{1}{ \big ( \xi _n^{(m+1)} \big )_{i,j,k} } \cdot | {\tilde{n}}^{m+1}_{i,j,k} |^2. \end{aligned}$$

(B.12)

Subsequently, a combination of (B.10) and (B.12) yields

$$\begin{aligned}{} & {} \Bigg \langle {\tilde{n}}^{m+1}, G_{\hat{{\textsf{N}}}^m}^1 \big ( \hat{{\textsf{N}}}^{m+1} \big ) - G_{n^m}^1 \big ( n^{m+1} \big ) \Bigg \rangle \nonumber \\{} & {} \quad \ge h^3 \sum _{i,j,k} \Big ( \frac{\frac{1}{2}}{ \big ( \xi _n^{(m+1)} \big )_{i,j,k} } \cdot \big | {\tilde{n}}^{m+1}_{i,j,k} \big |^2 - ( \epsilon _0^*)^{-1} \big | {\tilde{n}}^m_{i,j,k} \big | \cdot \big | {\tilde{n}}^{m+1}_{i,j,k} \big | \Big ). \end{aligned}$$

(B.13)

Moreover, its combination with (B.1)–(B.4) indicates that

$$\begin{aligned} \begin{aligned}&{\Delta t}\Bigg \langle {\tilde{n}}^{m+1}, \ln \hat{{\textsf{N}}}^{m+1} - \ln n^{m+1} - \big ( \ln \hat{{\textsf{N}}}^m - \ln n^m \big )\Bigg \rangle + \Bigg \langle {\tilde{n}}^{m+1}, G_{\hat{{\textsf{N}}}^m}^1 ( \hat{{\textsf{N}}}^{m+1} ) - G_{n^m}^1 ( n^{m+1} ) \Bigg \rangle \\&\quad \ge h^3 \sum _{i,j,k} \Bigg ( \frac{\frac{1}{2}}{ ( \xi _n^{(m+1)} )_{i,j,k} } \cdot | {\tilde{n}}^{m+1}_{i,j,k} |^2 - \frac{3}{2} ( \epsilon _0^*)^{-1} | {\tilde{n}}^m_{i,j,k} | \cdot | {\tilde{n}}^{m+1}_{i,j,k} | \Bigg ). \end{aligned} \end{aligned}$$

(B.14)

At each fixed grid point (i, j, k), if (i, j, k) is not in $\Lambda _n$, i.e., $0< n^{m+1}_{i,j,k} < 2C^* +1$, the following estimates are available:

$$\begin{aligned} \begin{aligned}&\frac{1}{ n^{m+1}_{i,j,k} } \ge \frac{1}{2 C^* +1}, \, \, \frac{1}{\hat{{\textsf{N}}}^{m+1}_{i,j,k} } \ge \frac{1}{C^*}, \, \, \frac{1}{ n^m_{i,j,k} } \ge \frac{1}{C^* +1}, \quad (\text{ by }~(4.16), (4.29)),\\&\text{ so } \text{ that } \, \, \, \frac{1}{\big (\xi _n^{(m+1)} \big )_{i,j,k} } \ge \min \Bigg ( \frac{1}{\hat{{\textsf{N}}}^{m+1}_{i,j,k} }, \frac{1}{n^{m+1}_{i,j,k} }, \frac{1}{n^m_{i,j,k} } \Bigg ) \ge \frac{1}{2 C^* +1}. \end{aligned} \end{aligned}$$

(B.15)

In turn, the following inequality is valid for $0< n^{m+1}_{i,j,k} < 2C^* +1$:

$$\begin{aligned} \begin{aligned}&\frac{\frac{1}{2}}{ \big ( \xi _n^{(m+1)} \big )_{i,j,k} } \cdot | {\tilde{n}}^{m+1}_{i,j,k} |^2 - \frac{3}{2} \big ( \epsilon _0^*\big )^{-1} | {\tilde{n}}^m_{i,j,k} | \cdot | {\tilde{n}}^{m+1}_{i,j,k} | \\&\quad \ge \, \frac{\frac{1}{2}}{2 C^* +1} | {\tilde{n}}^{m+1}_{i,j,k} |^2 - \frac{3}{2} \big ( \epsilon _0^*\big )^{-1} | {\tilde{n}}^m_{i,j,k} | \cdot | {\tilde{n}}^{m+1}_{i,j,k} |\\&\quad \ge \frac{1}{4 C^* +2} | {\tilde{n}}^{m+1}_{i,j,k} |^2 - \frac{1}{8 C^* +4} | {\tilde{n}}^{m+1}_{i,j,k} |^2\\&\qquad - \frac{9 (2 C^* +1)}{16} ( \epsilon _0^*)^{-2} | {\tilde{n}}^m_{i,j,k} |^2\\&\quad \ge \frac{1}{8 C^* +4} | {\tilde{n}}^{m+1}_{i,j,k} |^2 - \frac{9 (2 C^* +1)}{16} ( \epsilon _0^*)^{-2} | {\tilde{n}}^m_{i,j,k} |^2. \end{aligned} \end{aligned}$$

(B.16)

At each fixed grid point (i, j, k), if $(i, j, k) \in \Lambda _n$, i.e., $n^{m+1}_{i,j,k} \ge 2C^* +1$, we see that $n^{m+1}_{i,j,k} > \max ( \hat{{\textsf{N}}}^{m+1}_{i,j,k}, n^m_{i,j,k} )$, so that

$$\begin{aligned} \frac{1}{\big (\xi _n^{(m+1)} \big )_{i,j,k} } \ge \min \Big ( \frac{1}{\hat{{\textsf{N}}}^{m+1}_{i,j,k} }, \frac{1}{n^{m+1}_{i,j,k} }, \frac{1}{n^m_{i,j,k} } \Big ) = \frac{1}{n^{m+1}_{i,j,k} }. \end{aligned}$$

(B.17)

Meanwhile, since $n^{m+1}_{i,j,k} > \hat{{\textsf{N}}}^{m+1}_{i,j,k}$, we see that ${\tilde{n}}^{m+1}_{i,j,k} = \hat{{\textsf{N}}}^{m+1}_{i,j,k} - n^{m+1}_{i,j,k} < 0$, and the following fact is observed:

$$\begin{aligned} \begin{aligned}&\hat{{\textsf{N}}}^{m+1}_{i,j,k} \le C^* \le \frac{C^*}{2 C^* +1} ( 2 C^* +1) \le \frac{C^* n^{m+1}_{i,j,k} }{2 C^* +1},\\&| {\tilde{n}}^{m+1} | = | \hat{{\textsf{N}}}^{m+1}_{i,j,k} - n^{m+1}_{i,j,k} | \ge | n^{m+1}_{i,j,k} | - \frac{C^* n^{m+1}_{i,j,k} }{2 C^* +1} \ge \frac{C^* +1 }{2 C^* +1} n^{m+1}_{i,j,k},\\&\frac{1}{(\xi _n^{(m+1)} )_{i,j,k} } \cdot | {\tilde{n}}^{m+1}_{i,j,k} | \ge \frac{1}{n^{m+1}_{i,j,k} } \cdot \frac{C^* +1 }{2 C^* +1} n^{m+1}_{i,j,k} = \frac{C^* +1 }{2 C^* +1} = \frac{1}{2} + \frac{\frac{1}{2}}{2 C^* +1}. \end{aligned} \nonumber \\ \end{aligned}$$

(B.18)

Subsequently, the following inequality could be derived, for $n^{m+1}_{i,j,k} \ge 2 C^* +1$:

$$\begin{aligned} \begin{aligned}&\frac{\frac{1}{2}}{ \big ( \xi _n^{(m+1)} \big )_{i,j,k} } \cdot | {\tilde{n}}^{m+1}_{i,j,k} |^2 - \frac{3}{2} ( \epsilon _0^*)^{-1} | {\tilde{n}}^m_{i,j,k} | \cdot | {\tilde{n}}^{m+1}_{i,j,k} |\\&\quad \ge \frac{\frac{1}{2} (C^* + 1) }{2 C^* +1} \cdot | {\tilde{n}}^{m+1}_{i,j,k} | - \frac{3}{2} \big ( \epsilon _0^*\big )^{-1} \cdot {\Delta t}^2 \cdot | {\tilde{n}}^{m+1}_{i,j,k} | \ge \Big ( \frac{\frac{1}{2} (C^* + 1) }{2 C^* +1} - \frac{3}{2} ( \epsilon _0^*)^{-1} {\Delta t}^2 \Big ) | {\tilde{n}}^{m+1}_{i,j,k} |\\&\quad \ge \frac{1}{4} { | {\tilde{n}}^{m+1}_{i,j,k} | } \ge \frac{1}{4} ( 2 C^* +1) \ge \frac{1}{2} C^*, \end{aligned} \nonumber \\ \end{aligned}$$

(B.19)

in which the $\Vert \cdot \Vert _\infty$ estimate (4.28) has been applied in the second step, and the fact the $\frac{\frac{1}{2} (C^* + 1) }{2 C^* +1} - \frac{3}{2} ( \epsilon _0^*)^{-1} {\Delta t}^2 \ge \frac{1}{4}$ has been used in the fourth step.

Consequently, a substitution of the point-wise inequalities (B.16), (B.19) into (B.14) results in the desired estimate (4.33), by taking ${\tilde{C}}_4 = \frac{9 (2 C^* +1)}{16} ( \epsilon _0^*)^{-2}$. The other inequality in (4.33) could be derived in the same manner; the technical details are skipped for the sake of brevity.

Finally, if $K_n^* =0$ and $K_p^*=0$, i.e, both $\Lambda _n$ and $\Lambda _p$ are empty sets, we see that inequality (B.16) is valid at every grid points. In turn, a summation in space results in the improved nonlinear estimate (4.34). This finishes the proof of Lemma 4.1.

Appendix C: Proof of Lemma 4.2

The Taylor expansions (B.6)–(B.8), (B.11) are still valid, which in turn imply that

$$\begin{aligned} G_{\hat{{\textsf{N}}}^m}^1 \big ( \hat{{\textsf{N}}}^{m+1} \big ) - G_{n^m}^1 \big ( n^{m+1} \big ) = \mathcal{NLE}_1^m + \mathcal{NLE}_{2}^m = \frac{1}{2 \xi _n^{(m+1)} } {\tilde{n}}^{m+1} + \frac{1}{2 \xi _n^{(m)} } {\tilde{n}}^m. \end{aligned}$$

(C.1)

In particular, $\xi _n^{(m)}$ could be analyzed in a more precise way:

$$\begin{aligned} \begin{aligned}&\xi _n^{(m)} \text{ is } \text{ between } \text{ the } \text{ values } \text{ of } \hat{\mathsf N}^{m+1}, \hat{{\textsf{N}}}^m \text{ and } n^m, \text{ so } \text{ that }\\&\quad \max \big ( | \xi _n^{(m)} - \hat{{\textsf{N}}}^{m+1} |, | \xi _n^{(m)} - \hat{{\textsf{N}}}^m | \big ) \\&\qquad \le | \hat{{\textsf{N}}}^{m+1} - \hat{\mathsf N}^m | + \big | n^m - \hat{{\textsf{N}}}^m \big | \le C^* {\Delta t}+ {\Delta t}= (C^* +1 ) {\Delta t}, \end{aligned} \end{aligned}$$

(C.2)

in which the regularity requirement (4.16) for the constructed profile $\hat{{\textsf{N}}}$, as well as the $\Vert \cdot \Vert _\infty$ a-priori error estimate for ${\tilde{n}}^m$, has been applied. A similar error bound could also be derived for $\xi _n^{(m+1)}$, with the technical details skipped for the sake of brevity:

$$\begin{aligned} \max \big ( \big | \xi _n^{(m+1)} - \hat{{\textsf{N}}}^{m+1} \big |, \big | \xi _n^{(m+1)} - \hat{{\textsf{N}}}^m \big | \big )\le & {} \,| \hat{{\textsf{N}}}^{m+1} - \hat{\mathsf N}^m | + | n^{m+1} - \hat{{\textsf{N}}}^{m+1} | + | n^m - \hat{\mathsf N}^m | \nonumber \\\le & {} \, C^* {\Delta t}+ {\Delta t}= \big (C^* +1 \big ) {\Delta t}. \end{aligned}$$

(C.3)

Then we arrive at the point-wise error bounds:

$$\begin{aligned} \begin{aligned} \Big | \frac{1}{2 \xi _n^{(m)} } - \frac{1}{\hat{{\textsf{N}}}^{m+1} + \hat{{\textsf{N}}}^m } \Big | =&\, \frac{ | \hat{{\textsf{N}}}^{m+1} - \xi _n^{(m)} + \hat{{\textsf{N}}}^m - \xi _n^{(m)} | }{ 2 \xi _n^{(m)} ( \hat{{\textsf{N}}}^{m+1} + \hat{{\textsf{N}}}^m ) }\\ \le&\, \big ( \epsilon _0^* \big )^{-1} \cdot 2 \big ( \epsilon _0^* \big )^{-1} \cdot 2 (C^* +1 ) {\Delta t}= M^{(0)} {\Delta t},\\ \Big | \frac{1}{2 \xi _n^{(m+1)} } - \frac{1}{\hat{{\textsf{N}}}^{m+1} + \hat{{\textsf{N}}}^m } \Big | \le&\, M^{(0)} {\Delta t}, \, \, \, \text{(similar } \text{ analysis) }, \end{aligned} \end{aligned}$$

(C.4)

with $M^{(0)} = 4 ( \epsilon _0^* )^{-2} \cdot (C^* +1 )$, and the phase separation estimates (4.30), (4.31), (4.51) have been applied. As a further consequence, the following inequality could be derived

$$\begin{aligned} \begin{aligned}&\langle {\tilde{n}}^{m+1} - {\tilde{n}}^m, G_{\hat{{\textsf{N}}}^m}^1 \big ( \hat{{\textsf{N}}}^{m+1} \big ) - G_{n^m}^1 ( n^{m+1} ) \rangle = \Bigg \langle {\tilde{n}}^{m+1} - {\tilde{n}}^m, \frac{1}{2 \xi _n^{(m+1)} } {\tilde{n}}^{m+1} + \frac{1}{2 \xi _n^{(m)} } {\tilde{n}}^m \Bigg \rangle \\&\quad \ge \Bigg \langle {\tilde{n}}^{m+1} - {\tilde{n}}^m, \frac{1}{\hat{\mathsf N}^{m+1} + \hat{{\textsf{N}}}^m } \big ( {\tilde{n}}^{m+1} + {\tilde{n}}^m \big ) \Bigg \rangle - \Bigg \langle | {\tilde{n}}^{m+1} - {\tilde{n}}^m |, M^{(0)} {\Delta t}\big ( | {\tilde{n}}^{m+1} | + | {\tilde{n}}^m | \big ) \Bigg \rangle \\&\quad \ge \Big \langle \frac{1}{\hat{{\textsf{N}}}^{m+1} + \hat{{\textsf{N}}}^m }, \big ( {\tilde{n}}^{m+1} \big )^2 \Big \rangle - \Big \langle \frac{1}{\hat{{\textsf{N}}}^{m+1} + \hat{{\textsf{N}}}^m }, \big ( {\tilde{n}}^m \big )^2 \Big \rangle - 2 M^{(0)} {\Delta t}\big ( \Vert {\tilde{n}}^{m+1} \Vert _2^2 + \Vert {\tilde{n}}^m \Vert _2^2 \big ). \end{aligned} \nonumber \\ \end{aligned}$$

(C.5)

On the other hand, the following estimates are straightforward:

$$\begin{aligned} \begin{aligned} \Big | \frac{1}{\hat{{\textsf{N}}}^{m+1} + \hat{{\textsf{N}}}^m } - \frac{1}{2 \hat{{\textsf{N}}}^{m+1} } \Big | =&\,\frac{ | \hat{\mathsf N}^{m+1} - \hat{{\textsf{N}}}^m | }{ 2 \hat{{\textsf{N}}}^{m+1} ( \hat{{\textsf{N}}}^{m+1} + \hat{{\textsf{N}}}^m ) } \le C^* {\Delta t}\cdot 4 ( \epsilon _0^* )^{-2} = M^{(1)} {\Delta t},\\ \Big | \frac{1}{\hat{{\textsf{N}}}^{m+1} + \hat{{\textsf{N}}}^m } - \frac{1}{2 \hat{{\textsf{N}}}^m } \Big | \le&\, M^{(1)} {\Delta t}, \end{aligned} \end{aligned}$$

(C.6)

with $M^{(1)} = 4 C^* ( \epsilon _0^* )^{-2}$, in which the temporal regularity assumptions (4.17) for the constructed profile $\hat{{\textsf{N}}}$, as well as the phase separation estimate (4.30), have been applied. Subsequently, the following inequalities are obvious:

$$\begin{aligned} \begin{aligned} \Bigg \langle \frac{1}{\hat{{\textsf{N}}}^{m+1} + \hat{{\textsf{N}}}^m }, ( {\tilde{n}}^{m+1} )^2 \Bigg \rangle - \frac{1}{2} \Bigg \langle \frac{1}{\hat{{\textsf{N}}}^{m+1} }, ( {\tilde{n}}^{m+1} )^2 \Bigg \rangle \ge&- M^{(1)} {\Delta t}\Vert {\tilde{n}}^{m+1} \Vert _2^2,\\ \Bigg \langle \frac{1}{\hat{{\textsf{N}}}^{m+1} + \hat{{\textsf{N}}}^m }, \big ( {\tilde{n}}^m \big )^2 \Bigg \rangle - \frac{1}{2} \Bigg \langle \frac{1}{\hat{{\textsf{N}}}^m }, \big ( {\tilde{n}}^m \big )^2 \Bigg \rangle \le&\, M^{(1)} {\Delta t}\Vert {\tilde{n}}^m \Vert _2^2. \end{aligned} \end{aligned}$$

(C.7)

Finally, a combination of (C.5) and (C.7) results in the refined estimate (4.53), by taking ${\tilde{C}}_6 = 2\,M^{(0)} + M^{(1)}$. Since both $M^{(0)}$ and $M^{(1)}$ only depend on $C^*$, $\epsilon _0^*$, the same dependence preserves for ${\tilde{C}}_6$.

Inequalities (4.54), (4.55), and (4.56) could be derived in a similar manner; the technical details are skipped for the sake of brevity.