Error Analysis with Polynomial Dependence on $\varepsilon ^{-1}$ in SAV Methods for the Cahn-Hilliard Equation

Abstract

The optimal error estimate that depends only on the polynomial degree of $\varepsilon ^{-1}$ is established for the temporal semi-discrete scheme of the Cahn-Hilliard equation based on the scalar auxiliary variable (SAV) formulation. The key to our analysis is converting the structure of the SAV time-stepping scheme back to a form compatible with the original format of the Cahn-Hilliard equation, which makes it feasible to use spectral estimates to handle the nonlinear term. Based on the transformation of the SAV numerical scheme, the optimal error estimate for the temporal semi-discrete scheme which depends only on the low polynomial order of $\varepsilon ^{-1}$ instead of the exponential order, is derived by using mathematical induction, spectral arguments, and the superconvergence properties of some nonlinear terms. Numerical examples are provided to illustrate the discrete energy decay property and validate our theoretical convergence analysis.

1 Introduction

In this paper we consider the initial boundary value problem for the Cahn-Hilliard (CH) phase field equation

$$\begin{aligned} \partial _t u&= \Delta \big ( - \varepsilon \Delta u + \frac{1}{\varepsilon } f(u) \big ) & \qquad \text{ in }\,\,\,{\varOmega }\times (0,T], \end{aligned}$$

(1.1a)

$$\begin{aligned} \partial _{\varvec{n}} u&= \partial _{\varvec{n}} \big ( - \varepsilon \Delta u + \frac{1}{\varepsilon } f(u) \big ) = 0 & \qquad \text{ on }\,\,\, \partial {\varOmega }\times (0,T],\end{aligned}$$

(1.1b)

$$\begin{aligned} u(\cdot , 0)&= u_0, & \qquad \text{ in }\,\,\,{\varOmega }, \end{aligned}$$

(1.1c)

where ${\varOmega }\subset R^d$, $d = 2, 3$ is a bounded domain, ${\varvec{n}}$ is the outward normal, $\varepsilon$ is a small parameter, and f is the derivative of a non-negative potential function F with two local minima, i.e., $f = F'$. For instance, the Ginzburg-Landau energy function

$$ F(v) = \frac{1}{4} (v^2 - 1)^2 \quad \text{ and } \quad f(v) = v^3 - v. $$

In view of its wide application as a phase field model [10, 39], many numerical methods and analyses have been developed to approximate the Cahn-Hilliard equation (1.1). However, most of these methods have been developed for the Cahn-Hilliard equation with fixed $\varepsilon >0$. As pointed out in [4, 15, 17, 18], error estimates using the direct Gronwall inequality argument yield a constant factor $e^{T/\varepsilon }$, resulting in exponential growth of the error as $\varepsilon \rightarrow 0$. Such an estimate is clearly not useful for very small value of $\varepsilon$, especially in solving the problem of whether the computed numerical interface converges to the original sharp interface of the Hele-Shaw problem when $\varepsilon \rightarrow 0$ (see [14, 28] for details). To overcome this difficulty, Feng and Prohl [25] first established a priori error estimates with polynomial dependence on $\varepsilon ^{-1}$ for the time-discrete format of the Cahn-Hilliard equation. Subsequently, Feng and Wu obtained a posteriori error estimates with polynomial-order of $\varepsilon ^{-1}$ in [26] for the same time-discrete methods. The main idea of the polynomial-order error estimates of $\varepsilon ^{-1}$ is to use the spectral estimates given by Alikakos and Fusco [3] and Chen [13] for the linearized Cahn-Hilliard operator to handle the nonlinear term in the error analysis. After this, spectral estimates have been frequently used to eliminate the exponential dependence on $\varepsilon ^{-1}$ in the error analyses of other numerical methods for the Cahn-Hilliard equations (see [6, 19, 21, 29, 34] and references therein) and the related phase field equations, including the Allen-Cahn equations in [1, 6,7,8, 20, 23, 24, 27], the Ginzburg-Landau equations in [5], and the phase field models with nonlinear constitutive laws in [16].

On the other hand, it is well known that the Cahn-Hilliard equation (1.1) is the $H^{-1}$-gradient flow of the energy functional

$$\begin{aligned} E[u] := \int _{{\varOmega }} \Big ( \frac{\varepsilon }{2} |\nabla u |^2 + \frac{1}{\varepsilon } F(u) \Big ) \, {\mathrm d}x. \end{aligned}$$

(1.2)

As a result, the solution of the Cahn-Hilliard equation has decaying energy. Indeed, testing (1.1) by $- \varepsilon \Delta u + \frac{1}{\varepsilon } f(u)$ yields

$$\begin{aligned} E[u(\cdot , t_2)] - E[u(\cdot , t_1)] = - \int _{t_1}^{t_2} \Vert \partial _t u \Vert _{H^{-1}}^2 \, {\mathrm d}t \le 0, \end{aligned}$$

(1.3)

where $0 \le t_1 \le t_2 \le T$ and $\Vert \partial _t u\Vert _{H^{-1}} = \Vert \nabla w\Vert$ with $w = - \varepsilon \Delta u + \frac{1}{\varepsilon } f(u)$.

Accordingly, great efforts have been devoted to the construction of efficient and accurate numerical methods that preserve the energy decay properties at the discrete level. In particular, for widely used linear time-stepping methods with energy decay, including stabilized semi-implicit schemes [11, 12, 33], invariant energy quadratization (IEQ) methods [35,36,37], and the scalar auxiliary variable (SAV) approach in [30, 31] and [2], optimal error estimates of these schemes have been well established for the Cahn-Hilliard equations with fixed $\varepsilon ^{-1}$. However, the main challenge with error analysis in these methods is how to establish error bounds that depend only on the polynomial order of $\varepsilon ^{-1}$ rather than the exponential order for small $\varepsilon \rightarrow 0$. This difficulty arises because, unlike in Feng’s previous work [19, 21, 25, 26, 29], the numerical methods based on the IEQ/SAV formulation break the standard structure of the nonlinear term of the Cahn-Hilliard equation, which is crucial to the utilization of the spectral arguments. As a result, spectral estimates become ineffective in estimating errors for the IEQ/SAV approach. To the best of our knowledge, optimal error estimates that depend only on the polynomial order of $\varepsilon ^{-1}$ for the IEQ/SAV methods to the Cahn-Hilliard equation remain open.

The objective of this paper is to establish error bounds that depend only on low polynomial order of $\varepsilon ^{-1}$ for semi-discrete methods based on the SAV formulation of the Cahn-Hilliard equation. The SAV reformulation of the Cahn-Hilliard equation was introduced in [31, 32] as an enhanced version of the invariant energy quadratization (IEQ) approach [35,36,37,38], for developing energy-decay methods at the discrete level. By reconstructing the system based on the SAV reformulation, a linear and easy-to-implement time semi-discrete scheme is obtained. The SAV formulation introduces new difficulties to the error analysis for the Cahn-Hilliard equation due to the presence of a new scalar r in the nonlinear part (see equation (2.5b)), which alters the structure of the original Cahn-Hilliard equation and makes the spectral argument used in [19, 21, 25, 26, 29] not directly applicable. To overcome this problem, we transform the structure of the SAV scheme using the new scalar variable into a form that is compatible with the original format of the Cahn-Hilliard equation (1.1). This transformation enable us to use spectral estimates in our error analysis. However, this transformation will accordingly introduces a strong perturbation term (see equation (4.11)) that needs to be delicately controlled. To address this issue, an inductive argument is used to deal with the difficulties that arise in the structural transformation. Specifically, we need to establish error bounds of $\Vert \nabla e^{i-1}- \nabla e^{i-2}\Vert$ and $\Vert e^{i-1}- e^{i-2}\Vert _{H^{-1}}$ for $i\le n$. By integrating them into the perturbation terms, and utilizing the super-convergence characteristics of resulting nonlinear terms, we can complete the mathematical induction method. Therefore, by reconverting the scheme format using the new variable back to a form compatible with the original format of the Cahn-Hilliard equation, and using mathematical induction, the superconvergence characteristics of some nonlinear terms, as well as the spectrum argument, we derive optimal error estimates that depend only on polynomial degree of $\varepsilon ^{-1}$. To the best of our knowledge, our numerical analysis has the following properties that were not possessed in the existing literature:

(1) We have addressed the challenge of spectral estimates being ineffective in the SAV approach for the Cahn-Hilliard equation;

(2) This is the first error estimate with polynomial dependence on $\varepsilon ^{-1}$ for the IEQ/SAV-type schemes of the Cahn-Hilliard equation (1.1).

The techniques developed here can also be applied for other IEQ/SAV-type schemes of phase-field models.

The rest of this paper is organized as follows. In Sect. 2, we present the SAV reformulation of the Cahn-Hilliard equation and introduce an equivalent transformation of the SAV time-stepping scheme. In Sect. 3, we show the properties of energy decay and derive the consistency estimates for the proposed method. In Sect. 4, we present an error estimate of the semi-discrete SAV scheme to derive a convergence rate that does not depend on $\varepsilon ^{-1}$ exponentially. The spectrum estimate plays a crucial role in the proof. Finally, in Sect. 5, we present a few numerical experiments to validate the theoretical results.

2 Formulation of the Semi-Discrete SAV Scheme

In this section, we construct a backward Euler implicit-explicit type temporal semi-discrete numerical scheme based on the SAV reformulation of the CH equation (1.1), and also present an equivalent formulation of the SAV scheme.

2.1 Function Spaces

Let $W^{s, p}({\varOmega })$ denote the usual Sobolev spaces, and $H^s({\varOmega })$ denote the Hilbert spaces $W^{s, 2}({\varOmega })$ with norm $\Vert \cdot \Vert _{H^s}$. Let $\Vert \cdot \Vert$ and $(\cdot , \cdot )$ represent the $L^2$ norm and $L^2$ inner product, respectively. In addition, define for $p \ge 0$

$$\begin{aligned} H^{-p}({\varOmega }):=(H^p({\varOmega }))^*, \quad H_0^{-p}({\varOmega }):= \{ u \in H^{-p}({\varOmega }) \,| \,\langle u, 1\rangle _p = 0 \}, \end{aligned}$$

(2.1)

where $\langle \cdot , \cdot \rangle _p$ stands for the dual product between $H^p({\varOmega })$ and $H^{-p}({\varOmega })$. We denote $L_0^2({\varOmega }):= H_0^0({\varOmega })$. For $v \in L_0^2({\varOmega })$, let $- \Delta ^{-1} v:=v_1 \in H^1({\varOmega }) \cap L_0^2({\varOmega })$, where $v_1$ is the solution to

$$\begin{aligned} -\Delta v_1 = v \, \, \text{ in } \, \, {\varOmega }, \qquad \partial _{\varvec{n}} v_1 = 0 \,\, \text{ on } \,\, \partial {\varOmega }, \end{aligned}$$

(2.2)

and $\Vert v\Vert _{-1}:= \sqrt{(v, - \Delta ^{-1} v)}$.

For $v \in L_0^2({\varOmega }) \cap H^1({\varOmega })$, we have the following inequality

$$\begin{aligned} \Vert v\Vert ^2 = (\nabla v, \nabla (- \Delta )^{-1} v) \le \Vert \nabla v\Vert \Vert v\Vert _{-1}. \end{aligned}$$

(2.3)

We denote by C generic constant and $C_i$, $\tilde{C}_i$, $\tilde{C}$, $\kappa _i'$ and $\kappa _i$ specific constants, which are independent of $\tau$, h and $\varepsilon$, but may possibly depend on the domain ${\varOmega }$, T and the constants of Sobolev inequalities. We use notation $\lesssim$ in the sense that $f \lesssim g$ means that $f \le Cg$ with positive constant C independent of $\tau$, h and $\varepsilon$.

2.2 The SAV Reformulation

The SAV formulation of the CH equation (cf. [30, 31]) introduces a scalar auxiliary variable

$$\begin{aligned} r = \sqrt{\int _{\varOmega }F(u) {\mathrm d}x + c_0} \quad \text{ with }\quad g(u)=\frac{f(u)}{ \sqrt{\int _{\varOmega }F(u) {\mathrm d}x +c_0}}, \end{aligned}$$

(2.4)

with a positive $c_0$ (which guarantees that the function r has a positive lower bound), and reformulate (1.1) as

$$\begin{aligned} \partial _t u&= \Delta w & \qquad \text{ in }\,\,\,{\varOmega }\times (0,T], \end{aligned}$$

(2.5a)

$$\begin{aligned} w&= - \varepsilon \Delta u + \frac{1}{\varepsilon } r g(u) & \qquad \text{ in }\,\,\,{\varOmega }\times (0,T], \end{aligned}$$

(2.5b)

$$\begin{aligned} \frac{{\mathrm d}r}{{\mathrm d}t}&= \frac{1}{2} \big (g(u), \partial _tu\big ) & \qquad \text{ in }\,\,\,{\varOmega }\times (0,T], \end{aligned}$$

(2.5c)

$$\begin{aligned} \partial _{\varvec{n}} u&= \partial _{\varvec{n}} w = 0 & \qquad \text{ on }\,\,\, \partial {\varOmega }\times (0,T], \end{aligned}$$

(2.5d)

$$\begin{aligned} u(\cdot , 0)&= u_0 & \qquad \text{ in }\,\,\,{\varOmega }, \end{aligned}$$

(2.5e)

$$\begin{aligned} r(0)&= \sqrt{\int _{\varOmega }F(u_0) {\mathrm d}x+c_0}. & \end{aligned}$$

(2.5f)

We define an energy functional with respect to u and r:

$$\begin{aligned} E(u, r) = \frac{\varepsilon }{2} \Vert \nabla u\Vert ^2 + \frac{1}{\varepsilon } r^2, \end{aligned}$$

(2.6)

and taking the $L^2$ inner product of the first equation (2.5a) with w, of the second equation (2.5b) with $\partial _t u$, and of the third equation (2.5c) with $\frac{1}{\varepsilon } 2r$, performing integration by parts and summing up the two obtained equations, we find that

$$\begin{aligned} \frac{{\mathrm d}}{{\mathrm d}t} E(u, r) = - \Vert \nabla w\Vert ^2 = - \Vert \partial _t u\Vert _{-1}^2 \le 0. \end{aligned}$$

(2.7)

2.3 The Equivalent Formulation of the SAV Scheme

Let $\{t_n\}_{n=0}^{N+1}$ be a uniform partition of [0, T] with the time step size $\tau$, where N is a positive integer and hence $\tau =\frac{T}{N+1}$. We consider the following temporal semi-discrete SAV scheme for solving the system (2.5):

$$\begin{aligned} \left\{ \begin{aligned} \frac{u^{n+1} - u^n}{\tau }&= \Delta w^{n + 1}, \\ w^{n+1}&= - \varepsilon \Delta u^{n+1} + \frac{1}{\varepsilon } r^{n+1} g(u^n), \\ r^{n+1} - r^n&= \frac{1}{2} \big (g(u^n), u^{n+1} - u^n \big ), \\ \partial _{\varvec{n}} u^{n+1} |_{\partial {\varOmega }}&= \partial _{\varvec{n}} w^{n+1} |_{\partial {\varOmega }} = 0, \end{aligned} \right. \end{aligned}$$

(2.8)

with $u^0 = u_0$ and $r^0 = r(0)$ for $n = 0, 1, \dots , N$.

By taking the $L^2$ inner product of the first equation in (2.8) with $v = 1$, we have the following conservation property, which is important to the error estimates.

Lemma 2.1

The numerical solution of (2.8) satisfies

$$\begin{aligned} \frac{1}{|{\varOmega }|} \int _{{\varOmega }} u^n(x) \, {\mathrm d}x = \frac{1}{|{\varOmega }|} \int _{{\varOmega }} u^0(x) \, {\mathrm d}x, \qquad n = 1, \dots , N, \end{aligned}$$

(2.9)

and the error function $e^n:= u(t_n) - u^n$ satisfies

$$\begin{aligned} \int _{{\varOmega }} e^n \, {\mathrm d}x = 0, \qquad \qquad n = 1, \dots , N. \end{aligned}$$

(2.10)

Due to (2.9–2.10), both $u^{n+1} - u^n$ and $e^n$ belong to $L_0^2({\varOmega })$, allowing us to define their $\Vert \cdot \Vert _{-1}$ norm. Let $\hat{w}^{n+1} = w^{n+1} - \int _{{\varOmega }} w^{n+1} {\mathrm d}x$ such that $\hat{w}^{n+1} \in L_0^2({\varOmega })$. The first equation in (2.8) satisfies

$$\begin{aligned} \frac{u^{n+1} - u^n}{\tau }&= \Delta w^{n + 1} = \Delta \hat{w}^{n + 1}. \end{aligned}$$

(2.11)

Since both sides of the equation belong to $L_0^2({\varOmega })$, we can apply $\Delta ^{-1}$ to get

$$\begin{aligned} \Delta ^{-1} \frac{u^{n+1} - u^n}{\tau } = \hat{w}^{n + 1}, \end{aligned}$$

which combines with

$$\begin{aligned} (\hat{w}^{n + 1}, v) = (w^{n + 1}, v) \qquad \forall \,\, v \in L_0^2({\varOmega }), \end{aligned}$$

(2.12)

gives

$$\begin{aligned} \Big (\Delta ^{-1} \frac{u^{n+1} - u^n}{\tau }, v \Big ) = (w^{n + 1}, v)\qquad \forall \,\, v \in L_0^2({\varOmega }). \end{aligned}$$

(2.13)

Thus, by testing the second equation in (2.8) with $v \in L_0^2({\varOmega })$ and using (2.13), the semi-discrete SAV scheme (2.8) can be written as

$$\begin{aligned} \left\{ \begin{aligned} \Big (\Delta ^{-1} \frac{u^{n+1} - u^n}{\tau }, v\Big )&= \varepsilon \big ( \nabla u^{n+1}, \nabla v \big ) + \frac{1}{\varepsilon } r^{n+1} \big (g(u^n), v\big ) \qquad \forall \,\, v \in L_0^2({\varOmega }), \\ r^{n+1} - r^n&= \frac{1}{2} \big (g(u^n), u^{n+1} - u^n \big ), \qquad n = 0, 1, \dots , N. \end{aligned} \right. \end{aligned}$$

(2.14)

In order to avoid the exponentially dependence of the error bound on $\frac{1}{\varepsilon }$ induced by using the Gronwall inequality, we need to use a spectral estimate of the linearized Cahn-Hilliard operator, which is given in [3, 13, 25] and will be described in Sect. 4.

However, compared with the previous work [25], the SAV method (2.8) alters the structure of the CH equation such that the spectral argument can not be applied directly. To achieve the ideal error bound, we need to transform the structure into a form compatible with the CH equation (1.1) so that the spectral estimate of the linearized Cahn-Hilliard operator can be used. To this end, we define $A(v) = \sqrt{\int _{{\varOmega }} F(v) \, {\mathrm d}x + c_0}$, the Gateaux derivatives of A(v) can be defined as follows:

$$\begin{aligned} D A(v, w)&: = \frac{1}{2} \Big (\frac{f(v)}{\sqrt{\int _{{\varOmega }} F(v) \, {\mathrm d}x + c_0}}, w \Big ) = \frac{1}{2} \big ( g(v), w\big ),\end{aligned}$$

(2.15)

$$\begin{aligned} D^2 A(v, w)&: = \frac{1}{2} \frac{\int _{{\varOmega }} f'(v) w^2 \, {\mathrm d}x}{\sqrt{\int _{{\varOmega }} F(v) \, {\mathrm d}x + c_0}} - \frac{1}{4} \frac{ \big (\int _{{\varOmega }} f(v) w\, {\mathrm d}x \big )^2}{\Big (\int _{{\varOmega }} F(v) \, {\mathrm d}x + c_0\Big )^{\frac{3}{2}}}. \end{aligned}$$

(2.16)

By Taylor expansion, we derive

$$\begin{aligned} A(u^i) = A(u^{i-1}) + \frac{1}{2} \big ( g(u^{i - 1}), u^i - u^{i-1}\big ) + \frac{1}{2} D^2A(\xi _i\,; u^i - u^{i-1}), \end{aligned}$$

(2.17)

where $\xi _i = \theta u^i + (1 - \theta ) u^{i-1}$ with $\theta \in (0, 1)$. Thus we get

$$\begin{aligned} \frac{1}{2} \big (g(u^{i-1}), u^i - u^{i - 1}\big ) = A(u^i) - A(u^{i-1}) - \frac{1}{2} D^2 A (\xi _i\,; u^i - u^{i-1}), \end{aligned}$$

which together with the second equation in (2.14) implies

$$\begin{aligned} r^i - r^{i - 1} = A(u^i) - A(u^{i-1}) - \frac{1}{2} D^2 A (\xi _i\,; u^i - u^{i-1}). \end{aligned}$$

After summing up the above equation from $i = 1$ to n, we derive

$$\begin{aligned} r^n - r^0 = A(u^n) - A(u^0) - \frac{1}{2} \sum _{i = 1}^n D^2 A (\xi _i\,; u^i - u^{i-1}). \end{aligned}$$

Since $r^0 = A(u^0)$, we obtain

$$\begin{aligned} r^n = A(u^n) - \frac{1}{2} \sum _{i = 1}^n D^2 A (\xi _i\,; u^i - u^{i-1}). \end{aligned}$$

(2.18)

Then the SAV scheme (2.14) can be written as

$$\begin{aligned}&\Big (\Delta ^{-1} \frac{u^{n+1} - u^n}{\tau }, v\Big ) =\, \varepsilon \big (\nabla u^{n+1}, \nabla v \big ) + \frac{1}{\varepsilon } r^n \big (g(u^n), v \big ) + \frac{1}{\varepsilon } \big (g(u^n), v\big ) (r^{n+1} - r^n) \nonumber \\&\quad =\, \varepsilon \big (\nabla u^{n+1}, \nabla v \big ) + \frac{1}{\varepsilon } r^n \big (g(u^n), v \big ) + \frac{1}{2\varepsilon } \big (g(u^n), v\big ) \big (g(u^n) , u^{n+1} - u^n\big ), \end{aligned}$$

(2.19)

for all $v \in L_0^2({\varOmega })$, which together with (2.18) gives

$$\begin{aligned} \Big (\Delta ^{-1} \frac{u^{n+1} - u^n}{\tau }, v\Big ) =\,&\varepsilon \big (\nabla u^{n+1}, \nabla v \big ) + \frac{1}{\varepsilon } \big (f(u^n), v\big ) \nonumber \\&- \frac{1}{2\varepsilon } \big (g(u^n), v\big ) \sum _{i = 1}^n D^2 A(\xi _i\,; u^i - u^{i - 1}) \nonumber \\&+ \frac{1}{2\varepsilon } \big (g(u^n), v\big ) \big (g(u^n), u^{n+1} - u^n\big ) \qquad \forall \,\, v \in L_0^2({\varOmega }). \end{aligned}$$

(2.20)

The above equation (2.20) provides an equivalent formulation of the semi-discrete SAV scheme (2.8), which will be frequently used in the subsequent error analysis.

3 Energy Decay and Consistency Analysis

In this section, we present several inequalities related to the proposed numerical method.

3.1 Assumption and Regularity

Before presenting the detailed numerical analysis, we first make some assumptions. The quartic growth of the Ginzburg-Landau energy function $F(u) = \frac{1}{4} (u^2 - 1)^2$ at infinity poses various technical difficulties for the analysis and approximation of CH equations. Although the CH equation does not satisfy the maximum principle, if the maximum norm of the initial condition $u^0$ is bounded, it has been shown in [9] that the maximum norm of the solution of the CH equation for the truncation potential F(u) with quadratic growth rate at infinity is bounded. Therefore, it has been a common practice (cf. [33]) to consider the CH equations with a truncated F(u).

Assumption 3.1

We assume that the potential function F(u) whose derivative $f(u) = F'(u)$ satisfies the following condition:

• (1) $F \in C^4(\mathbb {R})$, $F(\pm 1) = 0$, and $F >0$ elsewhere.

• (2) $f(\pm 1) = 0$, $f'(\pm 1) > 0$, and there exists a non-negative constant L such that

  $$\begin{aligned} \max _{v \in \mathbb {R}}|f(v)| \le L, \quad \max _{v \in \mathbb {R}}|f'(v)| \le L \quad \text{ and } \quad \max _{v \in \mathbb {R}}|f''(v)| \le L. \end{aligned}$$

  (3.1)

In order to trace the dependence of the solution on the small parameter $\varepsilon > 0$, we assume that the solution of (1.1) satisfies the following conditions:

Assumption 3.2

Suppose there exist positive $\varepsilon$-independent constants $m_0$ and $\rho _j$ for $j = 1, 2, 3$ such that the solution of (1.1) satisfies

$$\begin{aligned}&\frac{1}{|{\varOmega }|} \int _{{\varOmega }} u\, {\mathrm d}x = m_0 \in (-1, 1), \end{aligned}$$

(3.2a)

$$\begin{aligned}&\mathop {\mathrm {ess\,sup}}\limits _{t \in [0, \infty ]} \Big \{ \frac{\varepsilon }{2} \Vert \nabla u\Vert ^2 + \frac{1}{\varepsilon } \int _{{\varOmega }} F(u) \, {\mathrm d}x \Big \} + \int _0^\infty \Vert u_t\Vert _{-1}^2 \, {\mathrm d}s \lesssim \varepsilon ^{-\rho _1}, \end{aligned}$$

(3.2b)

$$\begin{aligned}&\int _0^\infty \Vert u_t\Vert ^2 \, {\mathrm d}s \lesssim \varepsilon ^{- \rho _2}, \end{aligned}$$

(3.2c)

$$\begin{aligned}&\int _0^\infty \Vert \Delta ^{-1} u_{tt}\Vert _{-1}^2 \, {\mathrm d}t \lesssim \varepsilon ^{-\rho _3}. \end{aligned}$$

(3.2d)

Remark 3.1

(a) Note that the conditions (1) and (2) are satisfied by restricting the growth of F(v) for $|v| \ge M$. More precisely, for a given $M \ge 1$, we can replace $F(v) = \frac{1}{4} (v^2 - 1)^2$ by a cut-off function $\hat{F}(v) \in C^4(\mathbb {R})$ as follows:

$$\begin{aligned} \hat{F}(v) = \left\{ \begin{aligned}&((2M)^2 - 1)2M(v - 2M) + \frac{1}{4} ((2M)^2 - 1)^2 & \text{ for }\,\, v > 2M, \\&\Phi _+(v) & \text{ for }\,\, v \in (M, 2M], \\&\frac{1}{4} (v^2 - 1)^2 & \text{ for }\, \,v \in [-M, M], \\&\Phi _-(v) & \text{ for }\,\, v \in [-2M, -M), \\&-((2M)^2 - 1)2M(v + 2M) + \frac{1}{4} ((2M)^2 - 1)^2 & \text{ for }\, \,v < -2M, \end{aligned} \right. \end{aligned}$$

(3.3)

where $\Phi _+(v)$ and $\Phi _-(v) >0$ elsewhere between $M<|v|< 2M$ and satisfy the required conditions at $|v| = M$ and $|v| = 2M$, respectively. Then we replace $f(v) = (v^2 - 1)v$ by $\hat{F}'(v)$ which is

$$\begin{aligned} \hat{f}(v) = \hat{F}'(v) = \left\{ \begin{aligned}&((2M)^2 - 1)2M & \text{ for }\,\, v > 2M, \\&\Phi '_+(v) & \text{ for }\,\, v \in (M, 2M], \\&(v^2 - 1)v & \text{ for }\, \,v \in [-M, M], \\&\Phi '_-(v) & \text{ for }\,\, v \in [-2M, -M), \\&-((2M)^2 - 1)2M & \text{ for }\, \,v < -2M. \end{aligned} \right. \end{aligned}$$

(3.4)

In simplicity, we still denote the modified function $\hat{F}$ by F. It is then obvious that there exists L such that (3.1) are satisfied with f replaced by $\hat{f}$.

(b) The transformed SAV scheme (2.20) introduced a complicated term. With the condition (2) in the Assumption 3.1, we can get

$$\begin{aligned} D^2 A(v; w) \lesssim \big (\Vert f'(v)\Vert _{L^\infty } + \Vert f(v)\Vert ^2 \big ) \Vert w\Vert ^2 \lesssim \Vert w\Vert ^2, \end{aligned}$$

(3.5)

which will be frequently used to control the difficult term in the error analysis.

(c) Assumption 3.2 can be achieved in many cases. For example, suppose that f satisfies Assumption 3.1, $\partial {\varOmega }$ is of class $C^{2, 1}$, $u^0 \in H^3({\varOmega })$, and there exist positive $\varepsilon$-independent constants $\sigma _j$ for $j = 1, 2, 3$ such that

$$\begin{aligned}&E[u^0]=\frac{\varepsilon }{2}\Vert \nabla u^0\Vert ^2 + \frac{1}{\varepsilon } \int _{{\varOmega }} F(u^0)\, {\mathrm d}x \lesssim \varepsilon ^{-2\sigma _1},\end{aligned}$$

(3.6a)

$$\begin{aligned}&\Vert -\varepsilon \Delta u^0 + \frac{1}{\varepsilon } f(u^0)\Vert \lesssim \varepsilon ^{-2\sigma _2},\end{aligned}$$

(3.6b)

$$\begin{aligned}&\Vert -\varepsilon \Delta u^0 + \frac{1}{\varepsilon } f(u^0)\Vert _{H^1} \lesssim \varepsilon ^{-2\sigma _3}. \end{aligned}$$

(3.6c)

Then the estimates (3.2a)–(3.2d) can be derived by standard test function techniques and satisfy:

$$\begin{aligned} \rho _1 = 2\sigma _1, \quad \rho _2 = \max \{2\sigma _1 + 2, 2\sigma _2 - 1\}\quad \text{ and } \quad \rho _3 = \max \{2\sigma _1 + 4, 2\sigma _2 + 1, 2\sigma _3 - 1\}. \end{aligned}$$

We refer to [22, 25] for their detailed proof.

3.2 Energy Decay Structure

In this subsection, we prove the following energy decay property of the numerical solution, which comprise of the first theorem of this paper.

Theorem 3.1

(energy decay) The scheme (2.8) is unconditionally energy stable in the sense that

$$\begin{aligned} E(u^{n+1}, r^{n+1}) - E(u^n, r^n) \le - \frac{1}{\tau }\Vert u^{n+1} - u^n\Vert _{-1}^2 \le 0 \quad \text{ for } n\ge 1. \end{aligned}$$

(3.7)

Proof

Taking the inner product of the first equation in (2.8) with $-\Delta ^{-1}(u^{n+1} - u^n)$, and of the second equation with $u^{n+1} - u^n$, using (2.12) and multiplying the third equation in (2.8) by $\frac{2}{\varepsilon } r^{n+1}$, we derive that

$$\begin{aligned}&\frac{1}{\tau } \Vert u^{n+1} - u^n\Vert _{-1}^2 + \frac{\varepsilon }{2} \Big (\Vert \nabla u^{n+1}\Vert ^2 - \Vert \nabla u^n\Vert ^2 + \Vert \nabla u^{n+1} - \nabla u^n\Vert ^2 \Big )\nonumber \\&\qquad \qquad + \frac{1}{\varepsilon }r^{n+1} \big ( g(u^n), u^{n+1} - u^n \big ) = 0, \end{aligned}$$

(3.8)

$$\begin{aligned}&\frac{1}{\varepsilon } \Big ((r^{n+1})^2 - (r^n)^2 + (r^{n+1} - r^n)^2 \Big ) = \frac{1}{\varepsilon }r^{n+1} \big ( g(u^n), u^{n+1} - u^n \big ). \end{aligned}$$

(3.9)

Taking the summation of the above equations, we get

$$\begin{aligned} \frac{1}{\tau } \Vert u^{n+1} - u^n\Vert _{-1}^2 +&\frac{\varepsilon }{2} \Big (\Vert \nabla u^{n+1}\Vert ^2 - \Vert \nabla u^n\Vert ^2 + \Vert \nabla u^{n+1} - \nabla u^n\Vert ^2 \Big )\nonumber \\ +&\frac{1}{\varepsilon } \Big ((r^{n+1})^2 - (r^n)^2 + (r^{n+1} - r^n)^2 \Big ) = 0. \end{aligned}$$

(3.10)

which gives (3.7).$\square$

Remark 3.2

After summing up (3.10) from $n = 0$ to N, we get

$$\begin{aligned} \frac{\varepsilon }{2} \Vert \nabla u^{N+1}\Vert ^2 + \frac{1}{\varepsilon } (r^{N+1})^2&+ \frac{1}{\tau } \sum _{n = 0}^N\Vert u^{n+1} - u^n\Vert _{-1}^2 + \frac{\varepsilon }{2} \sum _{n = 0}^N \Vert \nabla u^{n+1} - \nabla u^n\Vert ^2 \nonumber \\&+ \frac{1}{\varepsilon }\sum _{n = 0}^N (r^{n+1} - r^n)^2 = \frac{\varepsilon }{2} \Vert \nabla u^0\Vert ^2 + \frac{1}{\varepsilon } (r^0)^2 \lesssim \varepsilon ^{-\rho _1}. \end{aligned}$$

(3.11)

It follows from $u^{n+1} - u^n \in L_0^2({\varOmega })$ and (2.3) that

$$\begin{aligned} \sum _{n = 0}^N\Vert u^{n+1} - u^n\Vert ^2 \le \,&\sum _{n = 0}^N\Vert u^{n+1} - u^n\Vert _{-1}\Vert \nabla u^{n+1} - \nabla u^n\Vert \nonumber \\ \le \,&\Big (\sum _{n = 0}^N\Vert u^{n+1} - u^n\Vert _{-1}^2 \Big )^{\frac{1}{2}} \Big (\sum _{n = 0}^N\Vert \nabla u^{n+1} - \nabla u^n\Vert ^2 \Big )^{\frac{1}{2}} \nonumber \\ \lesssim \,&\varepsilon ^{-(\rho _1 + \frac{1}{2})}\tau ^{\frac{1}{2}}. \end{aligned}$$

(3.12)

3.3 Consistency

It follows from (1.1b) that $\int _{{\varOmega }} \Delta u {\mathrm d}x = 0$ and $\int _{{\varOmega }} \partial _t u {\mathrm d}x = 0$. Since $\partial _t u \in L_0^2({\varOmega })$, we have $u(t_{n+1}) - u(t_n)\in L_0^2({\varOmega })$. For any $f \in L^2({\varOmega })$, let $\hat{f} = f - \int _{{\varOmega }} f{\mathrm d}x$ such that $\hat{f} \in L_0^2({\varOmega })$. Following the derivation of (2.11), the CH equation (1.1) can be written as

$$\begin{aligned} \Delta ^{-1} \partial _t u&= - \varepsilon \Delta u + \frac{1}{\varepsilon } \hat{f}(u). \end{aligned}$$

Similar to (2.13), we use the relation

$$\begin{aligned} (\hat{f}(u), v) = (f(u), v) \qquad \forall \,\, v \in L_0^2({\varOmega }), \end{aligned}$$

to get

$$\begin{aligned} \Big (\Delta ^{-1} \partial _t u(t_{n+1}) , v \Big )&= \varepsilon \big (\nabla u(t_{n+1}), \nabla v\big ) + \frac{1}{\varepsilon } \big (f(u(t_{n+1})), v \big ) \qquad \forall \,\, v \in L_0^2({\varOmega }). \end{aligned}$$

To derive the error estimates of the equivalent transformation (2.20) of the semi-discrete SAV scheme (2.8), we reformulate the CH equation (1.1) as the truncated form

$$\begin{aligned}&\Big (\Delta ^{-1} \frac{u(t_{n+1}) - u(t_n)}{\tau }, v \Big )\nonumber \\&= \varepsilon \big (\nabla u(t_{n+1}), \nabla v\big ) + \frac{1}{\varepsilon } \big (f(u(t_n)), v \big ) + \big (\mathcal {R}^{n+1}, v\big ) \quad \forall \,\, v \in L_0^2({\varOmega }), \end{aligned}$$

(3.13)

where the truncation error $\mathcal {R}^{n+1}$ is given by

$$\begin{aligned} \mathcal {R}^{n+1} {:}{=} \Big [ \Delta ^{-1} \frac{u(t_{n+1}) - u(t_n)}{\tau } - \Delta ^{-1} \partial _t u(t_{n+1}) \Big ] + \frac{1}{\varepsilon }\Big [f(u(t_{n+1})) - f(u(t_n))\Big ]. \end{aligned}$$

(3.14)

Lemma 3.1

(consistency estimate) Suppose that assumptions 3.1 and 3.2 hold, then we have the following consistency estimate:

$$\begin{aligned} \tau \sum _{n = 0}^N \Vert \mathcal {R}^{n+1} \Vert _{H^{-1}}^2 \le C \varepsilon ^{- \max \{\rho _2+2, \,\rho _3 \}} \tau ^2. \end{aligned}$$

(3.15)

Proof

For any $\varphi \in L_0^2({\varOmega })$, there holds $\varphi _1 = - \Delta ^{-1} \varphi \in L_0^2({\varOmega }) \cap H^1({\varOmega })$ with $\partial _{\varvec{n}} \varphi _1 = 0$ on $\partial {\varOmega }$ and $\Vert \varphi \Vert _{-1} = \Vert \nabla \varphi _1\Vert$, which gives

$$\begin{aligned} \Vert \varphi \Vert _{H^{-1}}&= \sup _{v \,\in H^1({\varOmega })} \frac{(\varphi , v )}{\Vert v\Vert _{H^1}} = \sup _{v \,\in H^1({\varOmega })} \frac{\big (-\Delta (-\Delta ^{-1})\varphi , v \big ) }{\Vert v\Vert _{H^1}} = \sup _{v \,\in H^1({\varOmega })} \frac{\big (\nabla \varphi _1, \nabla v \big ) }{\Vert v\Vert _{H^1}} \nonumber \\&\le \Vert \nabla \varphi _1\Vert = \Vert \varphi \Vert _{-1}. \end{aligned}$$

(3.16)

Since $\partial _t u \in L_0^2({\varOmega })$, it follows that $u_{tt} \in L_0^2({\varOmega })$. By the definition of $\Delta ^{-1}: L_0^2({\varOmega }) \rightarrow H^1({\varOmega }) \cap L_0^2({\varOmega })$, we have $\Delta ^{-1} u_{tt} \in L_0^2({\varOmega })$. By performing standard calculations, we obtain

$$\begin{aligned} \Big \Vert \Delta ^{-1} \frac{u(t_{n+1}) - u(t_n)}{\tau } - \Delta ^{-1} \partial _t u(t_{n+1}) \Big \Vert _{H^{-1}}^2 = \,&\Big \Vert \frac{1}{\tau }\int _{t_n}^{t_{n+1}} (s - t_n) \Delta ^{-1} u_{tt} \, {\mathrm d}s \Big \Vert _{H^{-1}}^2 \nonumber \\ \lesssim \,&\tau \int _{t_n}^{t_{n+1}} \Vert \Delta ^{-1} u_{tt} \Vert _{H^{-1}}^2 \, {\mathrm d}s \nonumber \\ \le \,&\tau \int _{t_n}^{t_{n+1}} \Vert \Delta ^{-1} u_{tt} \Vert _{-1}^2 \, {\mathrm d}s, \end{aligned}$$

(3.17)

and

$$\begin{aligned} \Big \Vert \frac{1}{\varepsilon } f(u(t_{n+1})) - \frac{1}{\varepsilon } f(u(t_n)) \Big \Vert _{H^{-1}}^2 = \,&\varepsilon ^{-2} \sup _{v \, \in H^1({\varOmega })} \frac{\big (f(u(t_{n+1})) - f(u(t_n)), v \big )^2 }{\Vert v\Vert _{H^1}^2} \nonumber \\ \le \,&\varepsilon ^{-2} \sup _{v \,\in H^1({\varOmega })} \frac{ \Vert f'(\xi ^n)\Vert _{L^3}^2 \Vert u(t_{n+1}) - u(t_n)\Vert ^2 \Vert v\Vert _{L^6}^2 }{\Vert v\Vert _{H^1}^2} \nonumber \\ \lesssim \,&\varepsilon ^{-2}\Vert u(t_{n+1}) - u(t_n)\Vert ^2 \nonumber \\ \lesssim \,&\varepsilon ^{-2} \tau \int _{t_n}^{t_{n+1}} \Vert u_t\Vert ^2\, {\mathrm d}s, \end{aligned}$$

(3.18)

where $\xi ^n$ is between $u(t_n)$ and $u(t_{n+1})$. Thus, we have

$$\begin{aligned} \tau \sum _{n = 0}^N \Vert \mathcal {R}^{n+1} \Vert _{H^{-1}}^2 \lesssim \,&\tau ^2 \int _{0}^{T} \Vert \Delta ^{-1} u_{tt} \Vert _{-1}^2 \, {\mathrm d}s + \varepsilon ^{-2} \tau ^2 \int _{0}^{T} \Vert u_t\Vert ^2\, {\mathrm d}s \nonumber \\ \lesssim \,&\varepsilon ^{- \rho _3}\tau ^2 + \varepsilon ^{- (\rho _2 + 2)} \tau ^2. \end{aligned}$$

(3.19)

The proof is completed. $\square$

4 Error Estimates

In this section, we will derive the error bound of the semi-discrete scheme (2.8), in which the focus is to obtain the polynomial type dependence of the error bound on $\varepsilon ^{-1}$. If we use the usual error estimate of the SAV numerical scheme (2.8), the error growth depends on $\varepsilon ^{-1}$ exponentially. To avoid the exponential dependence on $\varepsilon ^{-1}$ induced by using the Gronwall inequality, we need to use a spectral estimate of the linearized Cahn-Hilliard operator, which is given in [3, 13, 25].

Lemma 4.1

(spectral estimate) Suppose that Assumption 3.1 holds. Then there exist $0< \varepsilon _0<< 1$ and a positive constant $\lambda _0$ such that the principle eigenvalue of the linearized Cahn-Hilliard operator

$$\begin{aligned} \mathcal {L}_{CH} := \Delta (\varepsilon \Delta - \frac{1}{\varepsilon } f'(u) I) \end{aligned}$$

satisfies for all $t \in [0, T]$

$$\begin{aligned} \lambda _{CH} = \inf _{\begin{array}{c} 0 \ne v \in H^1({\varOmega }) \\ \Delta w = v \end{array}} \frac{\varepsilon \Vert \nabla v\Vert ^2 + \frac{1}{\varepsilon }\big (f'(u(\cdot , t)v, v)\big )}{\Vert \nabla w\Vert ^2} \ge -\lambda _0, \end{aligned}$$

(4.1)

for $\varepsilon \in (0, \varepsilon _0)$, where I denotes the identity operator and u is the solution of the Cahn-Hilliard problem (1.1).

We will now prove the following error estimates for the semi-discrete numerical scheme, which is the main result of this paper.

Theorem 4.1

(error estimate) We assume that assumptions 3.1 and 3.2 hold and that

$$\begin{aligned} \tau \le \tilde{C} \varepsilon ^{\beta _0} \quad \text{ with } \quad \beta _0 = \frac{4\alpha _0 + 32 + 4d}{4 - d}, \end{aligned}$$

(4.2)

then the discrete solution given by (2.8) satisfies the following error estimate for $e^n = u(t_n) - u^n$:

$$\begin{aligned} \max _{1 \le n \le m}\Vert e^n\Vert _{-1}^2&+ \frac{1}{2} \sum _{n = 1}^m \Vert e^n - e^{n-1}\Vert _{-1}^2 + \frac{1}{2} \varepsilon ^4\sum _{n = 1}^m \tau \Vert \nabla e^n\Vert ^2 + \max _{1 \le n \le m} \tau \varepsilon ^4 \Vert \nabla e^m\Vert ^2 \nonumber \\&+ \frac{1}{2} \varepsilon ^4 \sum _{n = 1}^m \tau \Vert \nabla e^n - \nabla e^{n-1}\Vert ^2 \le \kappa _0\varepsilon ^{- \alpha _0} \tau ^2. \end{aligned}$$

(4.3)

where $\alpha _0:= \max \{ \rho _1 + 3, \, 2\rho _2 + 4, \, \rho _2 + 6, \, \rho _3 + 4\}$ and the constant $\kappa _0$ is independent of $\tau$ when $\tau$ is sufficiently small. The specific values of $\tilde{C}$ and $\kappa _0$ will be given in the proof.

Proof

We use the mathematical induction as follows. The proof is split into four steps. The first step gives the error estimate for the first step $t = t^1$. Steps two and three use the spectral estimate (4.1) to avoid exponential blow-up in $\varepsilon ^{-1}$ of the error constants. In the last step, an inductive argument is used to conclude the proof.

Step 1: Estimation of $\Vert e^1\Vert _{-1}^2 + \tau \varepsilon \Vert \nabla e^1\Vert ^2$. For $n = 0$ in (2.14), we have

$$\begin{aligned} \left\{ \begin{aligned} \Big (\Delta ^{-1} \frac{u^1 - u^0}{\tau }, v\Big )&= \varepsilon \big (\nabla u^1, \nabla v \big ) + \frac{1}{\varepsilon } r^1 \big (g(u^0), v\big ) \qquad \forall \,\, v \in L_0^2({\varOmega }),\\ r^1&= r^0 + \frac{1}{2} \big (g(u^0), u^1 - u^0 \big ), \end{aligned} \right. \end{aligned}$$

(4.4)

After plugging the second equation into the first equation, we get

$$\begin{aligned} \Big (\Delta ^{-1} \frac{u^1 - u^0}{\tau }, v\Big )&= \varepsilon \big (\nabla u^1, \nabla v \big ) + \frac{1}{\varepsilon } \big (f(u^0), v\big ) + \frac{1}{2\varepsilon } \big (g(u^0), v\big )\big (g(u^0), u^1 - u^0 \big ), \end{aligned}$$

(4.5)

for all $v \in L_0^2({\varOmega })$. Subtracting (4.5) from (3.13), we get the corresponding error equation

$$\begin{aligned} \Big (\Delta ^{-1} \frac{e^1}{\tau }, v\Big )&= \varepsilon \big ( \nabla e^1, \nabla v \big ) + \big (\mathcal {R}^1, v\big ) - \frac{1}{2\varepsilon } \big (g(u^0), v\big )\big (g(u^0), u^1 - u^0 \big ) \qquad \forall \,\, v \in L_0^2({\varOmega }). \end{aligned}$$

(4.6)

Testing the above equation with $v = e^1 \in L_0^2({\varOmega })$, we have

$$\begin{aligned} \Vert e^1\Vert _{-1}^2 + \tau \varepsilon \Vert \nabla e^1\Vert ^2 = - \tau (\mathcal {R}^1, e^1) + \frac{\tau }{2\varepsilon } \big (g(u^0), u^1 - u^0 \big )\big (g(u^0), e^1 \big ). \end{aligned}$$

(4.7)

Using Poincaré’s inequality for $e^1 \in L_0^2({\varOmega })$ and Lemma 3.1, we have

$$\begin{aligned} \tau (\mathcal {R}^1, e^1)&\le \frac{\varepsilon \tau }{4}\Vert \nabla e^1\Vert ^2 + C\varepsilon ^{-1} \tau \Vert \mathcal {R}^1\Vert _{H^{-1}}^2 \nonumber \\&\le \frac{\varepsilon \tau }{4} \Vert \nabla e^1\Vert ^2 + C\varepsilon ^{- \max \{\rho _2+3, \,\rho _3 +1\}}\tau ^2. \end{aligned}$$

(4.8)

From (3.12), we have

$$\begin{aligned} \frac{\tau }{2\varepsilon } \big (g(u^0), u^1 - u^0\big ) \big (g(u^0), e^1\big ) \le \,&C \tau \varepsilon ^{-1} \Vert g(u^0)\Vert ^2\Vert u^1 - u^0\Vert \Vert e^1\Vert \nonumber \\ \le \,&C \tau \varepsilon ^{-1} \Vert u^1 - u^0\Vert \Vert e^1\Vert _{-1}^{\frac{1}{2}}\Vert \nabla e^1\Vert ^{\frac{1}{2}} \nonumber \\ \le \,&\frac{1}{2} \tau ^{\frac{1}{2}} \varepsilon ^{\frac{1}{2}}\Vert e^1\Vert _{-1} \Vert \nabla e^1\Vert + C \tau ^{\frac{3}{2}} \varepsilon ^{-\frac{5}{2}} \Vert u^1 - u^0\Vert ^2 \nonumber \\ \le \,&\frac{1}{4} \Vert e^1\Vert _{-1}^2 + \frac{\varepsilon \tau }{4} \Vert \nabla e^1\Vert ^2 + C\varepsilon ^{-(\rho _1 + 3)}\tau ^2 \end{aligned}$$

(4.9)

Thus, there exists a positive constant $\kappa _0'$ such that

$$\begin{aligned} \Vert e^1\Vert _{-1}^2 + \tau \varepsilon \Vert \nabla e^1\Vert ^2 \le \kappa _0' \varepsilon ^{- \max \{\rho _1 + 3, \,\rho _2+3, \,\rho _3 +1\}} \tau ^2. \end{aligned}$$

(4.10)

We use the mathematical induction as follows. For $m = 1$, (4.3) holds from (4.10). We suppose (4.3) holds for $m = 1, 2, \dots , N$, and show that (4.3) is also valid for $m = N+1$.

Subtracting (2.20) from (3.13), by denoting $e^n = u(t_n) - u^n$, we get the error equation

$$\begin{aligned} \Big (\Delta ^{-1} \frac{e^{n+1} - e^n}{\tau }, v\Big ) =\,&\varepsilon \big (\nabla e^{n+1}, \nabla v\big ) + \frac{1}{\varepsilon } \Big (f(u(t_n)) - f(u^n), v\Big ) \nonumber \\&+ \frac{1}{2\varepsilon } \big (g(u^n), v\big ) \sum _{i = 1}^n D^2 A(\xi _i\,; u^i - u^{i - 1}) \nonumber \\&- \frac{1}{2\varepsilon } \big (g(u^n), v\big ) \big (g(u^n), u^{n+1} - u^n\big ) + \big (\mathcal {R}^{n+1}, v\big ) \quad \forall \,\, v \in L_0^2({\varOmega }). \end{aligned}$$

(4.11)

Step 2: Estimation of $\displaystyle \Vert \nabla e^{N+1}\Vert ^2 + \sum _{ n = 0}^N \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2$. By testing (4.11) with $v = e^{n+1} - e^n \in L_0^2({\varOmega })$, we have

$$\begin{aligned} \frac{1}{2}\Vert \nabla e^{n+1}\Vert ^2 - \,&\frac{1}{2}\Vert \nabla e^n\Vert ^2 + \frac{1}{2}\Vert \nabla e^{n+1} - \nabla e^n\Vert ^2 + \frac{\tau }{\varepsilon } \Big \Vert \frac{e^{n+1} - e^n}{\tau } \Big \Vert _{-1}^2 \nonumber \\ =\,&- \frac{1}{\varepsilon ^2} \Big (f(u(t_n)) - f(u^n), e^{n+1} - e^n\Big ) \nonumber \\&- \frac{1}{2\varepsilon ^2} \big (g(u^n), e^{n+1} - e^n\big ) \sum _{i = 1}^n D^2 A(\xi _i\,; u^i - u^{i - 1}) \nonumber \\&+ \frac{1}{2\varepsilon ^2} \big (g(u^n), e^{n+1} - e^n\big ) (g(u^n), u^{n+1} - u^n) - \frac{1}{\varepsilon }\big (\mathcal {R}^{n+1}, e^{n+1} - e^n\big ) \nonumber \\ =: \,&K_1 + K_2 + K_3 + K_4. \end{aligned}$$

(4.12)

It follows from (2.3) that

$$\begin{aligned} K_1 \le \,&\varepsilon ^{-2} \Vert f'(w^n)\Vert _{L^\infty } \Vert e^n\Vert \Vert e^{n+1} - e^n\Vert \nonumber \\ \le \,&C\varepsilon ^{-2} \Vert \nabla e^n\Vert ^{\frac{1}{2}} \Vert e^n\Vert _{-1}^{\frac{1}{2}}\Vert \nabla e^{n+1} - \nabla e^n\Vert ^{\frac{1}{2}} \Vert e^{n+1} - e^n\Vert _{-1}^{\frac{1}{2}}\nonumber \\ \le \,&C\frac{1}{\gamma _0} \varepsilon ^{-2}\tau ^{\frac{1}{2}} \Vert \nabla e^n\Vert \Vert e^n\Vert _{-1} + \gamma _0\varepsilon ^{-2} \tau ^{-\frac{1}{2}} \Vert \nabla e^{n+1} - \nabla e^n\Vert \Vert e^{n+1} - e^n\Vert _{-1}\nonumber \\ \le \,&C\tau \Vert \nabla e^n\Vert ^2 + C \varepsilon ^{-4} \Vert e^n\Vert _{-1}^2 + \frac{1}{2} \gamma _0\tau ^{-1}\varepsilon ^{-4} \Vert e^{n+1} - e^n \Vert _{-1}^2 + \frac{1}{2} \gamma _0 \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2 \end{aligned}$$

(4.13)

where $w_n = \theta u^n + (1 - \theta ) u(t_n)$ with $\theta \in (0, 1)$, and $\gamma _0$ is a sufficiently small constant such that the two terms involving $\Vert e^{n+1} - e^n\Vert _{-1}^2$ and $\Vert \nabla e^{n+1} - \nabla e^n\Vert ^2$ can be absorbed by the left hand side in the Step 4.

By using $D^2 A(v; w) \lesssim \big (\Vert f'(v)\Vert _{L^\infty } + \Vert f(v)\Vert ^2 \big ) \Vert w\Vert ^2$, the term $K_2$ can be bounded as

$$\begin{aligned} K_2 \lesssim \,&\frac{1}{2\varepsilon ^2} \big |\big (g(u^n), e^{n+1} - e^n\big )\big | \sum _{i = 1}^n \Vert u^i - u^{i - 1}\Vert ^2 \nonumber \\ \lesssim \,&\frac{1}{\varepsilon ^2} \big |\big (g(u^n), e^{n+1} - e^n\big )\big | \sum _{i = 1}^n \Vert e^i - e^{i - 1}\Vert ^2 \nonumber \\&+ \frac{1}{\varepsilon ^2} \big |\big (g(u^n), e^{n+1} - e^n\big )\big | \sum _{i = 1}^n\Vert u(t_i) - u(t_{i - 1})\Vert ^2 \nonumber \\ =:\,&K_{21} + K_{22}. \end{aligned}$$

(4.14)

From the assumptions of induction, we have

$$\begin{aligned} \big (g(u^n), e^{n+1} - e^n\big )&= \frac{1}{ \sqrt{\int _{{\varOmega }} F(u^n) \, {\mathrm d}x + c_0}} \big (f(u^n), e^{n+1} - e^n\big ) \nonumber \\&\le C\Vert f(u^n)\Vert \Vert e^{n+1} - e^n\Vert \nonumber \\&\le C\Vert \nabla e^{n+1} - \nabla e^n\Vert ^{\frac{1}{2}} \Vert e^{n+1} - e^n\Vert _{-1}^{\frac{1}{2}}, \end{aligned}$$

(4.15)

$$\begin{aligned} \sum _{i = 1}^n \Vert e^i - e^{i-1}\Vert ^2&\lesssim \Big (\sum _{i = 1}^n \Vert e^i - e^{i-1}\Vert _{-1}^2 \Big )^{\frac{1}{2}} \Big (\sum _{i = 1}^n \Vert \nabla e^i - \nabla e^{i-1}\Vert ^2 \Big )^{\frac{1}{2}}\nonumber \\&\lesssim \kappa _0 \varepsilon ^{-(\alpha _0 + 2)}\tau ^{\frac{3}{2}}, \end{aligned}$$

(4.16)

$$\begin{aligned} \sum _{i = 1}^n \Vert u(t_i) - u(t_{i - 1})\Vert ^2&\lesssim \tau \sum _{i = 1}^n \int _{t_{i-1}}^{t_i} \Vert u_t\Vert ^2 \, {\mathrm d}s \lesssim \tau \int _{0}^{t_n} \Vert u_t\Vert ^2 \, {\mathrm d}s \nonumber \\&\lesssim \varepsilon ^{-\rho _2}\tau . \end{aligned}$$

(4.17)

Using the estimates above, we have

$$\begin{aligned} K_{21} \le \,&C \varepsilon ^{-2} \Vert \nabla e^{n+1} - \nabla e^n\Vert ^{\frac{1}{2}} \Vert e^{n+1} - e^n\Vert _{-1}^{\frac{1}{2}} \kappa _0 \varepsilon ^{-(\alpha _0 + 2)}\tau ^{\frac{3}{2}} \nonumber \\ \le \,&\gamma _0 \varepsilon ^{-2} \tau ^{-\frac{1}{2}} \Vert \nabla e^{n+1} - \nabla e^n\Vert \Vert e^{n+1} - e^n\Vert _{-1} + C \frac{1}{\gamma _0}\kappa _0^2 \varepsilon ^{- (2\alpha _0 + 6)}\tau ^{\frac{7}{2}} \nonumber \\ \le \,&\frac{1}{2} \gamma _0 \tau ^{-1} \varepsilon ^{-4} \Vert e^{n+1} - e^n \Vert _{-1}^2 + \frac{1}{2} \gamma _0 \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2 + C\kappa _0^2 \varepsilon ^{- (2\alpha _0 + 6)}\tau ^{\frac{7}{2}} . \end{aligned}$$

(4.18)

and

$$\begin{aligned} K_{22} \le \,&C\varepsilon ^{-2} \Vert e^{n+1} - e^n\Vert _{-1}^{\frac{1}{2}} \Vert \nabla e^{n+1} - \nabla e^n\Vert ^{\frac{1}{2}} \varepsilon ^{-\rho _2} \tau \nonumber \\ \le \,&\gamma _0 \varepsilon ^{-2}\tau ^{-\frac{1}{2}} \Vert \nabla e^{n+1} - \nabla e^n\Vert \Vert e^{n+1} - e^n\Vert _{-1} + C \frac{1}{\gamma _0} \varepsilon ^{- (2\rho _2 + 2)} \tau ^{\frac{5}{2}} \nonumber \\ \le \,&\frac{1}{2} \gamma _0 \tau ^{-1} \varepsilon ^{-4} \Vert e^{n+1} - e^n \Vert _{-1}^2 + \frac{1}{2} \gamma _0 \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2 + C \varepsilon ^{- (2\rho _2 + 2)} \tau ^{\frac{5}{2}} . \end{aligned}$$

(4.19)

In addition, we divided the term $K_3$ into two parts as

$$\begin{aligned} K_3 = \,&\frac{1}{2\varepsilon ^2} \big (g(u^n), e^{n+1} - e^n\big ) \big (g(u^n), u^{n+1} - u^n\big ) \nonumber \\ = \,&-\frac{1}{2\varepsilon ^2} \big (g(u^n), e^{n+1} - e^n\big ) \big (g(u^n), e^{n+1} - e^n\big ) \nonumber \\&+ \frac{1}{2\varepsilon ^2} \big (g(u^n), e^{n+1} - e^n\big ) \big (g(u^n), u(t_{n+1}) - u(t_n)\big ) \nonumber \\ =:\,&K_{31} + K_{32}, \end{aligned}$$

(4.20)

where

$$\begin{aligned} K_{31} \le \,&C \varepsilon ^{-2} \Vert g(u^n)\Vert ^2 \Vert e^{n+1} - e^n\Vert _{-1} \Vert \nabla e^{n+1} - \nabla e^n\Vert \nonumber \\ \le \,&C \frac{1}{\gamma _0}\varepsilon ^{-4} \Vert e^{n+1} - e^n \Vert _{-1}^2 + \frac{1}{2} \gamma _0\Vert \nabla e^{n+1} - \nabla e^n\Vert ^2, \end{aligned}$$

(4.21)

and

$$\begin{aligned} K_{32} \le \,&C \varepsilon ^{-2} \Vert e^{n+1} - e^n\Vert _{-1}^{\frac{1}{2}}\Vert \nabla e^{n+1} - \nabla e^n\Vert ^{\frac{1}{2}} \Vert u(t_{n+1}) - u(t_n)\Vert \nonumber \\ \le \,&\gamma _0 \varepsilon ^{-2} \tau ^{-\frac{1}{2}} \Vert e^{n+1} - e^n\Vert _{-1}\Vert \nabla e^{n+1} - \nabla e^n\Vert + C \gamma _0^{-1}\varepsilon ^{- 2}\tau ^{\frac{3}{2}} \int _{t_n}^{t_{n+1}}\Vert u_t\Vert ^2 \, {\mathrm d}s \nonumber \\ \le \,&\frac{1}{2} \gamma _0\tau ^{-1} \varepsilon ^{-4} \Vert e^{n+1} - e^n \Vert _{-1}^2 + \frac{1}{2} \gamma _0 \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2 + C\varepsilon ^{- 2}\tau ^{\frac{3}{2}} \int _{t_n}^{t_{n+1}}\Vert u_t\Vert ^2 \, {\mathrm d}s , \end{aligned}$$

(4.22)

By using Poincaré’s inequality for $e^{n+1} - e^n \in L_0^2({\varOmega })$, the last term on the right hand side of (4.12) can be bounded by

$$\begin{aligned} K_4 = \frac{1}{\varepsilon }\big (\mathcal {R}^{n+1}, e^{n+1} - e^n\big ) \le \frac{1}{2} \gamma _0 \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2 + C\frac{1}{\gamma _0} \varepsilon ^{-2} \Vert \mathcal {R}^{n+1}\Vert _{H^{-1}}^2. \end{aligned}$$

(4.23)

Combining these estimates (4.13)–(4.23) together with (4.12), and taking the summation for $n = 0$ to N, we have

$$\begin{aligned} \Vert \nabla e^{N+1}\Vert ^2&+ (1 - 6 \gamma _0) \sum _{n = 0}^N \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2 \le C\tau \sum _{n = 0}^N \Vert \nabla e^n\Vert ^2 + C \varepsilon ^{-4}\sum _{n = 0}^N \Vert e^n\Vert _{-1}^2 \nonumber \\&+ (C + 4 \gamma _0 \tau ^{-1})\varepsilon ^{-4}\sum _{n = 0}^N \Vert e^{n+1} - e^n \Vert _{-1}^2 + C\varepsilon ^{-2} \sum _{n = 0}^N \Vert \mathcal {R}^{n+1}\Vert _{H^{-1}}^2 \nonumber \\&+ C\kappa _0^2 \varepsilon ^{- (2\alpha _0 + 6)}\tau ^{\frac{5}{2}} + C \varepsilon ^{- (2\rho _2 + 2)} \tau ^{\frac{3}{2}} + C \varepsilon ^{-2}\tau ^{\frac{3}{2}} \int _0^T\Vert u_t\Vert ^2 \, {\mathrm d}s. \end{aligned}$$

(4.24)

Step 3: Estimation of $\displaystyle \Vert e^{N+1}\Vert _{-1}^2 + \sum _{n = 0}^N \Vert e^n - e^{n-1}\Vert _{-1}^2$. By testing (4.11) with $v = e^{n+1} \in L_0^2({\varOmega })$, we get

$$\begin{aligned} \frac{1}{2\tau } \Big (\Vert e^{n+1}\Vert _{-1}^2&- \Vert e^n\Vert _{-1}^2 + \Vert e^{n+1} - e^n\Vert _{-1}^2 \Big ) + \varepsilon \Vert \nabla e^{n+1}\Vert ^2 + \frac{1}{\varepsilon } \Big (f(u(t_n)) - f(u^n), e^{n+1}\Big ) \nonumber \\&+ \frac{1}{2\varepsilon } \big (g(u^n), e^{n+1} \big ) \sum _{i = 1}^n D^2 A(\xi _i\,; u^i - u^{i - 1}) \nonumber \\&- \frac{1}{2\varepsilon } \big (g(u^n), e^{n+1}\big ) \big (g(u^n), u^{n+1} - u^n \big ) + \big (\mathcal {R}^{n+1}, e^{n+1} \big ) = 0. \end{aligned}$$

(4.25)

We denote $w^n = \theta u^n + (1-\theta ) u(t_n)$ with $\theta \in (0, 1)$, and using the Taylor expansion to get

$$\begin{aligned} \frac{1}{\varepsilon } \Big (f(u(t_n)) - f(u^n), e^{n+1}\Big ) =\,&\frac{1}{\varepsilon } \Big (f'(u(t_n)) e^n - \frac{1}{2} f''(w^n)(e^n)^2, e^{n+1}\Big ) \nonumber \\ =\,&\frac{1}{\varepsilon } \Big (f'(u(t_n)) e^{n+1}, e^{n+1}\Big ) - \frac{1}{\varepsilon } \Big (f'(u(t_n)) (e^{n+1} - e^n), e^{n+1}\Big ) \nonumber \\&- \frac{1}{2\varepsilon } \Big (f''(w^n)(e^n)^2, e^{n+1}\Big ). \end{aligned}$$

Then, (4.25) becomes

$$\begin{aligned} \frac{1}{2\tau } \Big (\Vert e^{n+1}\Vert _{-1}^2 -\,&\Vert e^n\Vert _{-1}^2 + \Vert e^{n+1} - e^n\Vert _{-1}^2 \Big ) + \varepsilon \Vert \nabla e^{n+1}\Vert ^2 + \frac{1}{\varepsilon } \Big (f'(u(t_n)) e^{n+1}, e^{n+1}\Big ) \nonumber \\ = \,&\frac{1}{\varepsilon } \Big (f'(u(t_n)) (e^{n+1} - e^n), e^{n+1}\Big ) + \frac{1}{2\varepsilon } \Big (f''(w^n)(e^n)^2, e^{n+1}\Big ) \nonumber \\&- \frac{1}{2\varepsilon } \big (g(u^n), e^{n+1} \big ) \sum _{i = 1}^n D^2 A(\xi _i\,; u^i - u^{i - 1}) \nonumber \\&+ \frac{1}{2\varepsilon } \big (g(u^n), e^{n+1}\big ) \big (g(u^n), u^{n+1} - u^n \big ) - \big (\mathcal {R}^{n+1}, e^{n+1} \big ) \nonumber \\ =: \,&T_1 + T_2 + T_3 + T_4 + T_5. \end{aligned}$$

(4.26)

Using Lemma 4.1, there holds

$$\begin{aligned} \varepsilon \Vert \nabla e^{n+1}\Vert ^2 + \frac{1}{\varepsilon } \Big (f'(u(t_n)) e^{n+1}, e^{n+1}\Big ) \ge - \lambda _0\Vert e^{n+1}\Vert _{-1}^2. \end{aligned}$$

(4.27)

If the entire $\varepsilon \Vert \nabla e^{n+1}\Vert ^2$ term is used to control the the term $\varepsilon ^{-1}\big (f'(u(t_n)) e^{n+1}, e^{n+1}\big )$, we will not be able to control the $\Vert \nabla e^{n+1}\Vert ^2$ terms in $T_j$, $j = 1, \dots , 5$. So we apply (4.27) with a scaling factor $(1 - \eta )$ close to but smaller than 1, to get

$$\begin{aligned} - (1 - \eta )\frac{1}{\varepsilon } \Big (f'(u(t_n)) e^{n+1}, e^{n+1}\Big ) \le (1 - \eta )\lambda _0 \Vert e^{n+1}\Vert _{-1}^2 + (1 - \eta ) \varepsilon \Vert \nabla e^{n+1}\Vert ^2. \end{aligned}$$

(4.28)

On the other hand,

$$\begin{aligned} - \frac{\eta }{\varepsilon } \Big (f'(u(t_n)) e^{n+1}, e^{n+1}\Big ) \le \frac{C\eta }{\varepsilon } \Vert e^{n+1}\Vert ^2 \le \frac{C\eta }{\varepsilon ^2 \eta _1} \Vert e^{n+1}\Vert _{-1}^2 + \frac{\eta \eta _1}{4}\Vert \nabla e^{n+1}\Vert ^2. \end{aligned}$$

(4.29)

The first term $T_1$ on the right-hand side of (4.26) can be bounded by

$$\begin{aligned} T_1 \le \,&\varepsilon ^{-1} \Vert f'(u(t_n))\Vert _{L^\infty } \Vert e^{n+1} - e^n\Vert \Vert e^{n+1}\Vert \nonumber \\ \le \,&C\varepsilon ^{-1} \Vert e^{n+1} - e^n\Vert _{-1}^{\frac{1}{2}} \Vert \nabla e^{n+1} - \nabla e^n\Vert ^{\frac{1}{2}} \Vert e^{n+1}\Vert _{-1}^{\frac{1}{2}} \Vert \nabla e^{n+1}\Vert ^{\frac{1}{2}} \nonumber \\ \le \,&\tau ^{-\frac{1}{2}} \Vert e^{n+1} - e^n\Vert _{-1} \Vert e^{n+1}\Vert _{-1} + C \varepsilon ^{- 2}\tau ^{\frac{1}{2}} \Vert \nabla e^{n+1} - \nabla e^n\Vert \Vert \nabla e^{n+1}\Vert \nonumber \\ \le \,&\frac{1}{2} \gamma _0 \tau ^{-1} \Vert e^{n+1} - e^n\Vert _{-1}^2 + C\frac{1}{\gamma _0}\Vert e^{n+1}\Vert _{-1}^2 + \frac{1}{16} \varepsilon ^{\eta _2} \Vert \nabla e^{n+1}\Vert ^2 + T_1^*. \end{aligned}$$

(4.30)

where $\displaystyle T_1^*: = C\varepsilon ^{-(\eta _2 + 4)}\tau \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2$ and $\gamma _0$ is sufficiently small. To control the last term $T_1^*$ on the right-hand side of (4.30), we assume that $\tau \le \tilde{C}_1 \varepsilon ^{\eta _2 + 8}$ to get

$$\begin{aligned} \tau \sum _{n = 0}^N T_1^* := \,&C\varepsilon ^{-(\eta _2 + 4)}\tau ^2 \sum _{n = 0}^N \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2 \nonumber \\ = \,&C\varepsilon ^{-(\eta _2 + 8)}\tau \Big (\varepsilon ^4\tau \sum _{n = 0}^N \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2 \Big )\nonumber \\ \le \,&\frac{1}{16} \varepsilon ^4 \tau \sum _{n = 0}^N \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2. \end{aligned}$$

Using the Sobolev interpolation inequality, we have for $v \in L_0^2({\varOmega }) \cap H^1({\varOmega })$

$$\begin{aligned} \Vert v\Vert _{L^4} \lesssim \Vert \nabla v\Vert ^{\frac{d}{4}}\Vert v\Vert ^{1- \frac{d}{4}} \lesssim \Vert \nabla v\Vert ^{\frac{d}{4}} \Vert v\Vert _{-1}^{\frac{1}{2} - \frac{d}{8}}\Vert \nabla v\Vert ^{\frac{1}{2} - \frac{d}{8} }= \Vert \nabla v\Vert ^{\frac{1}{2} + \frac{d}{8}} \Vert v\Vert _{-1}^{\frac{1}{2} - \frac{d}{8}}, \end{aligned}$$

(4.31)

which together with the assumptions of induction yields

$$\begin{aligned} T_2 \lesssim \,&\varepsilon ^{-1} \Vert e^n\Vert _{L^4}^2\Vert e^{n+1}\Vert \nonumber \\ \le \,&C \varepsilon ^{-1} \Vert \nabla e^n\Vert ^{1 + \frac{d}{4}} \Vert e^n\Vert _{-1}^{ 1 - \frac{d}{4}} \Vert e^{n+1}\Vert _{-1}^{\frac{1}{2}} \Vert \nabla e^{n+1}\Vert ^{\frac{1}{2}} \nonumber \\ \le \,&C \kappa _0 \varepsilon ^{- ( \alpha _0+ 3 + \frac{d}{2})} \tau ^{\frac{3}{2} - \frac{d}{8}} \Vert e^{n+1}\Vert _{-1}^{\frac{1}{2}} \Vert \nabla e^{n+1}\Vert ^{\frac{1}{2}} \nonumber \\ \le \,&C \varepsilon ^{\frac{\eta _2}{2}}\Vert \nabla e^{n+1}\Vert \Vert e^{n+1}\Vert _{-1} + C\kappa _0^2 \varepsilon ^{- (2\alpha _0 + 6 + d + \frac{\eta _2}{2})} \tau ^{3 -\frac{d}{4}} \nonumber \\ \le \,&C\Vert e^{n+1}\Vert _{-1}^2 + \frac{1}{16} \varepsilon ^{\eta _2}\Vert \nabla e^{n+1}\Vert ^2 + C\kappa _0^2 \varepsilon ^{- (2\alpha _0 + 6 + d + \frac{\eta _2}{2})} \tau ^{3 -\frac{d}{4}}. \end{aligned}$$

(4.32)

It follows from (4.16)–(4.17) that

$$\begin{aligned} T_3 \lesssim \,&\frac{1}{2\varepsilon }\big | \big (g(u^n), e^{n+1} \big )\big | \sum _{i = 1}^n \Vert u^i - u^{i - 1}\Vert ^2 \nonumber \\ \le \,&\frac{1}{\varepsilon }\big | \big (g(u^n), e^{n+1} \big )\big | \sum _{i = 1}^n \Vert e^i - e^{i - 1}\Vert ^2 + \frac{1}{\varepsilon }\big | \big (g(u^n), e^{n+1} \big )\big | \sum _{i = 1}^n \Vert u(t_i) - u(t_{i - 1}) \Vert ^2 \nonumber \\ \le \,&C\varepsilon ^{-1} \Vert g(u^n)\Vert \Vert e^{n+1}\Vert \Big (\kappa _0\varepsilon ^{-(\alpha _0 + 2)} \tau ^{\frac{3}{2}} + \varepsilon ^{-\rho _2} \tau \Big ) \nonumber \\ \le \,&C\varepsilon ^{-1} \Vert e^{n+1}\Vert _{-1}^{\frac{1}{2}} \Vert \nabla e^{n+1}\Vert ^{\frac{1}{2}} \Big (\kappa _0\varepsilon ^{-(\alpha _0 + 2)} \tau ^{\frac{3}{2}} + \varepsilon ^{-\rho _2}\tau \Big ) \nonumber \\ \le \,&\varepsilon ^{\frac{\eta _2}{2}}\Vert e^{n+1}\Vert _{-1}\Vert \nabla e^{n+1}\Vert + C \kappa _0^2 \varepsilon ^{-(2\alpha _0 + 6 + \frac{\eta _2}{2})} \tau ^3 + C \varepsilon ^{-(2\rho _2 + 2 + \frac{\eta _2}{2}) }\tau ^2 \nonumber \\ \le \,&C\Vert e^{n+1}\Vert _{-1}^2 + \frac{1}{16} \varepsilon ^{\eta _2}\Vert \nabla e^{n+1}\Vert ^2 + C \kappa _0^2 \varepsilon ^{-(2\alpha _0 + 6 + \frac{\eta _2}{2}) }\tau ^3 + C \varepsilon ^{-(2\rho _2 + 2 + \frac{\eta _2}{2}) }\tau ^2 . \end{aligned}$$

(4.33)

In order to estimate the term $T_4$, we divided it into two parts as

$$\begin{aligned} T_4 = \,&- \frac{1}{2\varepsilon } \big (g(u^n), e^{n+1}\big ) \big (g(u^n), e^{n+1} - e^n \big ) \nonumber \\&+ \frac{1}{2\varepsilon } \big (g(u^n), e^{n+1}\big ) \big (g(u^n), u(t_{n+1}) - u(t_n) \big ) \nonumber \\ =:\,&T_{41} + T_{42}. \end{aligned}$$

(4.34)

Similarly to the estimate of $T_1$, the term $T_{41}$ can be bounded by

$$\begin{aligned} T_{41} \le \,&C\varepsilon ^{-1} \Vert g(u^n)\Vert ^2 \Vert e^{n+1}\Vert \Vert e^{n+1} - e^n\Vert \nonumber \\ \le \,&C\varepsilon ^{-1} \Vert e^{n+1} - e^n\Vert _{-1}^{\frac{1}{2}} \Vert \nabla e^{n+1} - \nabla e^n\Vert ^{\frac{1}{2}} \Vert e^{n+1}\Vert _{-1}^{\frac{1}{2}} \Vert \nabla e^{n+1}\Vert ^{\frac{1}{2}} \nonumber \\ \le \,&\tau ^{-\frac{1}{2}} \Vert e^{n+1} - e^n\Vert _{-1}\Vert e^{n+1}\Vert _{-1} + C \varepsilon ^{-2}\tau ^{\frac{1}{2}} \Vert \nabla e^{n+1} - \nabla e^n\Vert \Vert \nabla e^{n+1}\Vert \nonumber \\ \le \,&\frac{1}{2} \gamma _0\tau ^{-1} \Vert e^{n+1} - e^n\Vert _{-1}^2 + C \frac{1}{\gamma _0} \Vert e^{n+1}\Vert _{-1}^2 + \frac{1}{16} \varepsilon ^{\eta _2} \Vert \nabla e^{n+1}\Vert ^2 + T_1^*, \end{aligned}$$

(4.35)

which together with the condition $\tau \le \tilde{C}_1 \varepsilon ^{\eta _2 + 8}$ yields

$$\begin{aligned} \tau \sum _{n = 0}^N T_1^* = C\varepsilon ^{-(\eta _2 + 4)}\tau ^2 \sum _{n = 0}^N \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2 \le \frac{1}{16} \varepsilon ^4 \tau \sum _{n = 0}^N \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2. \end{aligned}$$

In addition, we obtain

$$\begin{aligned} T_{42} \le \,&C \varepsilon ^{-1} \Vert \nabla e^{n+1}\Vert ^{\frac{1}{2}} \Vert e^{n+1}\Vert _{-1}^{\frac{1}{2}} \Vert u(t_{n+1}) - u(t_n)\Vert \nonumber \\ \le \,&\varepsilon ^{\frac{\eta _2}{2}} \Vert \nabla e^{n+1} \Vert \Vert e^{n+1} \Vert _{-1} + C \varepsilon ^{- (2 + \frac{\eta _2}{2})}\tau \int _{t_n}^{t_{n+1}}\Vert u_t\Vert ^2 \, {\mathrm d}s \nonumber \\ \le \,&C\Vert e^{n+1}\Vert _{-1}^2 + \frac{1}{16}\varepsilon ^{\eta _2} \Vert \nabla e^{n+1}\Vert ^2 + C \varepsilon ^{- (2 + \frac{\eta _2}{2})}\tau \int _{t_n}^{t_{n+1}}\Vert u_t\Vert ^2 \, {\mathrm d}s. \end{aligned}$$

(4.36)

For $T_5$, using Cauchy-Schwartz inequality and Poincaré’s inequality for $e^{n+1} \in L_0^2({\varOmega })$, we have

$$\begin{aligned} T_5 \le \Vert \mathcal {R}^{n+1}\Vert _{H^{-1}}\Vert e^{n+1}\Vert _{H^1} \le 8\varepsilon ^{-\eta _2} \Vert \mathcal {R}^{n+1}\Vert _{H^{-1}}^2 + \frac{1}{16} \varepsilon ^{\eta _2} \Vert \nabla e^{n+1}\Vert ^2. \end{aligned}$$

(4.37)

Combining the estimates (4.30)–(4.37) together with (4.26), taking the summation for $n = 0$ to N, and assuming that $\tau \le \tilde{C}_1 \varepsilon ^{\eta _2 + 8}$, we have

$$\begin{aligned} \frac{1}{2} \Vert e^{N+1}\Vert _{-1}^2&+ \frac{1}{2} ( 1 - 2\gamma _0) \sum _{n = 0}^N \Vert e^{n+1} - e^n\Vert _{-1}^2 + \varepsilon \tau \sum _{n = 0}^N \Vert \nabla e^{n+1}\Vert ^2 \nonumber \\ \le \,&\Big [C + (1-\eta ) \lambda _0 + \frac{C\eta }{\varepsilon ^2 \eta _1}\Big ] \tau \sum _{n=0}^N \Vert e^{n+1}\Vert _{-1}^2 + 8\varepsilon ^{-\eta _2} \tau \sum _{n = 0}^{N} \Vert \mathcal {R}^{n+1}\Vert _{H^{-1}}^2 \nonumber \\&+ \Big [ (1-\eta )\varepsilon + \frac{\eta \eta _1}{4} + \frac{\varepsilon ^{\eta _2}}{2} \Big ] \tau \sum _{n=0}^N \Vert \nabla e^{n+1}\Vert ^2 + \frac{1}{8} \varepsilon ^4 \tau \sum _{n = 0}^N\Vert \nabla e^{n+1} - \nabla e^n\Vert ^2 \nonumber \\&+ C\kappa _0^2 \varepsilon ^{- (2\alpha _0 + 6 + d + \frac{\eta _2}{2})} \tau ^{3 -\frac{d}{4}} + C \kappa _0^2 \varepsilon ^{-(2\alpha _0 + 6 + \frac{\eta _2}{2}) }\tau ^3 \nonumber \\&+ C \varepsilon ^{-(2\rho _2 + 2 + \frac{\eta _2}{2}) }\tau ^2 + C \varepsilon ^{- (2 + \frac{\eta _2}{2})} \tau ^2 \int _{0}^{T}\Vert u_t\Vert ^2 \, {\mathrm d}s. \end{aligned}$$

(4.38)

By taking $\eta = \varepsilon ^3$, $\eta _1 = \varepsilon$ and $\eta _2 = 4$, we have

$$\begin{aligned} (1-\eta ) \varepsilon + \frac{\eta \eta _1}{4} + \frac{\varepsilon ^{\eta _2}}{2} = \varepsilon - \frac{1}{4} \varepsilon ^4, \end{aligned}$$

(4.39)

which together with (3.2c) gives

$$\begin{aligned} \Vert e^{N+1}\Vert _{-1}^2&+ (1 - 2\gamma _0) \sum _{n = 0}^N \Vert e^{n+1} - e^n\Vert _{-1}^2 + \frac{1}{2} \varepsilon ^4 \tau \sum _{n = 0}^N \Vert \nabla e^{n+1}\Vert ^2 \nonumber \\ \le \,&C\tau \sum _{n=0}^N \Vert e^{n+1}\Vert _{-1}^2 + \frac{1}{4} \tau \varepsilon ^4 \sum _{n = 0}^N \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2 + C\varepsilon ^{-4} \tau \sum _{n = 0}^{N} \Vert \mathcal {R}^{n+1}\Vert _{H^{-1}}^2 \nonumber \\&+ C\kappa _0^2 \varepsilon ^{- (2\alpha _0 + 8 + d )} \tau ^{3 -\frac{d}{4}} + C \kappa _0^2 \varepsilon ^{-(2\alpha _0 + 8) }\tau ^3 \nonumber \\&+ C\varepsilon ^{-(2\rho _2 + 4) }\tau ^2 + C\varepsilon ^{-(\rho _2 + 4) }\tau ^2 . \end{aligned}$$

(4.40)

Step 4: Completion of the proof. We now conclude the proof by the following induction argument which is based on the results from Step 1 to Step 3. By multiplying $\tau \varepsilon ^4$ on both sides of (4.24), combining the estimate (4.40), and together with Lemma 3.1, we obtain

$$\begin{aligned} \Vert e^{N+1}\Vert _{-1}^2&+ \tau \varepsilon ^4\Vert \nabla e^{N+1}\Vert ^2 + (1 - C \tau - 6\gamma _0)\sum _{n = 0}^N \Vert e^{n+1} - e^n\Vert _{-1}^2 \nonumber \\&+ (\frac{3}{4} - 6\gamma _0)\tau \varepsilon ^4 \sum _{n = 0}^N \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2 \nonumber + \frac{1}{2} \varepsilon ^4 \tau \sum _{n = 0}^N \Vert \nabla e^{n+1}\Vert ^2 \nonumber \\ \le \,&C_0\tau \sum _{n=0}^N \Big (\Vert e^n\Vert _{-1}^2 + \tau \varepsilon ^4 \Vert \nabla e^n\Vert ^2 \Big ) + C_1\kappa _0^2 \varepsilon ^{- (2\alpha _0 + 8 + d )} \tau ^{3 -\frac{d}{4}} \nonumber \\&+ C_2 \kappa _0^2 \varepsilon ^{-(2\alpha _0 + 8) }\tau ^3 + C_3\kappa _0^2 \varepsilon ^{- (2\alpha _0 + 2)}\tau ^{\frac{7}{2}} + C_4\varepsilon ^{-\max \{ \rho _1 + 3, \, 2\rho _2 + 4, \, \rho _2 + 6, \, \rho _3 + 4\} }\tau ^2. \end{aligned}$$

(4.41)

in which the term $\kappa _0' \varepsilon ^{- \max \{\rho _1 + 3, \,\rho _2+3, \,\rho _3 +1\}} \tau ^2$ is absorbed in $C_4\varepsilon ^{-\max \{ \rho _1 + 3, \, 2\rho _2 + 4, \, \rho _2 + 6, \, \rho _3 + 4\} }\tau ^2$.

Suppose that for sufficiently small constant $\gamma _0$ satisfying $\frac{3}{4} - 6\gamma _0 \ge \frac{1}{2}$ and sufficiently small $\tau$ satisfying

$$\begin{aligned} \tau \le \Big ( \frac{C_4}{C_1} \kappa _0^{-2} \Big )^{\frac{4}{4-d}} \varepsilon ^{\frac{4\alpha _0 + 32 + 4d}{4-d}}, \quad \tau \le \Big ( \frac{C_4}{C_2} \kappa _0^{-2} \Big ) \varepsilon ^{\alpha _0 + 8}, \quad \tau \le \Big ( \frac{C_4}{C_3} \kappa _0^{-2}\Big )^{\frac{2}{3}} \varepsilon ^{\frac{ (2\alpha _0 + 4)}{3}}, \end{aligned}$$

then, by denoting $\alpha _0:= \max \{ \rho _1 + 3, \, 2\rho _2 + 4, \, \rho _2 + 6, \, \rho _3 + 4\}$, we derive

$$\begin{aligned} \Vert e^{N+1}\Vert _{-1}^2&+ \tau \varepsilon ^4\Vert \nabla e^{N+1}\Vert ^2 + \frac{1}{2}\sum _{n = 0}^N \Vert e^{n+1} - e^n\Vert _{-1}^2 + \frac{1}{2} \varepsilon ^4 \sum _{n = 0}^N \tau \Vert \nabla e^{n+1}\Vert ^2 \nonumber \\&+ \frac{1}{2}\tau \varepsilon ^4 \sum _{n = 0}^N \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2 \le \, C_0\tau \sum _{n=0}^N \Big (\Vert e^n\Vert _{-1}^2 + \tau \varepsilon ^4 \Vert \nabla e^n\Vert ^2 \Big ) + 4 C_4\varepsilon ^{- \alpha _0 }\tau ^2. \end{aligned}$$

(4.42)

We denote $\kappa _0:= 4C_4 e^{(C_0 T)}$ and use the Gronwall’s inequality to get

$$\begin{aligned} \Vert e^{N+1}\Vert _{-1}^2&+ \frac{1}{2}\sum _{n = 0}^N \Vert e^{n+1} - e^n\Vert _{-1}^2 + \frac{1}{2} \varepsilon ^4 \sum _{n = 0}^N \tau \Vert \nabla e^{n+1}\Vert ^2 + \tau \varepsilon ^4\Vert \nabla e^{N+1}\Vert ^2 \nonumber \\&+ \frac{1}{2}\tau \varepsilon ^4 \sum _{n = 0}^m \Vert \nabla e^{n+1} - \nabla e^n\Vert ^2 \le \, 4C_4 e^{(C_0 T)} \varepsilon ^{-\alpha _0} \tau ^2 = \kappa _0 \varepsilon ^{-\alpha _0} \tau ^2, \end{aligned}$$

(4.43)

where $\kappa _0$ is independent of $\tau$ when $\tau$ is sufficiently small. The induction is completed.

In the above proof, we have used these conditions:

$$\begin{aligned} \tau \le \tilde{C}_1 \varepsilon ^{12}, \quad \tau \le \tilde{C}_2 \varepsilon ^{\frac{4\alpha _0 + 32 + 4d}{4-d}}, \quad \tau \le \tilde{C}_3\varepsilon ^{\alpha _0 + 8}, \quad \tau \le \tilde{C}_4 \varepsilon ^{\frac{ (2\alpha _0 + 4)}{3}}, \end{aligned}$$

(4.44)

where $\tilde{C}_2 = \Big ( \frac{C_4}{C_1} \kappa _0^{-2} \Big )^{\frac{4}{4-d}}$, $\tilde{C}_3 = \Big ( \frac{C_4}{C_2} \kappa _0^{-2} \Big )$, and $\tilde{C}_4 = \Big ( \frac{C_4}{C_3} \kappa _0^{-2}\Big )^{\frac{2}{3}}$. By denoting

$$\begin{aligned} \beta _0 = \frac{4\alpha _0 + 32 + 4d}{4-d} \qquad \text{ and } \qquad \tilde{C} = \min \{\tilde{C}_1, \tilde{C}_2, \tilde{C}_3, \tilde{C}_4 \}, \end{aligned}$$

(4.45)

we specify the final condition on $\tau$, that is, $\tau \le \tilde{C} \varepsilon ^{ \beta _0}$. $\square$

5 Numerical Experiments

In this section, we present a two-dimensional numerical test to validate the theoretical results on the energy decay properties proved in Theorem 3.1, as well as the convergence rates of the proposed method given in Theorem 4.1. All the computations are performed using the software package NGSolve (https://ngsolve.org).

We solve the Cahn-Hilliard equation (1.1) on the two-dimensional square ${\varOmega }= [0, 1] \times [0, 1]$ under Neumann boundary conditions by using the proposed scheme (2.8) with the following initial condition

$$\begin{aligned} u_0(x, y) =&\tanh \big (((x - 0.65)^2 + (y - 0.5)^2 - 0.1^2)/\varepsilon \big )\nonumber \\&\times \tanh \big (((x - 0.35)^2 + (y - 0.5)^2 - 0.125^2)/\varepsilon \big ), \end{aligned}$$

(5.1)

where $\tanh (x):= (e^x - e^{-x})/(e^x + e^{-x})$. This type of initial condition is also adopted in [21, 26], where the set of the zero-level of the initial function $u_0(x, y)$ encloses two circles of radius 0.1 and 0.125, respectively.

To obtain a $C^4$ potential function F(v) that satisfies the assumption 3.1, we modify the common double-well potential $F(v) = \frac{1}{4} (v^2 - 1)^2$ by setting $M = 2$ in (3.3) to get a cut-off function $\hat{F}(v) \in C^4(\mathbb {R})$. Correspondingly, the ninth-order polynomials $\Phi _+(v)$ and $\Phi _-(v)$ in (3.3) are determined with the following conditions

$$\begin{aligned} \left\{ \begin{aligned}&\Phi _+^{(i)}(M) = F^{(i)}(M) \quad \text{ and } \quad \Phi _-^{(i)}(-M) = F^{(i)}(-M) \qquad \text{ for } \,\,\, i = 0, 1, 2, 3, 4,\\&\Phi _+(2M) = \Phi _-(-2M) = \frac{1}{4} ((2M)^2 - 1)^2, \\&\Phi _+^{(1)}(2M) = \Phi _-^{(1)}(-2M) = (2M)^3 - 2M, \\&\Phi _+^{(i)}(2M) = \Phi _-^{(i)}(-2M) = 0 \qquad \text{ for } \,\,\, i = 2, 3, 4. \end{aligned} \right. \end{aligned}$$

(5.2)

Note that the truncation point $M = 2$ used here are for convenience only. For simplicity, we still denote the modified function $\hat{F}(u)$ by F(u).

The spatial discretization is done by using the Galerkin finite element method. Let $S_h$ denotes the $P_s$ conforming finite element space defined by

$$\begin{aligned} S_h := \{ v_h \in C(\bar{\varOmega }); \, \, v_h|_{K} \in P_s(K), \, \,\, \forall \, K \in \mathcal {T}_h\}, \end{aligned}$$

where $\mathcal {T}_h$ is a quasi-uniform triangulation of ${\varOmega }$. We introduce space notation $S_h^{\circ }:= \{ v_h \in S_h; \, \, (v_h, 1) = 0 \}$, and define the discrete inverse Laplace operator $- \Delta _h^{-1}: L_0^2({\varOmega }) \rightarrow S_h^{\circ }$ such that

$$\begin{aligned} \big (\nabla (- \Delta _h^{-1} ) v, \nabla \eta _h \big ) = (v, \eta _h) \qquad \forall \, \, \, \eta _h \in S_h. \end{aligned}$$

(5.3)

Since the exact solution of the considered problem is not known, we compute the orders of convergence by the formula

$$ \text{ order } \text{ of } \text{ convergence } = \log \Bigg (\frac{\Vert u_N^{(\tau )}-u_N^{(\tau /2)}\Vert _{h, -1}}{\Vert u_N^{(\tau /2)}-u_N^{(\tau /4)}\Vert _{h, -1}}\Bigg )/\log (2) $$

based on the finest three meshes, where $u_N^{(\tau )}$ denotes the numerical solution at $t_N=T$ computed by using a stepsize $\tau$, and $\Vert v\Vert _{h, -1}:= \sqrt{(v, - \Delta _h^{-1} v)}$ for $v \in L_0^2({\varOmega })$.

The time discretization errors in $\Vert \cdot \Vert _{h, -1}$-norm are presented in Fig. 1 (left) for four different $\varepsilon = 0.08, 0.06, 0.05, 0.04$ at $T = 0.005$, where we have used finite elements of degree $s = 3$ with a sufficiently spatial mesh $h = 1/64$ so that the error from spatial discretization is negligibly small in observing the temporal convergence rates. From Fig. 1 (left), we see that the error of time discretization is $O(\tau )$, which is consistent with the theoretical results proved in Theorem 4.1. In addition, Fig. 1 (right) shows the evolution in time of the discrete SAV energy for four different $\varepsilon$, which should be decreasing according to Theorem 3.1. This graph clearly confirms this decay property. Therefore, the numerical experiments are in accordance with our theoretical results.

Fig. 1

(left) Time discretization errors;    (right) Evolution of the SAV energy

Figure 2 shows snapshots of the numerical interface for four different $\varepsilon =0.08, 0.06, 0.05, 0.04$ at six fixed time points. They clearly indicate that at each time point, as $\varepsilon$ tends to zero, the numerical interface converges to the sharp interface of the Hele-haw flow, which is consistent with the phenomenon stated in [21, 26]. It also shows that for larger $\varepsilon$, the numerical interface evolves faster in time.

Fig. 2

Snapshots of the zero-level sets of the numerical solutions