∎

11institutetext: L. Ju 22institutetext: Department of Mathematics, University
of South Carolina, Columbia, SC 29208, USA
22email<EMAIL_ADDRESS>33institutetext: X. Li 44institutetext: Corresponding
author. Department of Applied Mathematics, The Hong Kong Polytechnic
University, Hung Hom, Kowloon, Hong Kong
44email<EMAIL_ADDRESS>55institutetext: Z. Qiao 66institutetext:
Department of Applied Mathematics & Research Institute for Smart Energy, The
Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong
66email<EMAIL_ADDRESS>
# Stabilized exponential-SAV schemes preserving energy dissipation law and
maximum bound principle for the Allen–Cahn type equations ††thanks: This work
is supported by the CAS AMSS-PolyU Joint Laboratory of Applied Mathematics. L.
Ju’s work is partially supported by US National Science Foundation grant
DMS-2109633 and US Department of Energy grant DE-SC0020270. X. Li’s work is
partially supported by National Natural Science Foundation of China grant
11801024 and the Hong Kong Polytechnic University grants 4-ZZMK and 1-BD8N. Z.
Qiao’s work is partially supported by the Hong Kong Research Council RFS grant
RFS2021-5S03 and GRF grants 15300417 and 15302919.

Lili Ju  Xiao Li  Zhonghua Qiao

(Received: date / Accepted: date)

###### Abstract

It is well-known that the Allen–Cahn equation not only satisfies the energy
dissipation law but also possesses the maximum bound principle (MBP) in the
sense that the absolute value of its solution is pointwise bounded for all
time by some specific constant under appropriate initial/boundary conditions.
In recent years, the scalar auxiliary variable (SAV) method and many of its
variants have attracted much attention in numerical solution for gradient flow
problems due to their inherent advantage of preserving certain discrete
analogues of the energy dissipation law. However, existing SAV schemes usually
fail to preserve the MBP when applied to the Allen–Cahn equation. In this
paper, we develop and analyze new first- and second-order stabilized
exponential-SAV schemes for a class of Allen–Cahn type equations, which are
shown to simultaneously preserve the energy dissipation law and MBP in
discrete settings. In addition, optimal error estimates for the numerical
solutions are rigorously obtained for both schemes. Extensive numerical tests
and comparisons are also conducted to demonstrate the performance of the
proposed schemes.

###### Keywords:

maximum bound principle energy dissipation stabilized method exponential
scalar auxiliary variable

###### MSC:

35K55 65M12 65M15 65F30

##  1 Introduction

Let us consider a class of reaction-diffusion equations taking the following
form

$u_{t}=\varepsilon^{2}\Delta u+f(u),\quad t>0,\ \bm{x}\in\Omega,$ (1.1)

where $\Omega\subset\mathbb{R}^{d}$ is a spatial domain,
$u=u(t,\bm{x}):[0,\infty)\times\overline{\Omega}\to\mathbb{R}$ is the unknown
function, $\varepsilon>0$ denotes an interfacial parameter, and $f(u)$ is a
nonlinear reaction term with $f$ being continuously differentiable. We also
impose the initial condition

$u(0,\cdot)=u_{\text{\rm init}}\quad\text{on }\overline{\Omega}$

and the periodic or homogeneous Neumann boundary conditions. The equation
(1.1) usually can be regarded as the $L^{2}$ gradient flow with respect to the
energy functional

$E(u)=\int_{\Omega}\Big{(}\frac{\varepsilon^{2}}{2}|\nabla
u(\bm{x})|^{2}+F(u(\bm{x}))\Big{)}\,\mathrm{d}\bm{x},$ (1.2)

where $F$ is a smooth potential function satisfying $F^{\prime}=-f$, and thus,
the solution to the equation (1.1) decreases the energy (1.2) along with the
time, i.e., $\frac{\mathrm{d}}{\mathrm{d}t}E(u(t))\leq 0$, which is often
called the energy dissipation law. In addition, we also assume that

$\text{there exists a constant $\beta>0$ such that $f(\beta)\leq 0\leq
f(-\beta)$}.$ (1.3)

It has been proved in DuJuLiQi21  that the equation (1.1) satisfies the
maximum bound principle (MBP) in the sense that if the absolute value of the
initial data is bounded pointwise by $\beta$, then the absolute value of the
solution is also bounded by $\beta$ pointwise for all time, i.e.,

$\max_{\bm{x}\in\overline{\Omega}}|u_{\text{\rm
init}}(\bm{x})|\leq\beta\quad\Longrightarrow\quad\max_{\bm{x}\in\overline{\Omega}}|u(t,\bm{x})|\leq\beta,\quad\forall\,t>0.$
(1.4)

An important and special case of (1.1) is the Allen–Cahn equation with
$f(u)=u-u^{3}$, which was originally introduced in AlCa79  to model the motion
of anti-phase boundaries in crystalline solids. The solution represents the
difference between the concentrations of two components of the alloy and thus
should be evaluated between $-1$ and $1$, which is guaranteed by the MBP. With
the corresponding double-well potential $F(u)=\frac{1}{4}(u^{2}-1)^{2}$, the
associated energy functional (1.2) decays in time, which reflects the energy
dissipation of the phase transition process. Both the MBP and the energy
stability are also satisfied by some variants of (1.1), such as the nonlocal
Allen–Cahn equation for phase separations within long-range interactions
Bates06 ; DuJuLiQi19  and the fractional Allen–Cahn equation used to describe
some anomalous diffusion processes DuYaZh20 ; GuiZh15 . To obtain stable
numerical simulations and avoid nonphysical solutions for these models, it is
highly desirable to design numerical schemes preserving effectively these two
basic physical properties, the MBP and the energy dissipation law in time
discrete settings.

In the past decades, there has been a large amount of research denoted to
energy-stable numerical schemes for time discretization of gradient flow
equations, such as convex splitting schemes GuWaWi14 ; ShWaWaWi12 ; WiWaLo09 ,
stabilized semi-implicit schemes FeTaYa13 ; ShYa10b ; XuTa06 , and exponential
time differencing (ETD) schemes DuJuLiQi19 ; JuLiQiZh18 ; JuZhDu15 . More
recently, invariant energy quadratization (IEQ) schemes XuYaZhXi19 ; Yang16 ;
YangZh20  and scalar auxiliary variable (SAV) schemes ShXu18 ; ShXuYa18 ;
ShXuYa19  were proposed to naturally provide energy-stable and linear
algorithms with second-order temporal accuracy. While the main idea for both
methods is to reformulate and split the energy functional (1.2) in the
quadratic form by introducing extra variables, the SAV approach is usually
more efficient in terms of computations. Many variants of SAV schemes were
developed later; see AkrivisLiLi19 ; ChenYa19 ; ChengLiSh20 ; ChengLiSh21 ;
HouAzXu19 ; HuangShYa20 ; LiuLi20  and the references therein. In practice, a
suitable stabilization term is also introduced in such splitting in order to
maintain numerical stability for highly stiff problems. On the other hand,
existing SAV-type schemes usually fail to preserve the MBP, and a special case
is the auxiliary variable proposed in HuangShYa20  which is shown to be
positivity-preserving.

The MBP preservation recently has also attracted increasingly attention in the
field of numerical methods for the Allen–Cahn type equations of the form
(1.1). The semi-implicit schemes were extensively studied in, e.g. HoLe20 ;
HoTaYa17 ; LiaoTaZh20 ; ShTaYa16 ; TaYa16 ; XiFeYu17 , for the classic,
fractional, or surface Allen–Cahn equations. The first- and second-order
stabilized ETD schemes were shown to preserve the MBP unconditionally for the
nonlocal Allen–Cahn equation DuJuLiQi19  and the conservative Allen–Cahn
equation LiJuCaFe21 . An abstract framework on MBP preservation of the ETD
schemes for a class of semilinear parabolic equations was established in
DuJuLiQi21 , where sufficient conditions for the linear and nonlinear
operators are presented in order to guarantee the MBP. A family of stabilized
integrating factor Runge-Kutta (IFRK) schemes, up to third order, were
developed in LiLiJuFe21 , which can unconditionally preserve the MBP. In
addition, a fourth-order (conditionally) MBP-preserving IFRK scheme was
presented in JuLiQiYa21 . So far, as one of the very popular methods, there is
still not much systematical study on MBP-preserving SAV schemes.

The main goal of this paper is to develop first- and second-order energy
dissipative and MBP-preserving SAV-type schemes for the Allen–Cahn type
equation (1.1) by using an appropriate stabilization technique. Specifically,
we propose new stabilized exponential-SAV (ESAV) schemes by introducing an
artificial stabilization term rather than basing on the splitting of the
energy functional suggested in ShXuYa19 . With the effect of such
stabilization, we show that the proposed first-order scheme preserves the MBP
unconditionally with an appropriate stabilizing parameter and the second-order
one does under a time step size constraint. A main difficulty for numerical
analysis of the two schemes lies in that the coefficients of the nonlinear
term and the stabilization term are varying rather than constant due to the
use of the SAV approach. With the help of the energy dissipation and MBP, we
are able to show that such variable coefficients are bounded from above and
below by certain positive constants, and consequently, optimal error estimate
are successfully obtained for the proposed schemes. To the best of our
knowledge, this is the first work in the direction of designing such SAV-type
methods. More importantly, the proposed stabilizing approaches can be easily
generalized to deal with many other type of gradient flow problems where the
SAV methods apply.

The rest of this paper is organized as follows. Section 2 is devoted to
spatial discretization of the equation (1.1) and a brief summary of the
classic SAV and ESAV schemes for the time integration. Then, our first- and
second-order stabilized ESAV schemes are presented in Section 3, together with
the energy dissipation law, MBP preservation, and convergence analysis of the
resulting fully discrete systems. In Section 4, extensive numerical tests and
comparisons are carried out to demonstrate the performance of the proposed
schemes. Some concluding remarks are finally given in Section 5.

##  2 Spatial discretization and SAV schemes for time integration

For simplicity, throughout this paper we consider the two-dimensional square
domain $\Omega=(0,L)\times(0,L)$ for the equation (1.1) equipped with periodic
boundary conditions. Note that the extensions to three-dimension problems and
homogeneous Neumann boundary condition are straightforward. For other feasible
spatial discretization, we refer to DuJuLiQi21  for more details. In this
section, we first present some notations related to the spatial discretization
by central finite difference, then briefly review the classic SAV schemes for
time integration.

###  2.1 Spatial discretization and the space-discrete problem

Given a positive integer $M$, we set $h=L/M$ to be the size of the uniform
mesh partitioning $\overline{\Omega}$. Denote by $\Omega_{h}$ the set of mesh
points $(x_{i},y_{j})=(ih,jh)$, $1\leq i,j\leq M$. For a grid function $v$
defined on $\Omega_{h}$, we write $v_{ij}=v(x_{i},y_{j})$ for simplicity. Let
$\mathcal{M}_{h}$ be the set of all periodic grid functions on $\Omega_{h}$,
i.e.,

$\mathcal{M}_{h}=\\{v:\Omega_{h}\to\mathbb{R}\,|\,v_{i+kM,j+lM}=v_{ij},\
k,l\in\mathbb{Z},\ 1\leq i,j\leq M\\}.$

The discrete inner product $\langle\cdot,\cdot\rangle$, discrete $L^{2}$ norm
$\|\cdot\|$, and discrete $L^{\infty}$ norm $\|\cdot\|_{\infty}$ can be
defined as usual, namely,

$\langle
v,w\rangle=h^{2}\sum_{i,j=1}^{M}v_{ij}w_{ij},\qquad\|v\|=\sqrt{\langle
v,v\rangle},\qquad\|v\|_{\infty}=\max_{1\leq i,j\leq M}|v_{ij}|$

for any $v,w\in\mathcal{M}_{h}$, and

$\langle\bm{v},\bm{w}\rangle=\langle v^{1},w^{1}\rangle+\langle
v^{2},w^{2}\rangle,\qquad\|\bm{v}\|=\sqrt{\langle\bm{v},\bm{v}\rangle}$

for any
$\bm{v}=(v^{1},v^{2})^{T},\bm{w}=(w^{1},w^{2})^{T}\in\mathcal{M}_{h}\times\mathcal{M}_{h}$.
We apply the second-order central finite difference to approximate spatial
differentiation operators. For any $v\in\mathcal{M}_{h}$, the discrete Laplace
operator $\Delta_{h}$ is defined by

$\Delta_{h}v_{ij}=\frac{1}{h^{2}}(v_{i+1,j}+v_{i-1,j}+v_{i,j+1}+v_{i,j-1}-4v_{ij}),\quad
1\leq i,j\leq M,$

and the discrete gradient operator $\nabla_{h}$ is defined by

$\nabla_{h}v_{ij}=\Big{(}\frac{v_{i+1,j}-v_{ij}}{h},\frac{v_{i,j+1}-v_{ij}}{h}\Big{)}^{T},\quad
1\leq i,j\leq M.$

By periodic boundary conditions, the summation-by-parts formula is easy to
verify:

$\langle
v,\Delta_{h}w\rangle=-\langle\nabla_{h}v,\nabla_{h}w\rangle=\langle\Delta_{h}v,w\rangle,\quad\forall\,v,w\in\mathcal{M}_{h}.$

Obviously, $\Delta_{h}$ is self-adjoint and negative semi-definite. For any
function $\varphi:\overline{\Omega}\to\mathbb{R}$, we denote by $I_{h}$ the
operator projecting $\varphi$ on the mesh as
$(I_{h}\varphi)_{ij}=\varphi(x_{i},y_{j})$ for $1\leq i,j\leq M$. For example,
we have

$\max_{1\leq i,j\leq
N}|\Delta_{h}(I_{h}\varphi)_{ij}-\Delta\varphi(x_{i},y_{j})|\leq
C_{\varphi}h^{2},\quad\forall\,\varphi\in C_{\text{\rm
per}}^{4}(\overline{\Omega}).$

For simplicity, we may directly omit the notation $I_{h}$ when there is no
ambiguity.

Since $\mathcal{M}_{h}$ is a finite-dimensional linear space, any grid
function $v\in\mathcal{M}_{h}$ and any linear operator
$Q:\mathcal{M}_{h}\to\mathcal{M}_{h}$ can be regarded as a vector in
$\mathbb{R}^{M^{2}}$ and a matrix in $\mathbb{R}^{M^{2}\times M^{2}}$,
respectively. We still use the notations $\|\cdot\|$ and $\|\cdot\|_{\infty}$
to denote the matrix induced-norms consistent with $\|\cdot\|$ and
$\|\cdot\|_{\infty}$ defined for vectors before, respectively. By regarding
$\Delta_{h}$ as a linear operator, we know that $\Delta_{h}$ is the generator
of a contraction semigroup on $\mathcal{M}_{h}$ DuJuLiQi21 . Instead, by
viewing $\Delta_{h}$ as a matrix, it is weakly diagonally dominant with all
diagonal entries negative. Moreover, we have the following useful estimate and
the proof can be found in TaYa16 .

###### Lemma 1

For any $a>0$, we have $\|(aI-\Delta_{h})^{-1}\|_{\infty}\leq a^{-1}$, where
$I$ represents the identity matrix.

We have assumed that $f$ is continuously differentiable, so
$\|f^{\prime}\|_{C[-\beta,\beta]}$ is always finite and then the following
result is valid DuJuLiQi21 .

###### Lemma 2

Under the assumption (1.3), if $\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$
holds for some positive constant $\kappa$, then we have
$|f(\xi)+\kappa\xi|\leq\kappa\beta$ for any $\xi\in[-\beta,\beta]$.

Next, let us introduce the space-discrete version of (1.1). The space-discrete
problem is to find a function $u_{h}:[0,\infty)\to\mathcal{M}_{h}$ satisfies

$\frac{\mathrm{d}u_{h}}{\mathrm{d}t}=\varepsilon^{2}\Delta_{h}u_{h}+f(u_{h})$
(2.1)

with $u_{h}(0)=u_{\text{\rm init}}$. It is easy to verify the energy
dissipation law for (2.1) in the sense that

$\frac{\mathrm{d}}{\mathrm{d}t}E_{h}(u_{h}(t))\leq 0,$

where $E_{h}$ is the spatially-discretized energy functional defined as

$E_{h}(v):=\frac{\varepsilon^{2}}{2}\|\nabla_{h}v\|^{2}+\langle
F(v),1\rangle,\quad\forall\,v\in\mathcal{M}_{h}.$ (2.2)

According to DuJuLiQi21 , the MBP also holds for $u_{h}$, i.e.,
$\|u_{h}(t)\|_{\infty}\leq\beta$ for any $t>0$ if $\|u_{\text{\rm
init}}\|_{\infty}\leq\beta$.

Let us partition the time interval into $\\{t_{n}=n{\tau}\\}_{n\geq 0}$ with
${\tau}>0$ being a uniform time step size. In the remaining part of the paper,
we will study time integration schemes for the space-discrete system (2.1).
For simplicity of representation, we denote by $u^{n}$ the fully discrete
approximate value of $u_{e}(t_{n})$ or $u_{h,e}(t_{n})$ with $u_{e}$ and
$u_{h,e}$ denoting the exact solutions to the original continuous problem
(1.1) and the space-discrete problem (2.1), respectively. In general, for a
sequence $\\{v^{n}\\}$, we define the following notations:

$\delta_{t}v^{n+1}=\frac{v^{n+1}-v^{n}}{{\tau}},\qquad
v^{n+\frac{1}{2}}=\frac{v^{n+1}+v^{n}}{2}.$

###  2.2 Classic SAV schemes and stabilization

Here we give a brief summary of the classic SAV schemes. The main idea of SAV
is to reformulate the energy functional (1.2) in the quadratic form by
introducing an appropriate SAV. The framework of the classic SAV schemes is
based on a linear splitting of the energy functional and, as shown in ShXuYa19
, a suitable stabilization term is usually also introduced in such splitting
so that the numerical simulations can provide satisfactory results for highly
stiff problems in practice.

Denoting by $\kappa\geq 0$ the stabilizing constant, the energy functional
(1.2) can be rewritten with a stabilization term as

$\displaystyle E(u)$
$\displaystyle=\int_{\Omega}\Big{(}\frac{\varepsilon^{2}}{2}|\nabla
u|^{2}+\frac{\kappa}{2}u^{2}+F(u)-\frac{\kappa}{2}u^{2}\Big{)}\,\mathrm{d}\bm{x}$
$\displaystyle=\frac{\varepsilon^{2}}{2}\|\nabla
u\|_{L^{2}}^{2}+\frac{\kappa}{2}\|u\|_{L^{2}}^{2}+\int_{\Omega}\Big{(}F(u)-\frac{\kappa}{2}u^{2}\Big{)}\,\mathrm{d}\bm{x}.$
(2.3)

Suppose the last term in (2.2) is bounded from below, that is,

$E_{2}(u):=\int_{\Omega}\Big{(}F(u)-\frac{\kappa}{2}u^{2}\Big{)}\,\mathrm{d}\bm{x}\geq-
C_{0}$ (2.4)

for some constant $C_{0}\geq 0$. Choosing $\delta>C_{0}$, let us define an
auxiliary variable $r(t)=\sqrt{E_{2}(u(t))+\delta}$, and reformulate the
original problem (1.1) to the following equivalent system:

$\displaystyle u_{t}$ $\displaystyle=\varepsilon^{2}\Delta u-\kappa
u+\frac{r}{\sqrt{E_{2}(u)+\delta}}(f(u)+\kappa u),$ $\displaystyle r_{t}$
$\displaystyle=-\frac{1}{2\sqrt{E_{2}(u)+\delta}}(f(u)+\kappa u,u_{t}).$

Then the first-order SAV scheme (SAV1) is given by ShXuYa19

$\displaystyle\delta_{t}u^{n+1}$
$\displaystyle=\varepsilon^{2}\Delta_{h}u^{n+1}-\kappa
u^{n+1}+\frac{r^{n+1}}{\sqrt{E_{2h}(u^{n})+\delta}}(f(u^{n})+\kappa u^{n}),$
(2.6a) $\displaystyle\delta_{t}r^{n+1}$
$\displaystyle=-\frac{1}{2\sqrt{E_{2h}(u^{n})+\delta}}\langle f(u^{n})+\kappa
u^{n},\delta_{t}u^{n+1}\rangle,$ (2.6b)

and the Crank–Nicolson type second-order SAV scheme (SAV2) reads as ShXuYa19

$\displaystyle\delta_{t}u^{n+1}$
$\displaystyle=\varepsilon^{2}\Delta_{h}u^{n+\frac{1}{2}}-\kappa
u^{n+\frac{1}{2}}+\frac{r^{n+1}+r^{n}}{2\sqrt{E_{2h}(\widehat{u}^{n+\frac{1}{2}})+\delta}}(f(\widehat{u}^{n+\frac{1}{2}})+\kappa\widehat{u}^{n+\frac{1}{2}}),$
(2.7a) $\displaystyle\delta_{t}r^{n+1}$
$\displaystyle=-\frac{1}{2\sqrt{E_{2h}(\widehat{u}^{n+\frac{1}{2}})+\delta}}\langle
f(\widehat{u}^{n+\frac{1}{2}})+\kappa\widehat{u}^{n+\frac{1}{2}},\delta_{t}u^{n+1}\rangle,$
(2.7b)

where $\widehat{u}^{n+\frac{1}{2}}$ is generated by solving the system

$\frac{\widehat{u}^{n+\frac{1}{2}}-u^{n}}{{\tau}/2}=\varepsilon^{2}\Delta_{h}\widehat{u}^{n+\frac{1}{2}}+f(u^{n})-\kappa(\widehat{u}^{n+\frac{1}{2}}-u^{n}).$

Both (2.6) and (2.7) are linear schemes and energy dissipative in the sense
that $\overline{E}_{h}(u^{n+1},r^{n+1})\leq\overline{E}_{h}(u^{n},r^{n})$ with
respect to the following modified energy

$\overline{E}_{h}(u^{n},r^{n}):=\frac{\varepsilon^{2}}{2}\|\nabla_{h}u^{n}\|^{2}+\frac{\kappa}{2}\|u^{n}\|^{2}+(r^{n})^{2}-\delta.$

Note that in the discrete settings, $\overline{E}_{h}(u^{n},r^{n})$ is only an
approximation of the original discrete energy $E_{h}(u^{n})$ defined in (2.2)
and they are not equal in general since
$r^{n}\not=\sqrt{E_{2h}(u^{n})+\delta}$ for $n\geq 1$. However, the MBP cannot
be theoretically preserved by the above classic SAV schemes (2.6) and (2.7)
(See the discussion in Remark 3).

###  2.3 Exponential-SAV schemes

A variant of the classic SAV approach, called the exponential-SAV (ESAV)
scheme, was studied in LiuLi20 . We below summarize the ESAV method, also with
a stabilization term based on the energy splitting (2.2). Define the auxiliary
variable by $r(t)=\exp\\{E_{2}(u(t))\\}$ and reformulate (1.1) as

$\displaystyle u_{t}=\varepsilon^{2}\Delta u-\kappa
u+\frac{r}{\exp\\{E_{2}(u)\\}}(f(u)+\kappa u),$ $\displaystyle(\ln
r)_{t}=-\frac{r}{\exp\\{E_{2}(u)\\}}(f(u)+\kappa u,u_{t}).$

Then the first-order ESAV scheme (ESAV1) reads as

$\displaystyle\delta_{t}u^{n+1}=\varepsilon^{2}\Delta_{h}u^{n+1}-\kappa
u^{n+1}+\frac{r^{n}}{\exp\\{E_{2h}(u^{n})\\}}(f(u^{n})+\kappa u^{n}),$ (2.9a)
$\displaystyle\frac{\ln r^{n+1}-\ln
r^{n}}{{\tau}}=-\frac{r^{n}}{\exp\\{E_{2h}(u^{n})\\}}\langle f(u^{n})+\kappa
u^{n},\delta_{t}u^{n+1}\rangle.$ (2.9b)

Setting $\kappa=0$, the scheme (2.9) reduces exactly to the original ESAV
scheme (without stabilization) presented in LiuLi20 . The Crank–Nicolson type
ESAV scheme (ESAV2) is given by

$\displaystyle\delta_{t}u^{n+1}=\varepsilon^{2}\Delta_{h}u^{n+\frac{1}{2}}-\kappa
u^{n+\frac{1}{2}}+\frac{\widehat{r}^{n+\frac{1}{2}}}{\exp\\{E_{2h}(\widehat{u}^{n+\frac{1}{2}})\\}}(f(\widehat{u}^{n+\frac{1}{2}})+\kappa\widehat{u}^{n+\frac{1}{2}}),$
(2.10a) $\displaystyle\frac{\ln r^{n+1}-\ln
r^{n}}{{\tau}}=-\frac{\widehat{r}^{n+\frac{1}{2}}}{\exp\\{E_{2h}(\widehat{u}^{n+\frac{1}{2}})\\}}\langle
f(\widehat{u}^{n+\frac{1}{2}})+\kappa\widehat{u}^{n+\frac{1}{2}},\delta_{t}u^{n+1}\rangle,$
(2.10b)

where the value $(\widehat{u}^{n+\frac{1}{2}},\widehat{r}^{n+\frac{1}{2}})$
can be generated by an extrapolation as suggested in LiuLi20  or predicted by
the first-order scheme (2.9) with half of the time step size:

$\displaystyle\frac{\widehat{u}^{n+\frac{1}{2}}-u^{n}}{{\tau}/2}=\varepsilon^{2}\Delta_{h}\widehat{u}^{n+\frac{1}{2}}-\kappa\widehat{u}^{n+\frac{1}{2}}+\frac{r^{n}}{\exp\\{E_{2h}(u^{n})\\}}(f(u^{n})+\kappa
u^{n}),$ (2.11a) $\displaystyle\ln\widehat{r}^{n+\frac{1}{2}}-\ln
r^{n}=-\frac{r^{n}}{\exp\\{E_{2h}(u^{n})\\}}\langle f(u^{n})+\kappa
u^{n},\widehat{u}^{n+\frac{1}{2}}-u^{n}\rangle.$ (2.11b)

We will adopt (2.11) in the numerical experiments for the comparison. Both
(2.9) and (2.10) are energy dissipative in the sense that
$\widetilde{E}_{h}(u^{n+1},r^{n+1})\leq\widetilde{E}_{h}(u^{n},r^{n})$ with
respect to the following modified energy

$\widetilde{E}_{h}(u^{n},r^{n}):=\frac{\varepsilon^{2}}{2}\|\nabla_{h}u^{n}\|^{2}+\frac{\kappa}{2}\|u^{n}\|^{2}+\ln
r^{n}.$

Similar to the classic SAV schemes, the above ESAV schemes (2.9) and (2.10)
also cannot preserve the MBP (See the discussion in Remark 3).

##  3 New stabilized exponential-SAV schemes

From now on, we always assume the initial value $u_{\text{\rm init}}$ has the
enough regularity as needed. By spatial discretization, there is a constant
$h_{0}>0$, depending on $u_{\text{\rm init}}$, $F$, $\varepsilon$, such that

$E_{h}(u_{\text{\rm init}})\leq E(u_{\text{\rm
init}})+1,\quad\|\Delta_{h}u_{\text{\rm init}}\|\leq\|\Delta u_{\text{\rm
init}}\|_{L^{2}}+1,\qquad\forall\,h\in(0,h_{0}].$ (3.1)

The continuity of $F$ implies that $F$ is bounded from below on
$[-\beta,\beta]$. Therefore, according to the MBP (1.4), it holds that

$E_{1}(u):=\int_{\Omega}F(u)\,\mathrm{d}\bm{x}\geq-C_{*}$

for some constant $C_{*}\geq 0$. Introducing $s(t)=E_{1}(u(t))$, we then have
the following energy which is equivalent to $E(u)$:

$\mathcal{E}(u,s)=\frac{\varepsilon^{2}}{2}\|\nabla u\|_{L^{2}}^{2}+s.$

Partially inspired by the idea of ESAV method LiuLi20 , we rewrite the
equation (1.1) as the following equivalent system:

$\displaystyle u_{t}$ $\displaystyle=\varepsilon^{2}\Delta
u+\frac{\exp\\{s\\}}{\exp\\{E_{1}(u)\\}}f(u),$ $\displaystyle s_{t}$
$\displaystyle=-\frac{\exp\\{s\\}}{\exp\\{E_{1}(u)\\}}(f(u),u_{t}).$

The corresponding space-discrete problem is to find
$u_{h}(t)\in\mathcal{M}_{h}$ and $s_{h}(t)$ for $t>0$ satisfies

$\displaystyle\frac{\mathrm{d}u_{h}}{\mathrm{d}t}$
$\displaystyle=\varepsilon^{2}\Delta_{h}u_{h}+g(u_{h},s_{h})f(u_{h}),$ (3.3a)
$\displaystyle\frac{\mathrm{d}s_{h}}{\mathrm{d}t}$
$\displaystyle=-g(u_{h},s_{h})\Big{\langle}f(u_{h}),\frac{\mathrm{d}u_{h}}{\mathrm{d}t}\Big{\rangle},$
(3.3b)

where

$g(u_{h},s_{h}):=\frac{\exp\\{s_{h}\\}}{\exp\\{E_{1h}(u_{h})\\}}>0,$ (3.4)

and $E_{1h}$ denotes the space-discrete version of $E_{1}$, i.e.,
$E_{1h}(v):=\langle F(v),1\rangle$ for any $v\in\mathcal{M}_{h}$. Based on
such an equivalent form, we will give the stabilized ESAV schemes in the fully
discrete version. This section is devoted to the first-order scheme and the
second-order one will be discussed in the next section. Recall that we use
$u^{n}$ to represent the fully discrete approximate value of $u_{e}(t_{n})$,
the exact solution to the problem (1.1).

###  3.1 First-order sESAV scheme

The first-order stabilized ESAV fully-discrete scheme (sESAV1) is given by

$\displaystyle\delta_{t}u^{n+1}$
$\displaystyle=\varepsilon^{2}\Delta_{h}u^{n+1}+g(u^{n},s^{n})f(u^{n})-\kappa
g(u^{n},s^{n})(u^{n+1}-u^{n}),$ (3.5a) $\displaystyle\delta_{t}s^{n+1}$
$\displaystyle=-g(u^{n},s^{n})\langle f(u^{n}),\delta_{t}u^{n+1}\rangle,$
(3.5b)

where $\kappa\geq 0$ is a stabilizing constant and $g(u^{n},s^{n})>0.$ The
scheme (3.5) is started by $u^{0}=u_{\text{\rm init}}$ and
$s^{0}=E_{1h}(u^{0})$. We can rewrite (3.5) equivalently as follows:

$\displaystyle\Big{[}\Big{(}\frac{1}{{\tau}}+\kappa
g(u^{n},s^{n})\Big{)}I-\varepsilon^{2}\Delta_{h}\Big{]}u^{n+1}$
$\displaystyle=\frac{u^{n}}{{\tau}}+g(u^{n},s^{n})f(u^{n})+\kappa
g(u^{n},s^{n})u^{n},$ (3.6a) $\displaystyle s^{n+1}$
$\displaystyle=s^{n}-g(u^{n},s^{n})\langle f(u^{n}),u^{n+1}-u^{n}\rangle.$
(3.6b)

Obviously, (3.6) is uniquely solvable for any ${\tau}>0$ since
$(\frac{1}{{\tau}}+\kappa g(u^{n},s^{n}))I-\varepsilon^{2}\Delta_{h}$ is self-
adjoint and positive definite, which makes $u^{n+1}$ linearly determined from
(3.6a) and then $s^{n+1}$ computed explicitly by (3.6b). If we take $\kappa=0$
and $r^{n}=\exp\\{s^{n}\\}$, i.e., $s^{n}=\ln r^{n}$, it is easy to verify
that the scheme (3.5) gives us exactly the ESAV scheme (2.9) with $\kappa=0$.
However, they differ when $\kappa>0$.

####  3.1.1 Energy dissipation and MBP

Now let us define a discrete energy as follows

$\mathcal{E}_{h}(u^{n},s^{n}):=\frac{\varepsilon^{2}}{2}\|\nabla_{h}u^{n}\|^{2}+s^{n},$
(3.7)

which is clearly again an approximation of the original discrete energy
$E_{h}(u^{n})$. We first show that the sESAV1 scheme (3.5) preserves the
energy dissipation law and the MBP unccondtionally. Then, as an application of
both properties, we also prove the uniform boundedness of the variable
coefficient $g(u^{n},s^{n})$.

######  Theorem 3.1 (Energy dissipation of sESAV1)

For any $\kappa\geq 0$ and ${\tau}>0$, the sESAV1 scheme (3.5) is energy
dissipative in the sense that
$\mathcal{E}_{h}(u^{n+1},s^{n+1})\leq\mathcal{E}_{h}(u^{n},s^{n})$.

###### Proof

Taking the inner product with (3.5) by $u^{n+1}-u^{n}$ yields

$\Big{(}\frac{1}{{\tau}}+\kappa
g(u^{n},s^{n})\Big{)}\|u^{n+1}-u^{n}\|^{2}=\varepsilon^{2}\langle\Delta_{h}u^{n+1},u^{n+1}-u^{n}\rangle+g(u^{n},s^{n})\langle
f(u^{n}),u^{n+1}-u^{n}\rangle.$ (3.8)

Combining (3.8), (3.5b), and the identity

$\langle\Delta_{h}u^{n+1},u^{n+1}-u^{n}\rangle=-\frac{1}{2}\|\nabla_{h}u^{n+1}\|^{2}+\frac{1}{2}\|\nabla_{h}u^{n}\|^{2}-\frac{1}{2}\|\nabla_{h}u^{n+1}-\nabla_{h}u^{n}\|^{2},$

we obtain

$\mathcal{E}_{h}(u^{n+1},s^{n+1})-\mathcal{E}_{h}(u^{n},s^{n})=-\Big{(}\frac{1}{{\tau}}+\kappa
g(u^{n},s^{n})\Big{)}\|u^{n+1}-u^{n}\|^{2}-\frac{\varepsilon^{2}}{2}\|\nabla_{h}u^{n+1}-\nabla_{h}u^{n}\|^{2},$

which completes the proof.

Theorem 3.1 implies that the scheme (3.5) is energy dissipative with respect
to the modified energy $\mathcal{E}_{h}(u^{n},s^{n})$ rather than the original
energy $E_{h}(u^{n})$. Note that $s^{n}\not=E_{1h}(u^{n})$ for $n\geq 1$ in
general, and thus $\mathcal{E}_{h}(u^{n},s^{n})\not=E_{h}(u^{n})$.

###### Corollary 1

For any $\kappa\geq 0$ and ${\tau}>0$, it holds $s^{n}\leq E_{h}(u_{\text{\rm
init}})$ for all $n$.

###### Proof

By Theorem 3.1, since $s^{0}=E_{1h}(u_{\text{\rm init}})$, we have

$\frac{\varepsilon^{2}}{2}\|\nabla_{h}u^{n}\|^{2}+s^{n}=\mathcal{E}_{h}(u^{n},s^{n})\leq\mathcal{E}_{h}(u^{n-1},s^{n-1})\leq\cdots\leq\mathcal{E}_{h}(u^{0},s^{0})=E_{h}(u_{\text{\rm
init}}).$

Dropping off the nonnegative term leads to the expected result.

######  Theorem 3.2 (MBP of sESAV1)

If $\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$, the sESAV1 scheme (3.5)
preserves the MBP for $\\{u^{n}\\}$, i.e., the discrete version of (1.4) is
valid as follows:

$\|u_{\text{\rm
init}}\|_{\infty}\leq\beta\quad\Longrightarrow\quad\|u^{n}\|_{\infty}\leq\beta,\quad\forall\,n.$
(3.9)

###### Proof

Suppose $(u^{n},s^{n})$ is given and $\|u^{n}\|_{\infty}\leq\beta$ for some
$n$. From (3.6a), we have

$u^{n+1}=\Big{[}\Big{(}\frac{1}{{\tau}}+\kappa
g(u^{n},s^{n})\Big{)}I-\varepsilon^{2}\Delta_{h}\Big{]}^{-1}\Big{[}\frac{1}{{\tau}}u^{n}+g(u^{n},s^{n})(f(u^{n})+\kappa
u^{n})\Big{]}.$

Since $g(u^{n},s^{n})>0$, by Lemma 1, we have

$\Big{\|}\Big{[}\Big{(}\frac{1}{{\tau}}+\kappa
g(u^{n},s^{n})\Big{)}I-\varepsilon^{2}\Delta_{h}\Big{]}^{-1}\Big{\|}_{\infty}\leq\Big{(}\frac{1}{{\tau}}+\kappa
g(u^{n},s^{n})\Big{)}^{-1}.$

Since $\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$ and
$\|u^{n}\|_{\infty}\leq\beta$, according to Lemma 2, it holds

$\Big{\|}\frac{1}{{\tau}}u^{n}+g(u^{n},s^{n})(f(u^{n})+\kappa
u^{n})\Big{\|}_{\infty}\leq\Big{(}\frac{1}{{\tau}}+\kappa
g(u^{n},s^{n})\Big{)}\beta.$ (3.10)

Therefore, we obtain

$\|u^{n+1}\|_{\infty}\leq\Big{(}\frac{1}{{\tau}}+\kappa
g(u^{n},s^{n})\Big{)}^{-1}\Big{(}\frac{1}{{\tau}}+\kappa
g(u^{n},s^{n})\Big{)}\beta=\beta.$

By induction, we have $\|u^{n}\|_{\infty}\leq\beta$ for all $n$.

###### Remark 1

The inequality (3.10) is valid if
$({\tau}g(u^{n},s^{n}))^{-1}+\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$. In
other words, when $\kappa=0$ (no stabilization), the MBP still holds for the
sESAV1 scheme if the time step size satisfies
${\tau}\leq(g(u^{n},s^{n})\|f^{\prime}\|_{C[-\beta,\beta]})^{-1}$ for all $n$.

###### Remark 2

For the sESAV1 scheme (3.5), we know that the extra term $-\kappa
g(u^{n},s^{n})(u^{n+1}-u^{n})$ stabilizes the time stepping and $\kappa
g(u^{n},s^{n})$ is indeed the stabilizing constant, which is an $n$-dependent
quantity. In the proof of Theorem 3.2, the key ingredients to preserve the MBP
for $\\{u^{n}\\}$ involve two aspects: the positivity of $\kappa
g(u^{n},s^{n})$ and the relation of $u^{n+1}$ and $u^{n}$. The former implies
that the extra term is really a good stabilization term and the latter
guarantees the balance between the linear and nonlinear parts so that the
stabilized linear operator is sufficient to dominate the nonlinear term in
order to preserve the MBP.

###### Remark 3

For the classic SAV1 scheme (2.6), the stabilization term in (2.6a) actually
takes the form

$-\kappa u^{n+1}+\frac{r^{n+1}}{\sqrt{E_{2h}(u^{n})+\delta}}\kappa u^{n}.$

The sign of $r^{n+1}$, and thus the sign of
$\frac{r^{n+1}}{\sqrt{E_{2h}(u^{n})+\delta}}\kappa$, is uncertain, which
violates the positivity of the stabilizing constant. Even though $r^{n+1}$ may
be positive in practical computations in some specific cases, such a
stabilization term leads to an imbalance between the linear and nonlinear
parts since $r^{n+1}\not=\sqrt{E_{2h}(u^{n})+\delta}$ for $n\geq 0$ in
general, so the scheme (2.6) cannot preserve the MBP theoretically, which will
be also observed later in our numerical experiments. Similarly, the
stabilization term in the ESAV scheme (2.9a) reads as

$-\kappa u^{n+1}+\frac{r^{n}}{\exp\\{E_{2h}(u^{n})\\}}\kappa u^{n},$

and the imbalance also exists between the linear and nonlinear parts since
$r^{n}\not=\exp\\{E_{2h}(u^{n})\\}$ for $n\geq 1$ in general, and thus the
ESAV scheme (2.9) also does not preserve the MBP theoretically. Nevertheless,
$r^{n}>0$ always holds due to the definition of the auxiliary variable, and
this is the reason why we consider the ESAV approach rather than the classic
one in this work.

Note that the coefficient $g(u^{n},s^{n})$ may vary step-by-step, which is
different from the continuous case that $g(u,s)\equiv 1$ exactly. Fortunately,
the change of $g(u^{n},s^{n})$ is controllable in the sense that it can be
bounded by some constants, which is illustrated in the following.

###### Corollary 2

If $h\leq h_{0}$, $\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$, and
$\|u_{\text{\rm init}}\|_{\infty}\leq\beta$, then there exists a constant
$G^{*}=G^{*}(u_{\text{\rm init}},C_{*})$ such that $0<g(u^{n},s^{n})\leq
G^{*}$ for all $n$.

###### Proof

The positivity of $g(u^{n},s^{n})$ comes from its definition. We know from
Corollary 1 and Theorem 3.2 that $g(u^{n},s^{n})\leq\exp\\{E_{h}(u_{\text{\rm
init}})+C_{*}\\}$ for all $n$. The $h$-dependence of upper bound can be
removed by (3.1), which completes the proof.

Actually, it also holds that $g(u^{n},s^{n})$ has a positive lower bound
uniformly in $n$ for any fixed terminal time $T>0$. To show it, we first prove
an estimate on the discrete $H^{2}$ semi-norm of the numerical solution.

###### Lemma 3

Given a fixed time $T>0$. If $h\leq h_{0}$,
$\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$, and $\|u_{\text{\rm
init}}\|_{\infty}\leq\beta$, there exists a constant $M>0$ depending on
$C_{*}$, $|\Omega|$, $T$, $u_{\text{\rm init}}$, $\kappa$, $\varepsilon$, and
$\|f\|_{C^{1}[-\beta,\beta]}$, such that

$\|\delta_{t}u^{n+1}\|+\|\Delta_{h}u^{n+1}\|\leq M,\quad 0\leq n\leq\lfloor
T/{\tau}\rfloor-1.$

###### Proof

Taking the discrete inner product of (3.5a) with
$2{\tau}\Delta_{h}^{2}u^{n+1}$, we obtain

$\displaystyle(1+\kappa
g(u^{n},s^{n}){\tau})\langle\Delta_{h}u^{n+1}-\Delta_{h}u^{n},2\Delta_{h}u^{n+1}\rangle+2\varepsilon^{2}{\tau}\|\nabla_{h}\Delta_{h}u^{n+1}\|^{2}$
$\displaystyle\qquad\quad=-2g(u^{n},s^{n}){\tau}\langle\nabla_{h}f(u^{n}),\nabla_{h}\Delta_{h}u^{n+1}\rangle.$

Using the facts that

$\displaystyle\langle\Delta_{h}u^{n+1}-\Delta_{h}u^{n},2\Delta_{h}u^{n+1}\rangle$
$\displaystyle=\|\Delta_{h}u^{n+1}\|^{2}-\|\Delta_{h}u^{n}\|^{2}+\|\Delta_{h}u^{n+1}-\Delta_{h}u^{n}\|^{2},$
$\displaystyle-2g(u^{n},s^{n}){\tau}\langle\nabla_{h}f(u^{n}),\nabla_{h}\Delta_{h}u^{n+1}\rangle$
$\displaystyle\leq\frac{(g(u^{n},s^{n}))^{2}}{2\varepsilon^{2}}{\tau}\|\nabla_{h}f(u^{n})\|^{2}+2\varepsilon^{2}{\tau}\|\nabla_{h}\Delta_{h}u^{n+1}\|^{2},$

we obtain

$(1+\kappa
g(u^{n},s^{n}){\tau})(\|\Delta_{h}u^{n+1}\|^{2}-\|\Delta_{h}u^{n}\|^{2})\leq\frac{(g(u^{n},s^{n}))^{2}}{2\varepsilon^{2}}{\tau}\|\nabla_{h}f(u^{n})\|^{2}.$
(3.11)

By Theorem 3.2, we have $\|u^{n}\|_{\infty}\leq\beta$, and thus,

$\|\nabla_{h}f(u^{n})\|\leq\|f^{\prime}\|_{C[-\beta,\beta]}\|\nabla_{h}u^{n}\|\leq\|f^{\prime}\|_{C[-\beta,\beta]}C_{\Omega}\|\Delta_{h}u^{n}\|,$
(3.12)

where the second step comes from the discrete Poincaré’s inequality with
$C_{\Omega}$ being a constant depending only on $|\Omega|$ (since
$\nabla_{h}u^{n}$ has a zero mean due to the periodic boundary condition).
Then, by Corollary 2, (3.11) and (3.12), we obtain

$\displaystyle\|\Delta_{h}u^{n+1}\|^{2}$ $\displaystyle\leq(1+\kappa
g(u^{n},s^{n}){\tau})\|\Delta_{h}u^{n+1}\|^{2}$
$\displaystyle\leq\Big{[}1+\Big{(}\kappa
G^{*}+\frac{(G^{*}\|f^{\prime}\|_{C[-\beta,\beta]}C_{\Omega})^{2}}{2\varepsilon^{2}}\Big{)}{\tau}\Big{]}\|\Delta_{h}u^{n}\|^{2}.$
(3.13)

By recursion, we obtain

$\displaystyle\|\Delta_{h}u^{n+1}\|^{2}$
$\displaystyle\leq\Big{[}1+\Big{(}\kappa
G^{*}+\frac{(G^{*}\|f^{\prime}\|_{C[-\beta,\beta]}C_{\Omega})^{2}}{2\varepsilon^{2}}\Big{)}{\tau}\Big{]}^{n+1}\|\Delta_{h}u^{0}\|^{2}$
$\displaystyle\leq\mathrm{e}^{\big{(}\kappa
G^{*}+\frac{(G^{*}\|f^{\prime}\|_{C[-\beta,\beta]}C_{\Omega})^{2}}{2\varepsilon^{2}}\big{)}T}\|\Delta_{h}u_{\text{\rm
init}}\|^{2}.$

Then, using Corollary 2 again, we derive from (3.5a) directly to get

$\displaystyle\|\delta_{t}u^{n+1}\|$ $\displaystyle\leq(1+\kappa
g(u^{n},s^{n}){\tau})\|\delta_{t}u^{n+1}\|$
$\displaystyle\leq\varepsilon^{2}\|\Delta_{h}u^{n+1}\|+g(u^{n},s^{n})\|f(u^{n})\|\leq\varepsilon^{2}\|\Delta_{h}u^{n+1}\|+G^{*}F_{0}|\Omega|^{\frac{1}{2}},$

where $F_{0}:=\|f\|_{C[-\beta,\beta]}$. This completes the proof.

###### Corollary 3

Given a fixed time $T>0$. If $h\leq h_{0}$,
$\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$, and $\|u_{\text{\rm
init}}\|_{\infty}\leq\beta$, there exists a constant $G_{*}>0$ such that
$g(u^{n},s^{n})\geq G_{*}$ for $0\leq n\leq\lfloor T/{\tau}\rfloor$, where
$G_{*}$ depends on $C_{*}$, $|\Omega|$, $T$, $u_{\text{\rm init}}$, $\kappa$,
$\varepsilon$, and $\|f\|_{C^{1}[-\beta,\beta]}$.

###### Proof

According to the definition of $g(u^{n},s^{n})$ in (3.4) and the MBP for
$\\{u^{n}\\}$, it suffices to show the existence of the lower bound of
$\\{s^{n}\\}$. Using Lemma 3, we have

$\langle
f(u^{n}),u^{n+1}-u^{n}\rangle\leq{\tau}\|f(u^{n})\|\|\delta_{t}u^{n+1}\|\leq
F_{0}|\Omega|^{\frac{1}{2}}M{\tau},$

where $M$ is the constant defined in Lemma 3. Then, from (3.6b), we have

$s^{n+1}\geq s^{n}-G^{*}F_{0}|\Omega|^{\frac{1}{2}}M{\tau}.$

By recursion, noting that $s^{0}=E_{1h}(u_{\text{\rm init}})\geq-C_{*}$, we
obtain

$s^{n}\geq s^{0}-G^{*}F_{0}|\Omega|^{\frac{1}{2}}Mn{\tau}\geq-
C_{*}-G^{*}F_{0}|\Omega|^{\frac{1}{2}}MT,$

which completes the proof.

The combination of Corollaries 2 and 3 implies that $0<G_{*}\leq
g(u^{n},s^{n})\leq G^{*}$ for any fixed terminal time $T>0$, which will play
an important role in error estimates of the sESAV1 scheme (3.5) in the next
subsection.

####  3.1.2 Error estimates

In the following error analysis, as well as that for the second-order scheme
presented later, we will use many generic constants, and for simplicity of
notations, we may denote the constants with the same dependence but different
values by the same notation.

If the exact solution $u_{e}$ to (1.1) is smooth sufficiently, letting
$s_{e}(t)=E_{1}(u_{e}(t))$, we have

$\displaystyle\frac{u_{e}(t_{n+1})-u_{e}(t_{n})}{{\tau}}$
$\displaystyle=\varepsilon^{2}\Delta_{h}u_{e}(t_{n+1})+g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n}))$
$\displaystyle\quad-\kappa
g(u_{e}(t_{n}),s_{e}(t_{n}))(u_{e}(t_{n+1})-u_{e}(t_{n}))+R_{1u}^{n},$ (3.14a)
$\displaystyle\frac{s_{e}(t_{n+1})-s_{e}(t_{n})}{{\tau}}$
$\displaystyle=-g(u_{e}(t_{n}),s_{e}(t_{n}))\Big{\langle}f(u_{e}(t_{n})),\frac{u_{e}(t_{n+1})-u_{e}(t_{n})}{{\tau}}\Big{\rangle}+R_{1s}^{n},$
(3.14b)

where the truncation errors $R_{1u}^{n}$ and $R_{1s}^{n}$ satisfy

$\|R_{1u}^{n}\|\leq C_{e}({\tau}+h^{2}),\qquad|R_{1s}^{n}|\leq
C_{e}({\tau}+h^{2})$ (3.15)

with $C_{e}>0$ depending only on $u_{e}$, $\kappa$, $\varepsilon$, and
$\|f\|_{C^{1}[-\beta,\beta]}$. Define the error functions as

$e_{u}^{n}=u^{n}-u_{e}(t_{n}),\qquad e_{s}^{n}=s^{n}-s_{e}(t_{n}).$ (3.16)

We first show a lemma on the error estimate for the nonlinear term.

###### Lemma 4

If $h\leq h_{0}$ and $\|u^{n}\|_{\infty}\leq\beta$, we have

$|g(u^{n},s^{n})-g(u_{e}(t_{n}),s_{e}(t_{n}))|\leq
C_{g}(\|e_{u}^{n}\|+|e_{s}^{n}|),$ (3.17)

and

$\|g(u^{n},s^{n})f(u^{n})-g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n}))\|\leq
C_{g}(\|e_{u}^{n}\|+|e_{s}^{n}|),$ (3.18)

where the constant $C_{g}>0$ depends on $C_{*}$, $|\Omega|$, $u_{\text{\rm
init}}$, and $\|f\|_{C^{1}[-\beta,\beta]}$.

###### Proof

For the exact solutions $u_{e}(t_{n})$ and $s_{e}(t_{n})$, we have
$\|u_{e}(t_{n})\|_{\infty}\leq\beta$ by the MBP and $s_{e}(t_{n})\leq
E(u_{\text{\rm init}})$ by the energy dissipation law. Some careful
calculations yield

$\displaystyle|g(u^{n},s^{n})-g(u^{n},s_{e}(t_{n}))|$
$\displaystyle=\frac{1}{\exp\\{E_{1h}(u^{n})\\}}|\exp\\{s^{n}\\}-\exp\\{s_{e}(t_{n})\\}|$
$\displaystyle\leq\frac{\exp\\{\xi^{n}\\}}{\exp\\{E_{1h}(u^{n})\\}}|s^{n}-s_{e}(t_{n})|$
$\displaystyle\leq G^{*}|s^{n}-s_{e}(t_{n})|$

with $\xi^{n}$ being a number between $s^{n}$ and $s_{e}(t_{n})$, and

$\displaystyle|g(u^{n},s_{e}(t_{n}))-g(u_{e}(t_{n}),s_{e}(t_{n}))|$
$\displaystyle\qquad\quad=\exp\\{s_{e}(t_{n})\\}\Big{|}\frac{1}{\exp\\{E_{1h}(u^{n})\\}}-\frac{1}{\exp\\{E_{1h}(u_{e}(t_{n}))\\}}\Big{|}$
$\displaystyle\qquad\quad\leq(\exp\\{E(u_{\text{\rm
init}})+C_{*}\\})|E_{1h}(u^{n})-E_{1h}(u_{e}(t_{n}))|$
$\displaystyle\qquad\quad\leq|\Omega|^{\frac{1}{2}}(\exp\\{E(u_{\text{\rm
init}})+C_{*}\\})\|F(u^{n})-F(u_{e}(t_{n}))\|$ $\displaystyle\qquad\quad\leq
F_{0}|\Omega|^{\frac{1}{2}}(\exp\\{E(u_{\text{\rm
init}})+C_{*}\\})\|u^{n}-u_{e}(t_{n})\|.$

By combining both of the above inqualities, we obtain (3.17). In addition, we
have

$\displaystyle\|g(u^{n},s^{n})f(u_{e}(t_{n}))-g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n}))\|$
$\displaystyle\qquad\quad\leq\|f(u_{e}(t_{n}))\||g(u^{n},s^{n})-g(u_{e}(t_{n}),s_{e}(t_{n}))|$
$\displaystyle\qquad\quad\leq
F_{0}|\Omega|^{\frac{1}{2}}C(\|e_{u}^{n}\|+|e_{s}^{n}|).$

According to Corollary 2, it holds

$\displaystyle\|g(u^{n},s^{n})f(u^{n})-g(u^{n},s^{n})f(u_{e}(t_{n}))\|$
$\displaystyle\leq G^{*}\|f(u^{n})-f(u_{e}(t_{n}))\|$ $\displaystyle\leq
G^{*}\|f^{\prime}\|_{C[-\beta,\beta]}\|u^{n}-u_{e}(t_{n})\|.$

Then, we obtain (3.18) with the help of the triangular inequality to the above
two inequalities.

######  Theorem 3.3 (Error estimate of sESAV1)

Given a fixed time $T>0$ and suppose the exact solution $u_{e}$ is smooth
enough on $[0,T]\times\overline{\Omega}$. Assume that
$\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$ and $\|u_{\text{\rm
init}}\|_{\infty}\leq\beta$. If ${\tau}$ and $h$ are small sufficiently, then
we have the error estimate for the sESAV1 scheme (3.5) as follows:

$\|e_{u}^{n}\|+\|\nabla_{h}e_{u}^{n}\|+|e_{s}^{n}|\leq C({\tau}+h^{2}),\qquad
0\leq n\leq\lfloor T/{\tau}\rfloor,$

where the constant $C>0$ depends on $C_{*}$, $|\Omega|$, $T$, $u_{e}$,
$\kappa$, $\varepsilon$, and $\|f\|_{C^{1}[-\beta,\beta]}$ but is independent
of ${\tau}$ and $h$.

###### Proof

The difference between (3.5) and (3.14) leads to

$\displaystyle\delta_{t}e_{u}^{n+1}$
$\displaystyle=\varepsilon^{2}\Delta_{h}e_{u}^{n+1}+g(u^{n},s^{n})f(u^{n})-g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n}))-\kappa
g(u^{n},s^{n})(e_{u}^{n+1}-e_{u}^{n})$
$\displaystyle\quad+\kappa(g(u_{e}(t_{n}),s_{e}(t_{n}))-g(u^{n},s^{n}))(u_{e}(t_{n+1})-u_{e}(t_{n}))-R_{1u}^{n},$
(3.19a) $\displaystyle\delta_{t}e_{s}^{n+1}$
$\displaystyle=\Big{\langle}g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n}))-g(u^{n},s^{n})f(u^{n}),\frac{u_{e}(t_{n+1})-u_{e}(t_{n})}{{\tau}}\Big{\rangle}$
$\displaystyle\quad-g(u^{n},s^{n})\langle
f(u^{n}),\delta_{t}e_{u}^{n+1}\rangle-R_{1s}^{n}.$ (3.19b)

Taking the discrete inner product of (3.19a) with
$2{\tau}\delta_{t}e_{u}^{n+1}$ and rearranging the terms give us

$\displaystyle 2{\tau}\|\delta_{t}e_{u}^{n+1}\|^{2}+2\kappa
g(u^{n},s^{n})\|e_{u}^{n+1}-e_{u}^{n}\|^{2}-2\varepsilon^{2}\langle\Delta_{h}e_{u}^{n+1},e_{u}^{n+1}-e_{u}^{n}\rangle$
$\displaystyle\qquad=2{\tau}\langle
g(u^{n},s^{n})f(u^{n})-g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n})),\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\qquad\quad+2\kappa{\tau}(g(u_{e}(t_{n}),s_{e}(t_{n}))-g(u^{n},s^{n}))\langle
u_{e}(t_{n+1})-u_{e}(t_{n}),\delta_{t}e_{u}^{n+1}\rangle-2{\tau}\langle
R_{1u}^{n},\delta_{t}e_{u}^{n+1}\rangle.$

Since $g(u^{n},s^{n})\geq G_{*}>0$ by Corollary 3, using the identities

$\displaystyle\langle\Delta_{h}e_{u}^{n+1},e_{u}^{n+1}-e_{u}^{n}\rangle$
$\displaystyle=-\frac{1}{2}\|\nabla_{h}e_{u}^{n+1}\|^{2}+\frac{1}{2}\|\nabla_{h}e_{u}^{n}\|^{2}-\frac{1}{2}{\tau}^{2}\|\nabla_{h}\delta_{t}e_{u}^{n+1}\|^{2},$
$\displaystyle\|e_{u}^{n+1}-e_{u}^{n}\|^{2}$
$\displaystyle=\|e_{u}^{n+1}\|^{2}-\|e_{u}^{n}\|^{2}-2{\tau}\langle
e_{u}^{n},\delta_{t}e_{u}^{n+1}\rangle,$ (3.20)

we obtain

$\displaystyle
2G_{*}\kappa\|e_{u}^{n+1}\|^{2}-2G_{*}\kappa\|e_{u}^{n}\|^{2}+\varepsilon^{2}\|\nabla_{h}e_{u}^{n+1}\|^{2}-\varepsilon^{2}\|\nabla_{h}e_{u}^{n}\|^{2}+2{\tau}\|\delta_{t}e_{u}^{n+1}\|^{2}$
$\displaystyle\qquad\quad\leq 2{\tau}\langle
g(u^{n},s^{n})f(u^{n})-g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n})),\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\qquad\qquad+2\kappa{\tau}(g(u_{e}(t_{n}),s_{e}(t_{n}))-g(u^{n},s^{n}))\langle
u_{e}(t_{n+1})-u_{e}(t_{n}),\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\qquad\qquad+4G_{*}\kappa{\tau}\langle
e_{u}^{n},\delta_{t}e_{u}^{n+1}\rangle-2{\tau}\langle
R_{1u}^{n},\delta_{t}e_{u}^{n+1}\rangle.$ (3.21)

For the first term in the right-hand side of (3.21), by Lemma 4 we have

$\displaystyle 2{\tau}\langle
g(u^{n},s^{n})f(u^{n})-g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n})),\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\qquad\quad\leq
2C_{g}{\tau}(\|e_{u}^{n}\|+|e_{s}^{n}|)\|\delta_{t}e_{u}^{n+1}\|\leq
4C_{g}^{2}{\tau}(\|e_{u}^{n}\|^{2}+|e_{s}^{n}|^{2})+\frac{{\tau}}{2}\|\delta_{t}e_{u}^{n+1}\|^{2},$
(3.22)

where $C_{g}>0$ is the constant in Lemma 4. For the second term in the right-
hand side of (3.21), we have

$\displaystyle
2\kappa{\tau}(g(u_{e}(t_{n}),s_{e}(t_{n}))-g(u^{n},s^{n}))\langle
u_{e}(t_{n+1})-u_{e}(t_{n}),\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\qquad\quad\leq
2\kappa{\tau}|g(u_{e}(t_{n}),s_{e}(t_{n}))-g(u^{n},s^{n})|\|u_{e}(t_{n+1})-u_{e}(t_{n})\|\|\delta_{t}e_{u}^{n+1}\|$
$\displaystyle\qquad\quad\leq
2C_{g}\kappa{\tau}(\|u_{e}(t_{n+1})\|+\|u_{e}(t_{n})\|)(\|e_{u}^{n}\|+|e_{s}^{n}|)\|\delta_{t}e_{u}^{n+1}\|$
$\displaystyle\qquad\quad\leq
C_{1}\kappa^{2}{\tau}(\|e_{u}^{n}\|^{2}+|e_{s}^{n}|^{2})+\frac{{\tau}}{2}\|\delta_{t}e_{u}^{n+1}\|^{2},$
(3.23)

where $C_{1}>0$ depends on $C_{*}$, $|\Omega|$, $u_{e}$, and
$\|f\|_{C^{1}[-\beta,\beta]}$. Using the Young’s inequality, the third and
fourth terms in the right-hand side of (3.21) can be bounded respectively as

$\displaystyle 4G_{*}\kappa{\tau}\langle
e_{u}^{n},\delta_{t}e_{u}^{n+1}\rangle\leq
4G_{*}\kappa{\tau}\|e_{u}^{n}\|\|\delta_{t}e_{u}^{n+1}\|\leq
8G_{*}^{2}\kappa^{2}{\tau}\|e_{u}^{n}\|^{2}+\dfrac{{\tau}}{2}\|\delta_{t}e_{u}^{n+1}\|^{2},\qquad$
(3.24) $\displaystyle-2{\tau}\langle
R_{1u}^{n},\delta_{t}e_{u}^{n+1}\rangle\leq
2{\tau}\|R_{1u}^{n}\|\|\delta_{t}e_{u}^{n+1}\|\leq
4{\tau}\|R_{1u}^{n}\|^{2}+\dfrac{{\tau}}{4}\|\delta_{t}e_{u}^{n+1}\|^{2}.\quad$
(3.25)

Then, substituting (3.22)–(3.25) into (3.21) leads to

$\displaystyle
2G_{*}\kappa\|e_{u}^{n+1}\|^{2}-2G_{*}\kappa\|e_{u}^{n}\|^{2}+\varepsilon^{2}\|\nabla_{h}e_{u}^{n+1}\|^{2}-\varepsilon^{2}\|\nabla_{h}e_{u}^{n}\|^{2}+\frac{{\tau}}{4}\|\delta_{t}e_{u}^{n+1}\|^{2}$
$\displaystyle\qquad\leq(4C_{g}^{2}+C_{1}\kappa^{2}+8G_{*}^{2}\kappa^{2}){\tau}\|e_{u}^{n}\|^{2}+(4C_{g}^{2}+C_{1}\kappa^{2}){\tau}|e_{s}^{n}|^{2}+4{\tau}\|R_{1u}^{n}\|^{2}.$
(3.26)

Multiplying (3.19b) by $2{\tau}e_{s}^{n+1}$ yields

$\displaystyle|e_{s}^{n+1}|^{2}-|e_{s}^{n}|^{2}+|e_{s}^{n+1}-e_{s}^{n}|^{2}$
$\displaystyle\qquad\quad=2e_{s}^{n+1}\langle
g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n}))-g(u^{n},s^{n})f(u^{n}),u_{e}(t_{n+1})-u_{e}(t_{n})\rangle$
$\displaystyle\qquad\qquad-2{\tau}e_{s}^{n+1}g(u^{n},s^{n})\langle
f(u^{n}),\delta_{t}e_{u}^{n+1}\rangle-2{\tau}R_{1s}^{n}e_{s}^{n+1}.$ (3.27)

For the first term in the right-hand side of (3.27), by Lemma 4 we have

$\displaystyle 2e_{s}^{n+1}\langle
g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n}))-g(u^{n},s^{n})f(u^{n}),u_{e}(t_{n+1})-u_{e}(t_{n})\rangle$
$\displaystyle\qquad\quad\leq
2|e_{s}^{n+1}|\|g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n}))-g(u^{n},s^{n})f(u^{n})\|\|u_{e}(t_{n+1})-u_{e}(t_{n})\|$
$\displaystyle\qquad\quad\leq
2C_{g}{\tau}|e_{s}^{n+1}|(\|e_{u}^{n}\|+|e_{s}^{n}|)\|(u_{e})_{t}(\theta_{n})\|\qquad(\mbox{for
some }t_{n}<\theta_{n}<t_{n+1})$ $\displaystyle\qquad\quad\leq
C_{2}{\tau}(\|e_{u}^{n}\|^{2}+|e_{s}^{n}|^{2}+|e_{s}^{n+1}|^{2}),$ (3.28)

where $C_{2}>0$ depends on $C_{*}$, $|\Omega|$, $u_{e}$, and
$\|f\|_{C^{1}[-\beta,\beta]}$. For the second term in the right-hand side of
(3.27), using Corollary 2, we obtain

$\displaystyle-2{\tau}e_{s}^{n+1}g(u^{n},s^{n})\langle
f(u^{n}),\delta_{t}e_{u}^{n+1}\rangle$ $\displaystyle\leq
2G^{*}{\tau}\|f(u^{n})\||e_{s}^{n+1}|\|\delta_{t}e_{u}^{n+1}\|$
$\displaystyle\leq
C_{3}{\tau}|e_{s}^{n+1}|^{2}+\frac{{\tau}}{4}\|\delta_{t}e_{u}^{n+1}\|^{2},$
(3.29)

where $C_{3}>0$ depends on $C_{*}$, $|\Omega|$, $u_{\text{\rm init}}$, and
$\|f\|_{C[-\beta,\beta]}$. For the third term in the right-hand side of
(3.27), we have

$-2{\tau}R_{1s}^{n}e_{s}^{n+1}\leq{\tau}|R_{1s}^{n}|^{2}+{\tau}|e_{s}^{n+1}|^{2}.$
(3.30)

Then, substituting (3.28)–(3.30) into (3.27) leads to

$|e_{s}^{n+1}|^{2}-|e_{s}^{n}|^{2}\leq
C_{2}{\tau}\|e_{u}^{n}\|^{2}+C_{2}{\tau}|e_{s}^{n}|^{2}+(1+C_{2}+C_{3}){\tau}|e_{s}^{n+1}|^{2}+\frac{{\tau}}{4}\|\delta_{t}e_{u}^{n+1}\|^{2}+{\tau}|R_{1s}^{n}|^{2}.$
(3.31)

Adding (3.26) and (3.31), we obtain

$\displaystyle
2G_{*}\kappa(\|e_{u}^{n+1}\|^{2}-\|e_{u}^{n}\|^{2})+\varepsilon^{2}(\|\nabla_{h}e_{u}^{n+1}\|^{2}-\|\nabla_{h}e_{u}^{n}\|^{2})+(|e_{s}^{n+1}|^{2}-|e_{s}^{n}|^{2})$
$\displaystyle\qquad\quad\leq(4C_{g}^{2}+C_{1}\kappa^{2}+8G_{*}^{2}\kappa^{2}+C_{2}){\tau}\|e_{u}^{n}\|^{2}+(4C_{g}^{2}+C_{1}\kappa^{2}+C_{2}){\tau}|e_{s}^{n}|^{2}$
$\displaystyle\qquad\qquad+(1+C_{2}+C_{3}){\tau}|e_{s}^{n+1}|^{2}+4{\tau}\|R_{1u}^{n}\|^{2}+{\tau}|R_{1s}^{n}|^{2}.$

Then, using (3.15), we reach

$\displaystyle\quad~{}2G_{*}\kappa(\|e_{u}^{n+1}\|^{2}-\|e_{u}^{n}\|^{2})+\varepsilon^{2}(\|\nabla_{h}e_{u}^{n+1}\|^{2}-\|\nabla_{h}e_{u}^{n}\|^{2})+(|e_{s}^{n+1}|^{2}-|e_{s}^{n}|^{2})$
$\displaystyle\qquad\qquad\leq
C_{4}{\tau}(\|e_{u}^{n}\|^{2}+|e_{s}^{n}|^{2}+|e_{s}^{n+1}|^{2})+5C_{e}^{2}{\tau}({\tau}+h^{2})^{2},$

where the constant $C_{4}$ depends on $C_{*}$, $|\Omega|$, $T$, $u_{e}$,
$\kappa$, $\varepsilon$, and $\|f\|_{C^{1}[-\beta,\beta]}$.

Letting
$W^{n}:=2G_{*}\kappa\|e_{u}^{n}\|^{2}+\varepsilon^{2}\|\nabla_{h}e_{u}^{n}\|^{2}+|e_{s}^{n}|^{2}$,
we have

$W^{n+1}-W^{n}\leq\widetilde{C}_{4}{\tau}(W^{n}+W^{n+1})+5C_{e}^{2}{\tau}({\tau}+h^{2})^{2},$

where $\widetilde{C}_{4}$ depends on $C_{4}$ and $\kappa$. When
${\tau}\leq\frac{1}{2\widetilde{C}_{4}}$, noting that
$\frac{1+\widetilde{C}_{4}{\tau}}{1-\widetilde{C}_{4}{\tau}}\leq
1+4\widetilde{C}_{4}{\tau}$, we obtain

$W^{n+1}\leq(1+4\widetilde{C}_{4}{\tau})W^{n}+10C_{e}^{2}{\tau}({\tau}+h^{2})^{2}.$

Using the discrete Gronwall’s inequality, we obtain

$2G_{*}\kappa\|e_{u}^{n}\|^{2}+\varepsilon^{2}\|\nabla_{h}e_{u}^{n}\|^{2}+|e_{s}^{n}|^{2}=W^{n}\leq
10C_{e}^{2}\mathrm{e}^{4\widetilde{C}_{4}T}({\tau}+h^{2})^{2},$

which completes the proof.

###### Remark 4

For any fixed $h>0$, let us recall the space-discrete problem (3.3) and denote
by $u_{h,e}(t)$ the exact solution. By similar analysis as Theorem 3.3, one
can obtain the error estimates for sufficiently small ${\tau}$ as follows:

$\|u^{n}-u_{h,e}(t_{n})\|+\|\nabla_{h}u^{n}-\nabla_{h}u_{h,e}(t_{n})\|+|s^{n}-E_{1h}(u_{h,e}(t))|\leq
C_{h}{\tau},$

where the constant $C_{h}>0$ depends on $C_{*}$, $|\Omega|$, $T$, $u_{h,e}$,
$\kappa$, $\varepsilon$, and $\|f\|_{C^{1}[-\beta,\beta]}$ but is independent
of ${\tau}$.

###  3.2 Second-order sESAV scheme

For the space-discrete system (3.3), the second-order stabilized ESAV scheme
(sESAV2) is given by

$\displaystyle\delta_{t}u^{n+1}$
$\displaystyle=\varepsilon^{2}\Delta_{h}u^{n+\frac{1}{2}}+g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})f(\widehat{u}^{n+\frac{1}{2}})-\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})(u^{n+\frac{1}{2}}-\widehat{u}^{n+\frac{1}{2}}),$
(3.32a) $\displaystyle\delta_{t}s^{n+1}$
$\displaystyle=-g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
f(\widehat{u}^{n+\frac{1}{2}}),\delta_{t}u^{n+1}\rangle+\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
u^{n+\frac{1}{2}}-\widehat{u}^{n+\frac{1}{2}},\delta_{t}u^{n+1}\rangle,$
(3.32b)

where $g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})>0$ with
$(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})$ being generated by
the first-order scheme (3.5) with the time step size ${\tau}/2$, i.e.,

$\displaystyle\frac{\widehat{u}^{n+\frac{1}{2}}-u^{n}}{{\tau}/2}$
$\displaystyle=\varepsilon^{2}\Delta_{h}\widehat{u}^{n+\frac{1}{2}}+g(u^{n},s^{n})f(u^{n})-\kappa
g(u^{n},s^{n})(\widehat{u}^{n+\frac{1}{2}}-u^{n}),$ (3.33a)
$\displaystyle\widehat{s}^{n+\frac{1}{2}}-s^{n}$
$\displaystyle=-g(u^{n},s^{n})\langle
f(u^{n}),\widehat{u}^{n+\frac{1}{2}}-u^{n}\rangle.$ (3.33b)

The scheme (3.32) is started by $u^{0}=u_{\text{\rm init}}$ and
$s^{0}=E_{1h}(u^{0})$. By the definition of $u^{n+\frac{1}{2}}$, the last term
in (3.32a) is actually $-\frac{1}{2}\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})(u^{n+1}-2\widehat{u}^{n+\frac{1}{2}}+u^{n})$,
which provides a second-order truncation error in time. We can rewrite (3.32)
in the following form:

$\displaystyle\Big{[}\Big{(}\frac{2}{{\tau}}+\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\Big{)}I-\varepsilon^{2}\Delta_{h}\Big{]}u^{n+1}$
$\displaystyle\qquad=\Big{[}\Big{(}\frac{2}{{\tau}}-\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\Big{)}I+\varepsilon^{2}\Delta_{h}\Big{]}u^{n}+2g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})[f(\widehat{u}^{n+\frac{1}{2}})+\kappa\widehat{u}^{n+\frac{1}{2}}],$
(3.34a) $\displaystyle
s^{n+1}=s^{n}-g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
f(\widehat{u}^{n+\frac{1}{2}})-\kappa(u^{n+\frac{1}{2}}-\widehat{u}^{n+\frac{1}{2}}),u^{n+1}-u^{n}\rangle.$
(3.34b)

It is then easy to see that the system (3.34) is linear and uniquely solvable
for any ${\tau}>0$ since $(\frac{2}{{\tau}}+\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}}))I-\varepsilon^{2}\Delta_{h}$
is self-adjoint and positive definite.

####  3.2.1 Energy dissipation and MBP

The energy dissipation law and the MBP preservation of the sESAV2 scheme
(3.32) are stated below.

######  Theorem 3.4 (Energy dissipation of sESAV2)

For any $\kappa\geq 0$ and ${\tau}>0$, the sESAV2 scheme (3.32) is energy
dissipative in the sense that
$\mathcal{E}_{h}(u^{n+1},s^{n+1})\leq\mathcal{E}_{h}(u^{n},s^{n})$, where
$\mathcal{E}_{h}(u^{n},s^{n})$ is given by (3.7). Moreover, it holds that
$s^{n}\leq E_{h}(u_{\text{\rm init}})$ and $\widehat{s}^{n+\frac{1}{2}}\leq
E_{h}(u_{\text{\rm init}})$ for all $n$.

###### Proof

Taking the inner product of (3.32) with $u^{n+1}-u^{n}$ yields

$\displaystyle\frac{1}{{\tau}}\|u^{n+1}-u^{n}\|^{2}$
$\displaystyle=-\frac{\varepsilon^{2}}{2}\|\nabla_{h}u^{n+1}\|^{2}+\frac{\varepsilon^{2}}{2}\|\nabla_{h}u^{n}\|^{2}+g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
f(\widehat{u}^{n+\frac{1}{2}}),u^{n+1}-u^{n}\rangle$
$\displaystyle\quad-\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
u^{n+\frac{1}{2}}-\widehat{u}^{n+\frac{1}{2}},u^{n+1}-u^{n}\rangle.$ (3.35)

Combining (3.2.1) and (3.32b), we obtain

$\mathcal{E}_{h}(u^{n+1},s^{n+1})-\mathcal{E}_{h}(u^{n},s^{n})=-\frac{1}{{\tau}}\|u^{n+1}-u^{n}\|^{2}\leq
0.$

Similar to the proof of Corollary 1, the uniform boundedness of $\\{s^{n}\\}$
is a direct result of the energy stability. Since
$\widehat{s}^{n+\frac{1}{2}}$ is generated by the sESAV1 scheme (3.33),
according to Theorem 3.1, we have
$\widehat{s}^{n+\frac{1}{2}}\leq\mathcal{E}_{h}(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\leq\mathcal{E}_{h}(u^{n},s^{n})$,
and thus, we have $\widehat{s}^{n+\frac{1}{2}}\leq E_{h}(u_{\text{\rm
init}})$.

######  Theorem 3.5 (MBP of sESAV2)

If $h\leq h_{0}$, $\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$, and

${\tau}\leq\Big{(}\frac{\kappa
G^{*}}{2}+\frac{\varepsilon^{2}}{h^{2}}\Big{)}^{-1},$ (3.36)

where $G^{*}$ is the positive constant defined in Corollary 2, then the sESAV2
scheme (3.32) preserves the MBP for $\\{u^{n}\\}$, i.e., (3.9) is valid.

###### Proof

Suppose $(u^{n},s^{n})$ is given and $\|u^{n}\|_{\infty}\leq\beta$ for some
$n$. By Theorems 3.2 and 3.4, we have
$\|\widehat{u}^{n+\frac{1}{2}}\|_{\infty}\leq\beta$ and
$\widehat{s}^{n+\frac{1}{2}}\leq E_{h}(u_{\text{\rm init}})$. Then, we know
that $0<g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\leq G^{*}$
by the similar analysis as Corollary 2. The condition (3.36) implies

$\frac{2}{{\tau}}-\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\geq\frac{2\varepsilon^{2}}{h^{2}}.$

According to the definition of the matrix $\infty$-norm, we have

$\Big{\|}\Big{(}\frac{2}{{\tau}}-\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\Big{)}I+\varepsilon^{2}\Delta_{h}\Big{\|}_{\infty}=\frac{2}{{\tau}}-\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}}).$

Since $\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$ and
$\|\widehat{u}^{n+\frac{1}{2}}\|_{\infty}\leq\beta$, according to Lemma 2, we
have

$\|f(\widehat{u}^{n+\frac{1}{2}})+\kappa\widehat{u}^{n+\frac{1}{2}}\|_{\infty}\leq\kappa\beta.$

Therefore, using Lemma 1, we obtain from (3.34a) that

$\|u^{n+1}\|_{\infty}\leq\Big{(}\frac{2}{{\tau}}+\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\Big{)}^{-1}\Big{[}\Big{(}\frac{2}{{\tau}}-\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\Big{)}\beta+2\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\beta\Big{]}=\beta.$

By induction, we have $\|u^{n}\|_{\infty}\leq\beta$ for all $n$.

###### Remark 5

Theorem 3.5 implies that $0<g(u^{n},s^{n})\leq G^{*}$ and
$0<g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\leq G^{*}$ hold
for all $n$.

###### Remark 6

The condition (3.36) on the time step size implies
${\tau}=\mathcal{O}(h^{2}/\varepsilon^{2})$, which is the same as those
enforced in HoLe20 ; HoTaYa17 . This restriction comes essentially from the
explicit term $\Delta_{h}u^{n}$ due to the use of the Crank–Nicolson
approximation, which also means that the second-order scheme (3.32) cannot
preserve the MBP unconditionally even though we introduce the stabilization
term. In practical computations,
$g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\approx 1$ so that
the requirement for the time step size can be set to be
${\tau}\leq(\frac{\kappa}{2}+\frac{\varepsilon^{2}}{h^{2}})^{-1}$ in order to
preserve the MBP, which is later used in our numerical experiments.

Similar to the analysis for the sESAV1 scheme, we can show that both
$\\{g(u^{n},s^{n})\\}$ and
$\\{g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\\}$ have
uniform positive lower bounds.

###### Lemma 5

Given a fixed time $T>0$. If $h\leq h_{0}$,
$\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$, $\|u_{\text{\rm
init}}\|_{\infty}\leq\beta$, and ${\tau}\leq 1$ satisfying (3.36), there
exists a constant $M>0$ depending on $C_{*}$, $|\Omega|$, $T$, $u_{\text{\rm
init}}$, $\kappa$, $\varepsilon$, and $\|f\|_{C^{1}[-\beta,\beta]}$ such that

$\displaystyle{\tau}^{-1}\|\widehat{u}^{n+\frac{1}{2}}-u^{n}\|+\|\Delta_{h}\widehat{u}^{n+\frac{1}{2}}\|\leq
M,$ $\displaystyle{\tau}^{-1}\|u^{n+1}-u^{n}\|+\|\Delta_{h}u^{n+1}\|\leq M,$

for $0\leq n\leq\lfloor T/{\tau}\rfloor-1$.

###### Proof

Since $\widehat{u}^{n+\frac{1}{2}}$ is the solution to the sESAV1 substep
(3.33), according to (3.1.1), we have

$\|\Delta_{h}\widehat{u}^{n+\frac{1}{2}}\|^{2}\leq\Big{(}1+\frac{G^{*}\kappa}{2}+\frac{(G^{*}\|f^{\prime}\|_{C[-\beta,\beta]}C_{\Omega})^{2}}{4\varepsilon^{2}}\Big{)}\|\Delta_{h}u^{n}\|^{2},$
(3.37)

where we used ${\tau}\leq 1$.

Taking the discrete inner product of (3.32a) with
$2{\tau}\Delta_{h}^{2}u^{n+\frac{1}{2}}$, using the fact
$0<g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\leq G^{*}$, and
conducting the similar analysis as the proof of Lemma 3, we can obtain

$\displaystyle\|\Delta_{h}u^{n+1}\|^{2}\leq\|\Delta_{h}u^{n}\|^{2}+\Big{(}G^{*}\kappa+\frac{(G^{*}\|f^{\prime}\|_{C[-\beta,\beta]}C_{\Omega})^{2}}{2\varepsilon^{2}}\Big{)}{\tau}\|\Delta_{h}\widehat{u}^{n+\frac{1}{2}}\|^{2}.$

Substituting (3.37) into the above inequality, we have

$\displaystyle\|\Delta_{h}u^{n+1}\|^{2}$
$\displaystyle\leq\Big{[}1+\Big{(}G^{*}\kappa+\frac{(G^{*}\|f^{\prime}\|_{C[-\beta,\beta]}C_{\Omega})^{2}}{2\varepsilon^{2}}\Big{)}$
$\displaystyle\qquad\cdot\Big{(}1+\frac{G^{*}\kappa}{2}+\frac{(C_{\Omega}\|f^{\prime}\|_{C[-\beta,\beta]}G^{*})^{2}}{4\varepsilon^{2}}\Big{)}{\tau}\Big{]}\|\Delta_{h}u^{n}\|^{2}.$

By recursion, we can obtain a uniform upper bound for $\|\Delta_{h}u^{n+1}\|$.
Then, by (3.37) we also can get the upper bound for
$\|\Delta_{h}\widehat{u}^{n+\frac{1}{2}}\|$.

Finally, as a consequence of the above analysis, using (3.33a) and (3.32a), we
also get the boundedness of ${\tau}^{-1}\|\widehat{u}^{n+\frac{1}{2}}-u^{n}\|$
and ${\tau}^{-1}\|u^{n+1}-u^{n}\|$, and the proof is completed.

###### Corollary 4

Given a fixed time $T>0$. If $h\leq h_{0}$,
$\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$, $\|u_{\text{\rm
init}}\|_{\infty}\leq\beta$, and ${\tau}\leq 1$ satisfying (3.36), there
exists a constant ${\widetilde{G}}_{*}>0$ such that
$g(u^{n},s^{n})\geq{\widetilde{G}}_{*}$ and
$g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\geq{\widetilde{G}}_{*}$,
where ${\widetilde{G}}_{*}$ depends on $C_{*}$, $|\Omega|$, $T$, $u_{\text{\rm
init}}$, $\kappa$, $\varepsilon$, and $\|f\|_{C^{1}[-\beta,\beta]}$.

###### Proof

It suffices to show the existence of the lower bounds of $\\{s^{n}\\}$ and
$\\{\widehat{s}^{n+\frac{1}{2}}\\}$. By Lemma 5, we have
$\|u^{n+1}-u^{n}\|\leq M{\tau}$. From (3.32b), the similar analysis as
Corollary 3 leads to the lower boundedness of $\\{s^{n}\\}$. Then, since
$\|\widehat{u}^{n+\frac{1}{2}}-u^{n}\|\leq M{\tau}\leq M$, we can derive from
(3.33b) to give

$\widehat{s}^{n+\frac{1}{2}}\geq s^{n}-G^{*}F_{0}|\Omega|^{\frac{1}{2}}M,$

which completes the proof.

####  3.2.2 Error estimates

It is easy to check that the exact solution $u_{e}$ to (1.1) with
$s_{e}(t)=E_{1}(u_{e}(t))$ satisfies

$\displaystyle\frac{u_{e}(t_{n+1})-u_{e}(t_{n})}{{\tau}}=\frac{\varepsilon^{2}}{2}\Delta_{h}(u_{e}(t_{n+1})+u_{e}(t_{n}))+g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))f(u_{e}(t_{n+\frac{1}{2}}))$
$\displaystyle\qquad-\kappa
g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))\Big{(}\frac{u_{e}(t_{n+1})+u_{e}(t_{n})}{2}-u_{e}(t_{n+\frac{1}{2}})\Big{)}+R_{2u}^{n},$
(3.38a)
$\displaystyle\frac{s_{e}(t_{n+1})-s_{e}(t_{n})}{{\tau}}=-g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))\Big{\langle}f(u_{e}(t_{n+\frac{1}{2}})),\frac{u_{e}(t_{n+1})-u_{e}(t_{n})}{{\tau}}\Big{\rangle}$
$\displaystyle\qquad+\kappa
g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))\Big{\langle}\frac{u_{e}(t_{n+1})+u_{e}(t_{n})}{2}-u_{e}(t_{n+\frac{1}{2}}),\frac{u_{e}(t_{n+1})-u_{e}(t_{n})}{{\tau}}\Big{\rangle}+R_{2s}^{n},$
(3.38b)

where the truncation errors $R_{2u}^{n}$ and $R_{2s}^{n}$ satisfy

$\|R_{2u}^{n}\|\leq C_{e}({\tau}^{2}+h^{2}),\qquad|R_{2s}^{n}|\leq
C_{e}({\tau}^{2}+h^{2}).$ (3.39)

Apart from the numerical error functions $e_{u}^{n}$ and $e_{s}^{n}$ defined
by (3.16), let us also define

$\widehat{e}_{u}^{n+\frac{1}{2}}=\widehat{u}^{n+\frac{1}{2}}-u_{e}(t_{n+\frac{1}{2}}),\qquad\widehat{e}_{s}^{n+\frac{1}{2}}=\widehat{s}^{n+\frac{1}{2}}-s_{e}(t_{n+\frac{1}{2}}).$

We first present an estimate for $\widehat{e}_{u}^{n+\frac{1}{2}}$ and
$\widehat{e}_{s}^{n+\frac{1}{2}}$, which will be used in the proof of the
error estimate for the sESAV2 scheme (3.32). Recalling the proof of Theorem
3.3 for the sESAV1 scheme, the error equations with respect to
$\widehat{e}_{u}^{n+\frac{1}{2}}$ and $\widehat{e}_{s}^{n+\frac{1}{2}}$ read
as

$\displaystyle\frac{\widehat{e}_{u}^{n+\frac{1}{2}}-e_{u}^{n}}{{\tau}/2}$
$\displaystyle=\varepsilon^{2}\Delta_{h}\widehat{e}_{u}^{n+\frac{1}{2}}+g(u^{n},s^{n})f(u^{n})-g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n}))-\kappa
g(u^{n},s^{n})(\widehat{e}_{u}^{n+\frac{1}{2}}-e_{u}^{n})$
$\displaystyle\quad+\kappa(g(u_{e}(t_{n}),s_{e}(t_{n}))-g(u^{n},s^{n}))(u_{e}(t_{n+\frac{1}{2}})-u_{e}(t_{n}))-\widehat{R}_{1u}^{n},$
(3.40a) $\displaystyle\widehat{e}_{s}^{n+\frac{1}{2}}-e_{s}^{n}$
$\displaystyle=\langle
g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n}))-g(u^{n},s^{n})f(u^{n}),u_{e}(t_{n+\frac{1}{2}})-u_{e}(t_{n})\rangle$
$\displaystyle\quad-g(u^{n},s^{n})\langle
f(u^{n}),\widehat{e}_{u}^{n+\frac{1}{2}}-e_{u}^{n}\rangle-\frac{{\tau}}{2}\widehat{R}_{1s}^{n},$
(3.40b)

where

$\|\widehat{R}_{1u}^{n}\|\leq
C_{e}({\tau}+h^{2}),\qquad|\widehat{R}_{1s}^{n}|\leq C_{e}({\tau}+h^{2}).$
(3.41)

###### Lemma 6

Suppose that $h\leq h_{0}$ and $\|u^{n}\|_{\infty}\leq\beta$. If ${\tau}$ is
small sufficiently, we have

$\|\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2}+|\widehat{e}_{s}^{n+\frac{1}{2}}|^{2}\leq\widehat{C}(\|e_{u}^{n}\|^{2}+|e_{s}^{n}|^{2})+\widehat{C}C_{e}^{2}({\tau}^{2}+{\tau}h^{2})^{2},$
(3.42)

where the constant $\widehat{C}>0$ depends on $C_{*}$, $|\Omega|$, $u_{e}$,
$\kappa$, and $\|f\|_{C^{1}[-\beta,\beta]}$.

###### Proof

Taking the discrete inner product of (3.40a) with
${\tau}\widehat{e}_{u}^{n+\frac{1}{2}}$ and rearranging the terms, we have

$\displaystyle\Big{(}1+\frac{\kappa
g(u^{n},s^{n})}{2}{\tau}\Big{)}\big{(}\|\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2}-\|e_{u}^{n}\|^{2}+\|\widehat{e}_{u}^{n+\frac{1}{2}}-e_{u}^{n}\|^{2}\big{)}+\varepsilon^{2}{\tau}\|\nabla_{h}\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2}$
$\displaystyle\qquad={\tau}\langle
g(u^{n},s^{n})f(u^{n})-g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n})),\widehat{e}_{u}^{n+\frac{1}{2}}\rangle$
$\displaystyle\qquad\quad+\kappa{\tau}\langle(g(u_{e}(t_{n}),s_{e}(t_{n}))-g(u^{n},s^{n}))(u_{e}(t_{n+\frac{1}{2}})-u_{e}(t_{n})),\widehat{e}_{u}^{n+\frac{1}{2}}\rangle-{\tau}\langle\widehat{R}_{1u}^{n},\widehat{e}_{u}^{n+\frac{1}{2}}\rangle.$

Using Young’s inequality, we then get

$-{\tau}\langle\widehat{R}_{1u}^{n},\widehat{e}_{u}^{n+\frac{1}{2}}\rangle\leq{\tau}\|\widehat{R}_{1u}^{n}\|\|\widehat{e}_{u}^{n+\frac{1}{2}}\|\leq{\tau}^{2}\|\widehat{R}_{1u}^{n}\|^{2}+\frac{1}{4}\|\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2}.$

Similar to the deductions of (3.22) and (3.23), applying Lemma 4 leads to

$\displaystyle{\tau}\langle
g(u^{n},s^{n})f(u^{n})-g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n})),\widehat{e}_{u}^{n+\frac{1}{2}}\rangle\leq
2C_{g}^{2}{\tau}^{2}(\|e_{u}^{n}\|^{2}+|e_{s}^{n}|^{2})+\frac{1}{4}\|\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2},$
$\displaystyle\kappa{\tau}\langle(g(u_{e}(t_{n}),s_{e}(t_{n}))-g(u^{n},s^{n}))(u_{e}(t_{n+\frac{1}{2}})-u_{e}(t_{n})),\widehat{e}_{u}^{n+\frac{1}{2}}\rangle\leq
C_{1}\kappa^{2}{\tau}^{2}(\|e_{u}^{n}\|^{2}+|e_{s}^{n}|^{2})+\frac{1}{4}\|\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2}.$

Then, we have

$\displaystyle\Big{(}\frac{1}{4}+\frac{\kappa
g(u^{n},s^{n})}{2}{\tau}\Big{)}\|\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2}+\Big{(}1+\frac{\kappa
g(u^{n},s^{n})}{2}{\tau}\Big{)}\|\widehat{e}_{u}^{n+\frac{1}{2}}-e_{u}^{n}\|^{2}$
$\displaystyle\qquad\quad\leq\Big{(}1+\frac{\kappa
g(u^{n},s^{n})}{2}{\tau}\Big{)}\|e_{u}^{n}\|^{2}+(2C_{g}^{2}+C_{1}\kappa^{2}){\tau}^{2}(\|e_{u}^{n}\|^{2}+|e_{s}^{n}|^{2})+{\tau}^{2}\|\widehat{R}_{1u}^{n}\|^{2}.$

By Remark 5, we then can simplify the above equation to get

$\displaystyle\|\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2}+4\|\widehat{e}_{u}^{n+\frac{1}{2}}-e_{u}^{n}\|^{2}$
$\displaystyle\leq(4+2G^{*}\kappa{\tau})\|e_{u}^{n}\|^{2}+(8C_{g}^{2}+4C_{1}\kappa^{2}){\tau}^{2}(\|e_{u}^{n}\|^{2}+|e_{s}^{n}|^{2})+4{\tau}^{2}\|\widehat{R}_{1u}^{n}\|^{2}.$

When ${\tau}\leq 1$, using (3.41), we obtain

$\displaystyle\|\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2}+4\|\widehat{e}_{u}^{n+\frac{1}{2}}-e_{u}^{n}\|^{2}$
$\displaystyle\leq(4+8C_{g}^{2}+2G^{*}\kappa+4C_{1}\kappa^{2})\|e_{u}^{n}\|^{2}$
$\displaystyle\quad+(8C_{g}^{2}+4C_{1}\kappa^{2})|e_{s}^{n}|^{2}+4C_{e}^{2}{\tau}^{2}({\tau}+h^{2})^{2}.$
(3.43)

Multiplying (3.40b) by $2\widehat{e}_{s}^{n+\frac{1}{2}}$ yields

$\displaystyle|\widehat{e}_{s}^{n+\frac{1}{2}}|^{2}-|e_{s}^{n}|^{2}+|\widehat{e}_{s}^{n+\frac{1}{2}}-e_{s}^{n}|^{2}$
$\displaystyle\qquad\quad=2\widehat{e}_{s}^{n+\frac{1}{2}}\langle
g(u_{e}(t_{n}),s_{e}(t_{n}))f(u_{e}(t_{n}))-g(u^{n},s^{n})f(u^{n}),u_{e}(t_{n+\frac{1}{2}})-u_{e}(t_{n})\rangle$
$\displaystyle\qquad\qquad-2\widehat{e}_{s}^{n+\frac{1}{2}}g(u^{n},s^{n})\langle
f(u^{n}),\widehat{e}_{u}^{n+\frac{1}{2}}-e_{u}^{n}\rangle-{\tau}\widehat{R}_{1s}^{n}\widehat{e}_{s}^{n+\frac{1}{2}}.$

We can bound the third and second terms in the right-hand side of the above
equation repectively as

$\displaystyle\displaystyle-{\tau}\widehat{R}_{1s}^{n}\widehat{e}_{s}^{n+\frac{1}{2}}\leq{\tau}^{2}|\widehat{R}_{1s}^{n}|^{2}+\frac{1}{4}|\widehat{e}_{s}^{n+\frac{1}{2}}|^{2},$
$\displaystyle\displaystyle-2\widehat{e}_{s}^{n+\frac{1}{2}}g(u^{n},s^{n})\langle
f(u^{n}),\widehat{e}_{u}^{n+\frac{1}{2}}-e_{u}^{n}\rangle\leq\frac{1}{4}|\widehat{e}_{s}^{n+\frac{1}{2}}|^{2}+C_{5}\|\widehat{e}_{u}^{n+\frac{1}{2}}-e_{u}^{n}\|^{2},$

where $C_{5}>0$ depends on $C_{*}$, $|\Omega|$, $u_{\text{\rm init}}$, and
$\|f\|_{C[-\beta,\beta]}$. By estimating the first term in the similar way as
(3.28), we obtain

$\displaystyle|\widehat{e}_{s}^{n+\frac{1}{2}}|^{2}-|e_{s}^{n}|^{2}+|\widehat{e}_{s}^{n+\frac{1}{2}}-e_{s}^{n}|^{2}$
$\displaystyle\qquad\quad\leq
C_{2}{\tau}(\|e_{u}^{n}\|^{2}+|e_{s}^{n}|^{2}+|\widehat{e}_{s}^{n+\frac{1}{2}}|^{2})+C_{5}\|\widehat{e}_{u}^{n+\frac{1}{2}}-e_{u}^{n}\|^{2}+\frac{1}{2}|\widehat{e}_{s}^{n+\frac{1}{2}}|^{2}+{\tau}^{2}|\widehat{R}_{1s}^{n}|^{2},$

and thus,

$(1-2C_{2}{\tau})|\widehat{e}_{s}^{n+\frac{1}{2}}|^{2}\leq
2|e_{s}^{n}|^{2}+2C_{2}{\tau}(\|e_{u}^{n}\|^{2}+|e_{s}^{n}|^{2})+2C_{5}\|\widehat{e}_{u}^{n+\frac{1}{2}}-e_{u}^{n}\|^{2}+2{\tau}^{2}|\widehat{R}_{1s}^{n}|^{2}.$

When ${\tau}\leq\frac{1}{4C_{2}}$, using (3.41), we can get

$|\widehat{e}_{s}^{n+\frac{1}{2}}|^{2}\leq\|e_{u}^{n}\|^{2}+5|e_{s}^{n}|^{2}+4C_{5}\|\widehat{e}_{u}^{n+\frac{1}{2}}-e_{u}^{n}\|^{2}+4C_{e}^{2}{\tau}^{2}({\tau}+h^{2})^{2}.$
(3.44)

The sum of (3.2.2) multiplied by $C_{5}$ and (3.44) leads to (3.42).

######  Theorem 3.6 (Error estimate of sESAV2)

Given a fixed time $T>0$ and suppose the exact solution $u_{e}$ is smooth
enough on $[0,T]\times\overline{\Omega}$. Assume that
$\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$ and $\|u_{\text{\rm
init}}\|_{\infty}\leq\beta$. If ${\tau}$ and $h$ are small sufficiently and
satisfy (3.36), then we have the error estimate for the sESAV2 scheme (3.32)
as follows:

$\|e_{u}^{n}\|+\|\nabla_{h}e_{u}^{n}\|+|e_{s}^{n}|\leq
C({\tau}^{2}+h^{2}),\qquad 0\leq n\leq\lfloor T/{\tau}\rfloor,$

where the constant $C>0$ depends on $C_{*}$, $|\Omega|$, $T$, $u_{e}$,
$\kappa$, $\varepsilon$, and $\|f\|_{C^{1}[-\beta,\beta]}$ but is independent
of ${\tau}$ and $h$.

###### Proof

The difference between (3.32) and (3.38) leads to

$\displaystyle\delta_{t}e_{u}^{n+1}=\varepsilon^{2}\Delta_{h}e_{u}^{n+\frac{1}{2}}+g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})f(\widehat{u}^{n+\frac{1}{2}})-g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))f(u_{e}(t_{n+\frac{1}{2}}))$
$\displaystyle\qquad\qquad+\kappa\big{(}g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))-g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\big{)}\Big{(}\frac{u_{e}(t_{n+1})+u_{e}(t_{n})}{2}-u_{e}(t_{n+\frac{1}{2}})\Big{)}$
$\displaystyle\qquad\qquad-\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})(e_{u}^{n+\frac{1}{2}}-\widehat{e}_{u}^{n+\frac{1}{2}})-R_{2u}^{n},$
(3.45a)
$\displaystyle\delta_{t}e_{s}^{n+1}=\Big{\langle}g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))f(u_{e}(t_{n+\frac{1}{2}}))-g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})f(\widehat{u}^{n+\frac{1}{2}}),\frac{u_{e}(t_{n+1})-u_{e}(t_{n})}{{\tau}}\Big{\rangle}$
$\displaystyle\qquad\qquad+\kappa\big{(}g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})-g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))\big{)}\Big{\langle}\frac{u_{e}(t_{n+1})+u_{e}(t_{n})}{2}-u_{e}(t_{n+\frac{1}{2}}),\frac{u_{e}(t_{n+1})-u_{e}(t_{n})}{{\tau}}\Big{\rangle}$
$\displaystyle\qquad\qquad-g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
f(\widehat{u}^{n+\frac{1}{2}}),\delta_{t}e_{u}^{n+1}\rangle+\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
u^{n+\frac{1}{2}}-\widehat{u}^{n+\frac{1}{2}},\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\qquad\qquad+\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\Big{\langle}e_{u}^{n+\frac{1}{2}}-\widehat{e}_{u}^{n+\frac{1}{2}},\frac{u_{e}(t_{n+1})-u_{e}(t_{n})}{{\tau}}\Big{\rangle}-R_{2s}^{n}.$
(3.45b)

Taking the discrete inner product of (3.45a) with
$2{\tau}\delta_{t}e_{u}^{n+1}$ and rearranging the term yield

$\displaystyle\varepsilon^{2}\|\nabla_{h}e_{u}^{n+1}\|^{2}-\varepsilon^{2}\|\nabla_{h}e_{u}^{n}\|^{2}+2{\tau}\|\delta_{t}e_{u}^{n+1}\|^{2}$
$\displaystyle\qquad=2{\tau}\langle
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})f(\widehat{u}^{n+\frac{1}{2}})-g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))f(u_{e}(t_{n+\frac{1}{2}})),\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\qquad\quad+2\kappa{\tau}\big{(}g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))-g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\big{)}\Big{\langle}\frac{u_{e}(t_{n+1})+u_{e}(t_{n})}{2}-u_{e}(t_{n+\frac{1}{2}}),\delta_{t}e_{u}^{n+1}\Big{\rangle}$
$\displaystyle\qquad\quad-2\kappa{\tau}g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
e_{u}^{n+\frac{1}{2}}-\widehat{e}_{u}^{n+\frac{1}{2}},\delta_{t}e_{u}^{n+1}\rangle-2{\tau}\langle
R_{2u}^{n},\delta_{t}e_{u}^{n+1}\rangle.$ (3.46)

Since
$g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\geq{\widetilde{G}}_{*}>0$,
we get

$\displaystyle
2\kappa{\tau}g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
e_{u}^{n+\frac{1}{2}}-\widehat{e}_{u}^{n+\frac{1}{2}},\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\qquad=2\kappa{\tau}g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\Big{\langle}\frac{e_{u}^{n+1}-e_{u}^{n}}{2}+e_{u}^{n}-\widehat{e}_{u}^{n+\frac{1}{2}},\delta_{t}e_{u}^{n+1}\Big{\rangle}$
$\displaystyle\qquad=\kappa
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\|e_{u}^{n+1}-e_{u}^{n}\|^{2}+2\kappa{\tau}g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
e_{u}^{n}-\widehat{e}_{u}^{n+\frac{1}{2}},\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\qquad\geq\kappa{\widetilde{G}}_{*}\|e_{u}^{n+1}-e_{u}^{n}\|^{2}+2\kappa{\tau}g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
e_{u}^{n}-\widehat{e}_{u}^{n+\frac{1}{2}},\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\qquad=\kappa{\widetilde{G}}_{*}\|e_{u}^{n+1}\|^{2}-\kappa{\widetilde{G}}_{*}\|e_{u}^{n}\|^{2}-2\kappa{\widetilde{G}}_{*}{\tau}\langle
e_{u}^{n},\delta_{t}e_{u}^{n+1}\rangle+2\kappa{\tau}g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
e_{u}^{n}-\widehat{e}_{u}^{n+\frac{1}{2}},\delta_{t}e_{u}^{n+1}\rangle,$

where we have used (3.20) in the last step. Then, we obtain from (3.46) that

$\displaystyle{\widetilde{G}}_{*}\kappa\|e_{u}^{n+1}\|^{2}-{\widetilde{G}}_{*}\kappa\|e_{u}^{n}\|^{2}+\varepsilon^{2}\|\nabla_{h}e_{u}^{n+1}\|^{2}-\varepsilon^{2}\|\nabla_{h}e_{u}^{n}\|^{2}+2{\tau}\|\delta_{t}e_{u}^{n+1}\|^{2}$
$\displaystyle\qquad=2{\tau}\langle
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})f(\widehat{u}^{n+\frac{1}{2}})-g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))f(u_{e}(t_{n+\frac{1}{2}})),\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\qquad\quad+2\kappa{\tau}\big{(}g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))-g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\big{)}\Big{\langle}\frac{u_{e}(t_{n+1})+u_{e}(t_{n})}{2}-u_{e}(t_{n+\frac{1}{2}}),\delta_{t}e_{u}^{n+1}\Big{\rangle}$
$\displaystyle\qquad\quad+2{\widetilde{G}}_{*}\kappa{\tau}\langle
e_{u}^{n},\delta_{t}e_{u}^{n+1}\rangle+2\kappa{\tau}g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle\widehat{e}_{u}^{n+\frac{1}{2}}-e_{u}^{n},\delta_{t}e_{u}^{n+1}\rangle-2{\tau}\langle
R_{2u}^{n},\delta_{t}e_{u}^{n+1}\rangle.$ (3.47)

For the last three terms in the right-hand side of (3.47), we have
respectively

$\displaystyle
2\kappa{\tau}g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle\widehat{e}_{u}^{n+\frac{1}{2}}-e_{u}^{n},\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\leq
2G^{*}\kappa{\tau}(\|\widehat{e}_{u}^{n+\frac{1}{2}}\|+\|e_{u}^{n}\|)\|\delta_{t}e_{u}^{n+1}\|$
$\displaystyle\leq
6{G^{*}}^{2}\kappa^{2}{\tau}(\|\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2}+\|e_{u}^{n}\|^{2})+\frac{{\tau}}{3}\|\delta_{t}e_{u}^{n+1}\|^{2},$
(3.48) $\displaystyle 2{\widetilde{G}}_{*}\kappa{\tau}\langle
e_{u}^{n},\delta_{t}e_{u}^{n+1}\rangle$ $\displaystyle\leq
3{\widetilde{G}}_{*}^{2}\kappa^{2}{\tau}\|e_{u}^{n}\|^{2}+\frac{{\tau}}{3}\|\delta_{t}e_{u}^{n+1}\|^{2},$
(3.49) $\displaystyle-2{\tau}\langle R_{2u}^{n},\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\leq
3{\tau}\|R_{2u}^{n}\|^{2}+\frac{{\tau}}{3}\|\delta_{t}e_{u}^{n+1}\|^{2}.$
(3.50)

By the energy dissipation and MBP of the sESAV1 substep (3.33), we know that
$\widehat{u}^{n+\frac{1}{2}}$ and $\widehat{s}^{n+\frac{1}{2}}$ are bounded
uniformly. By conducting the similar deductions to the proof of Lemma 4, we
can obtain

$\displaystyle|g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})-g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))|\leq
C_{g}(\|\widehat{e}_{u}^{n+\frac{1}{2}}\|+|\widehat{e}_{s}^{n+\frac{1}{2}}|),$
(3.51a)
$\displaystyle\|g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})f(\widehat{u}^{n+\frac{1}{2}})-g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))f(u_{e}(t_{n+\frac{1}{2}}))\|\leq
C_{g}(\|\widehat{e}_{u}^{n+\frac{1}{2}}\|+|\widehat{e}_{s}^{n+\frac{1}{2}}|),$
(3.51b)

where $C_{g}>0$ is the same constant defined in Lemma 4. Then, the first and
second terms in the right-hand side of (3.47) can be bounded respectively as

$\displaystyle 2{\tau}\langle
g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})f(\widehat{u}^{n+\frac{1}{2}})-g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))f(u_{e}(t_{n+\frac{1}{2}})),\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\qquad\leq
2C_{g}{\tau}(\|\widehat{e}_{u}^{n+\frac{1}{2}}\|+|\widehat{e}_{s}^{n+\frac{1}{2}}|)\|\delta_{t}e_{u}^{n+1}\|$
$\displaystyle\qquad\leq
6C_{g}^{2}{\tau}(\|\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2}+|\widehat{e}_{s}^{n+\frac{1}{2}}|^{2})+\frac{{\tau}}{3}\|\delta_{t}e_{u}^{n+1}\|^{2},$
(3.52)

and

$\displaystyle
2\kappa{\tau}\big{(}g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))-g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\big{)}\Big{\langle}\frac{u_{e}(t_{n+1})+u_{e}(t_{n})}{2}-u_{e}(t_{n+\frac{1}{2}}),\delta_{t}e_{u}^{n+1}\Big{\rangle}$
$\displaystyle\qquad\leq
C_{g}\kappa{\tau}(\|\widehat{e}_{u}^{n+\frac{1}{2}}\|+|\widehat{e}_{s}^{n+\frac{1}{2}}|)(\|u_{e}(t_{n+1})\|+\|u_{e}(t_{n})\|+2\|u_{e}(t_{n+\frac{1}{2}})\|)\|\delta_{t}e_{u}^{n+1}\|$
$\displaystyle\qquad\leq
C_{1}\kappa^{2}{\tau}(\|\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2}+|\widehat{e}_{s}^{n+\frac{1}{2}}|^{2})+\frac{{\tau}}{3}\|\delta_{t}e_{u}^{n+1}\|^{2},$
(3.53)

where $C_{1}>0$ has the same dependence as the constant $C_{1}$ used in (3.23)
but may have a different value. Substituting (3.48)–(3.53) into (3.47) leads
to

$\displaystyle{\widetilde{G}}_{*}\kappa\|e_{u}^{n+1}\|^{2}-{\widetilde{G}}_{*}\kappa\|e_{u}^{n}\|^{2}+\varepsilon^{2}\|\nabla_{h}e_{u}^{n+1}\|^{2}-\varepsilon^{2}\|\nabla_{h}e_{u}^{n}\|^{2}+\frac{{\tau}}{3}\|\delta_{t}e_{u}^{n+1}\|^{2}$
$\displaystyle\qquad\leq(6{G^{*}}^{2}\kappa^{2}+6C_{g}^{2}+C_{1}\kappa^{2}){\tau}\|\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2}+(6C_{g}^{2}+C_{1}\kappa^{2}){\tau}|\widehat{e}_{s}^{n+\frac{1}{2}}|^{2}$
$\displaystyle\qquad\quad+(3{\widetilde{G}}_{*}^{2}+6{G^{*}}^{2})\kappa^{2}{\tau}\|e_{u}^{n}\|^{2}+3{\tau}\|R_{2u}^{n}\|^{2}.$
(3.54)

Multiplying (3.45b) by $2{\tau}e_{s}^{n+1}$ yields

$\displaystyle|e_{s}^{n+1}|^{2}-|e_{s}^{n}|^{2}+|e_{s}^{n+1}-e_{s}^{n}|^{2}$
$\displaystyle\quad=2e_{s}^{n+1}\langle
g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))f(u_{e}(t_{n+\frac{1}{2}}))-g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})f(\widehat{u}^{n+\frac{1}{2}}),u_{e}(t_{n+1})-u_{e}(t_{n})\rangle$
$\displaystyle\quad\quad+\kappa
e_{s}^{n+1}\big{(}g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})-g(u_{e}(t_{n+\frac{1}{2}}),s_{e}(t_{n+\frac{1}{2}}))\big{)}\langle
u_{e}(t_{n+1})+u_{e}(t_{n})-2u_{e}(t_{n+\frac{1}{2}}),u_{e}(t_{n+1})-u_{e}(t_{n})\rangle$
$\displaystyle\quad\quad-2{\tau}e_{s}^{n+1}g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
f(\widehat{u}^{n+\frac{1}{2}}),\delta_{t}e_{u}^{n+1}\rangle+2\kappa{\tau}e_{s}^{n+1}g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
u^{n+\frac{1}{2}}-\widehat{u}^{n+\frac{1}{2}},\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\quad\quad+2\kappa
e_{s}^{n+1}g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
e_{u}^{n+\frac{1}{2}}-\widehat{e}_{u}^{n+\frac{1}{2}},u_{e}(t_{n+1})-u_{e}(t_{n})\rangle-2{\tau}R_{2s}^{n}e_{s}^{n+1}.$
(3.55)

The last term in the right-hand side of (3.55) can be estimated by

$-2{\tau}R_{2s}^{n}e_{s}^{n+1}\leq{\tau}|e_{s}^{n+1}|^{2}+{\tau}|R_{2s}^{n}|^{2}.$
(3.56)

By the boundedness of $u^{n+\frac{1}{2}}$, $\widehat{u}^{n+\frac{1}{2}}$, and
$\widehat{s}^{n+\frac{1}{2}}$, the sum of the third and fourth terms in the
right-hand side of (3.55) can be estimated similarly to (3.1.2) as follows:

$\displaystyle-2{\tau}e_{s}^{n+1}g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
f(\widehat{u}^{n+\frac{1}{2}}),\delta_{t}e_{u}^{n+1}\rangle+2\kappa{\tau}e_{s}^{n+1}g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})\langle
u^{n+\frac{1}{2}}-\widehat{u}^{n+\frac{1}{2}},\delta_{t}e_{u}^{n+1}\rangle$
$\displaystyle\qquad\leq
2G^{*}\kappa{\tau}(\|u^{n+\frac{1}{2}}\|+\|\widehat{u}^{n+\frac{1}{2}}\|)|e_{s}^{n+1}|\|\delta_{t}e_{u}^{n+1}\|+2G^{*}{\tau}\|f(\widehat{u}^{n+\frac{1}{2}})\||e_{s}^{n+1}|\|\delta_{t}e_{u}^{n+1}\|$
$\displaystyle\qquad\leq
C_{6}{\tau}|e_{s}^{n+1}|^{2}+\frac{{\tau}}{3}\|\delta_{t}e_{u}^{n+1}\|^{2},$
(3.57)

where $C_{6}>0$ depends on $C_{*}$, $|\Omega|$, $u_{\text{\rm init}}$,
$\kappa$, and $\|f\|_{C[-\beta,\beta]}$. Then, using the facts that
$\|u_{e}(t_{n+1})-u_{e}(t_{n})\|\leq C{\tau}$,
$\|u_{e}(t_{n+1})+u_{e}(t_{n})-2u_{e}(t_{n+\frac{1}{2}})\|\leq C\tau^{2}$
(where $C>0$ is a constant due to smoothness of $u_{e}$), and the inequalities
(3.51), in the similar spirit of deriving (3.28), the sum of the first, second
and fifth terms in the right-hand side of (3.55) can be bounded above by

${\tau}\big{(}\|\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2}+|\widehat{e}_{s}^{n+\frac{1}{2}}|^{2}+\|e_{u}^{n}\|^{2}+\|e_{u}^{n+1}\|^{2}+|e_{s}^{n+1}|^{2}\big{)}$
(3.58)

multiplied with a positive constant depending on $C_{*}$, $|\Omega|$, $u_{e}$,
$\kappa$, and $\|f\|_{C^{1}[-\beta,\beta]}$. Combining (3.55) with
(3.56)–(3.58), we obtain

$\displaystyle|e_{s}^{n+1}|^{2}-|e_{s}^{n}|^{2}$ $\displaystyle\leq
C_{7}{\tau}\big{(}\|\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2}+|\widehat{e}_{s}^{n+\frac{1}{2}}|^{2}+\|e_{u}^{n}\|^{2}+\|e_{u}^{n+1}\|^{2}+|e_{s}^{n+1}|^{2}\big{)}$
$\displaystyle\quad+\frac{{\tau}}{3}\|\delta_{t}e_{u}^{n+1}\|^{2}+{\tau}|R_{2s}^{n}|^{2}$
(3.59)

with $C_{7}$ depending on $C_{*}$, $|\Omega|$, $u_{e}$, $\kappa$, and
$\|f\|_{C^{1}[-\beta,\beta]}$.

Adding (3.54) and (3.59), we obtain

$\displaystyle{\widetilde{G}}_{*}\kappa(\|e_{u}^{n+1}\|^{2}-\|e_{u}^{n}\|^{2})+\varepsilon^{2}(\|\nabla_{h}e_{u}^{n+1}\|^{2}-\|\nabla_{h}e_{u}^{n}\|^{2})+(|e_{s}^{n+1}|^{2}-|e_{s}^{n}|^{2})$
$\displaystyle\qquad\leq
C_{8}{\tau}\big{(}\|\widehat{e}_{u}^{n+\frac{1}{2}}\|^{2}+|\widehat{e}_{s}^{n+\frac{1}{2}}|^{2}+\|e_{u}^{n}\|^{2}+\|e_{u}^{n+1}\|^{2}+|e_{s}^{n+1}|^{2}\big{)}+3{\tau}\|R_{2u}^{n}\|^{2}+{\tau}|R_{2s}^{n}|^{2},$
(3.60)

where $C_{8}>0$ depends on $C_{*}$, $|\Omega|$, $u_{e}$, $\kappa$, and
$\|f\|_{C^{1}[-\beta,\beta]}$. Substituting (3.42) into (3.60) and using the
estimate (3.39), when ${\tau}\leq 1$, we have

$\displaystyle{\widetilde{G}}_{*}\kappa(\|e_{u}^{n+1}\|^{2}-\|e_{u}^{n}\|^{2})+\varepsilon^{2}(\|\nabla_{h}e_{u}^{n+1}\|^{2}-\|\nabla_{h}e_{u}^{n}\|^{2})+(|e_{s}^{n+1}|^{2}-|e_{s}^{n}|^{2})$
$\displaystyle\qquad\leq
C_{8}(\widehat{C}+1){\tau}(\|e_{u}^{n}\|^{2}+\|e_{u}^{n+1}\|^{2}+|e_{s}^{n+1}|^{2})+(C_{8}\widehat{C}+4)C_{e}^{2}{\tau}({\tau}^{2}+h^{2})^{2}.$

When ${\tau}$ is small sufficiently, similar to the last paragraph in the
proof of Theorem 3.3, applying the discrete Gronwall’s inequality yields

${\widetilde{G}}_{*}\kappa\|e_{u}^{n}\|^{2}+\varepsilon^{2}\|\nabla_{h}e_{u}^{n}\|^{2}+|e_{s}^{n}|^{2}\leq
C({\tau}^{2}+h^{2})^{2},$

which completes the proof.

##  4 Numerical experiments

This section is devoted to numerical tests and comparisons between the
proposed sESAV schemes and existing SAV schemes listed in Sections 2.2 and
2.3. We consider the Allen–Cahn equation (1.1) in two-dimensional spatial
domain $\Omega=(0,1)\times(0,1)$ equipped with periodic boundary conditions,
so that the schemes can be solved efficiently by the fast Fourier transform.
We take two types of commonly-used nonlinear functions $f(u)$. One is given by

$f(u)=-F^{\prime}(u)=u-u^{3}$ (4.1)

with $F$ being the double-well potential

$F(u)=\frac{1}{4}(u^{2}-1)^{2}.$

In this case, one has $\beta=1$ and $\|f^{\prime}\|_{C[-1,1]}=2$. The constant
$C_{0}$ in (2.4) is $C_{0}=\frac{1}{4}(\kappa^{2}+2\kappa)$. The other one is
determined by the Flory–Huggins potential

$F(u)=\frac{\theta}{2}[(1+u)\ln(1+u)+(1-u)\ln(1-u)]-\frac{\theta_{c}}{2}u^{2},$

which gives

$f(u)=-F^{\prime}(u)=\frac{\theta}{2}\ln\frac{1-u}{1+u}+\theta_{c}u,$ (4.2)

where $\theta_{c}>\theta>0$. In the following experiments, we set $\theta=0.8$
and $\theta_{c}=1.6$, then the positive root of $f(\rho)=0$ gives us
$\beta\approx 0.9575$, and $\|f^{\prime}\|_{C[-\beta,\beta]}\approx 8.02$. The
constant $C_{0}$ in (2.4) is then determined by
$C_{0}=-F(\alpha)+\frac{\kappa}{2}\alpha^{2}$, where $\alpha>0$ solves
$f(\alpha)+\kappa\alpha=0$.

###  4.1 Convergence in time

We first verify the convergence order in time for the proposed sESAV schemes.
Let us set $\varepsilon=0.01$ in (1.1) and take a smooth initial value

$u_{\text{\rm init}}(x,y)=0.1\sin(2\pi x)\sin(2\pi y).$

The temporal convergence tests are conducted by fixing the spatial mesh size
$h=1/512$. As requested by the stabilizing condition
$\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$, we set $\kappa=2$ for the
double-well potential case (i.e., $f(u)$ given by (4.1)) and $\kappa=8.02$ for
the Flory–Huggins potential case (i.e., $f(u)$ given by (4.2)). We compute the
numerical solutions at $t=2$ using the sESAV1 and sESAV2 schemes with various
time step sizes ${\tau}=2^{-k}$, $k=4,5,\dots,12$. To compute the numerical
errors, we treat the sESAV2 solution obtained by ${\tau}=0.1\times 2^{-12}$ as
the benchmark solution. Figure 1 shows the relation between the $L^{2}$-norm
error and the time step size, where the left picture corresponds to the
double-well potential case and the right one for the Flory–Huggins potential
case. The first-order temporal accuracy for sESAV1 and the second-order for
sESAV2 are observed for both cases as expected.

Figure 1: The $L^{2}$-norm errors vs. the time step size produced by the
proposed sESAV1 and sESAV2 schemes for the double-well potential case (4.1)
(left) and the Flory–Huggins potential case (4.2) (right).

###  4.2 Comparisons with existing SAV schemes

In the following numerical experiments, we compare the proposed sESAV schemes
with classic SAV and ESAV schemes by focusing on the MBP and energy
dissipation law. While various modified energies are introduced as
approximations of the original energy in discrete settings in order to
facilitate the proof of energy dissipation law, the original one possesses the
most accurate physical meaning for the model problem. Therefore, we are
concerned about the behavior of the original (discrete) energy $E_{h}(u)$
defined in (2.2) for reflecting the phase transition process. The dynamic
process considered usually needs a long-time evolution to reach the steady
state; here we conduct simulations in a short time interval for the comparison
among these schemes.

Let us still consider the problem (1.1) with $\varepsilon=0.01$. We adopt the
uniform spatial mesh with $h=1/512$ and give the initial value by random
numbers between $-0.8$ and $0.8$ on each mesh point. We then set the time step
size ${\tau}=0.01$, and compute the numerical solutions by using the sESAV
schemes, the classic SAV schemes (SAV1 and SAV2), and the ESAV schemes (ESAV1
and ESAV2). Note that we set $\delta=C_{0}+0.01$ for the classic SAV schemes
(2.6) and (2.7). For all comparison experiments, we will consider two settings
for the stabilizing parameter: $\kappa=\|f^{\prime}\|_{C[-\beta,\beta]}$ and
$\kappa=\frac{1}{2}\|f^{\prime}\|_{C[-\beta,\beta]}$, where the former one
satisfies the requirement for the MBP preservation for the sESAV schemes and
the latter one was adopted in ShXuYa19  for the classic SAV schemes. In
addition, we take the numerical results obtained by the IFRK4 scheme
JuLiQiYa21  with the small time step size $10^{-4}$ as the benchmark solution.

First, we test the double-well potential case (4.1), and correspondingly, set
the stabilizing parameter $\kappa=1$ and $\kappa=2$ respectively to carry out
the experiments. Figure 2 shows the evolutions of the supremum norms and the
energies of simulated solutions computed by the sESAV1, SAV1, and ESAV1
schemes. For either $\kappa=1$ or $\kappa=2$, the sESAV1 scheme preserves the
MBP, while the supremum norms of the SAV1 and ESAV1 solutions obviously evolve
beyond $1$, which means that the MBP is violated. The energy dissipation are
observed for these three schemes, where the sESAV1 scheme provides the most
accurate result. In addition, the larger $\kappa$ leads to larger errors in
the results, especially for the ESAV1 scheme. Figure 3 plots corresponding
results computed by the second-order schemes. Again, only the sESAV2 scheme
preserves the MBP and the energy dissipation perfectly. The SAV2 and ESAV2
solutions evolve beyond $1$ but closer to $1$ than their first-order results
due to the higher-order temporal accuracy.

Figure 2: Evolutions of the supremum norms and the energies of simulated
solutions computed by the sESAV1, SAV1, and ESAV1 schemes with ${\tau}=0.01$
and $\kappa=1$ (top row) or $\kappa=2$ (bottom row) for the double-well
potential case.

Figure 3: Evolutions of the supremum norms and the energies of simulated
solutions computed by the sESAV2, SAV2, and ESAV2 schemes with ${\tau}=0.01$
and $\kappa=1$ (top row) or $\kappa=2$ (bottom row) for the double-well
potential case.

Next, we test the Flory–Huggins potential case (4.2) and correspondingly set
$\kappa=4.01$ and $\kappa=8.02$ respectively. Figures 4 and 5 present the
evolutions of the supremum norms and the energies of simulated solutions
obtained by the first- and second-order schemes, respectively. Similar to the
double-well potential case, only the sESAV schemes preserve the MBP and the
energy dissipation law as expected. The SAV1, ESAV1, and ESAV2 schemes, having
the supremum norms beyond the theoretical bound $0.9575$, lead to inaccurate
dynamic processes. Especially, the ESAV1 solution with $\kappa=8.02$ evolves
beyond $1$, which yields complex numbers due to the existence of the
logarithmic term and gives the completely wrong dynamics. For the SAV2
solutions, the dynamic processes look moderately correct according to the
energy evolutions. Moreover, it is interesting that the supremum norm goes
larger than the desired bound for $\kappa=8.02$ while it does not exceed for
$\kappa=4.01$, but both results are still a bit away from the expected value
$0.9575$.

Figure 4: Evolutions of the supremum norms and the energies of simulated
solutions computed by the sESAV1, SAV1, and ESAV1 schemes with ${\tau}=0.01$
and $\kappa=4.01$ (top row) or $\kappa=8.02$ (bottom row) for the
Flory–Huggins potential case.

Figure 5: Evolutions of the supremum norms and the energies of simulated
solutions computed by the sESAV2, SAV2, and ESAV2 schemes with ${\tau}=0.01$
and $\kappa=4.01$ (top row) or $\kappa=8.02$ (bottom row) for the
Flory–Huggins potential case.

###  4.3 Long-time coarsening dynamics simulations

Now we study the coarsening dynamics driven by the Allen–Cahn equation (1.1)
with $\varepsilon=0.01$. The spatial mesh size is $h=1/512$ and the initial
state is given by random numbers between $-0.8$ and $0.8$. We adopt the sESAV2
scheme with ${\tau}=0.01$ to simulate the long-time coarsening process. By the
comparisons shown above, we know that ${\tau}=0.01$ is sufficient to provide
accurate numerical results. The steady state of the coarsening dynamics is a
constant state $u\equiv\beta$ or $u\equiv-\beta$. When the absolute difference
between the energies at the two consecutive moments is smaller than the
tolerance value $10^{-8}$, we regard the dynamics as reaching its steady
state.

For the double-well potential case $f(u)$, we set $\kappa=2$ and the phase
structures captured at some moments are presented in Figure 6, and the
constant steady state $u\equiv-1$ is reached at around $t=604$. The left
picture given in Figure 7 implies the preservation of the MBP during the whole
phase transition process. The energy evolution is plotted in the right graph
of Figure 7, which states the energy dissipation of the process. For the
Flory–Huggins potential case, we set $\kappa=8.02$ and the simulated results
are shown in Figures 8 and 9. We observe that the steady state is reached at
around $t=602$ and the whole process of phase separation is similar to that of
the double-well potential case. Those results are almost identical to those
produced using the IFRK4 scheme in JuLiQiYa21 .

Figure 6: Simulated phase structures at $t=4$, $6$, $10$, $30$, $100$, and
$300$, respectively (left to right and top to bottom) by the sESAV2 scheme
with ${\tau}=0.01$ and $\kappa=2$ for the coarsening dynamics of the double-
well potential case.

Figure 7: Evolutions of the supremum norm (left) and the energy (right) for
the coarsening dynamics of the double-well potential case.

Figure 8: Simulated phase structures at $t=4$, $6$, $10$, $30$, $100$, and
$300$, respectively (left to right and top to bottom) by the sESAV2 scheme
with ${\tau}=0.01$ and $\kappa=8.02$ for the coarsening dynamics of the
Flory–Huggins potential case.

Figure 9: Evolutions of the supremum norm (left) and the energy (right) for
the coarsening dynamics of the Flory–Huggins potential case.

##  5 Conclusion

In this paper, we study MBP-preserving and energy dissipative schemes for the
Allen–Cahn type equations by combining the ESAV approach with stabilizing
technique. We present first- and second-order sESAV schemes and prove their
MBP preservation, energy dissipation, and error estimates. The main results
and observations include two aspects. First, we choose the ESAV approach
rather than the classic SAV approach, since the coefficient ($g(u^{n},s^{n})$
or $g(\widehat{u}^{n+\frac{1}{2}},\widehat{s}^{n+\frac{1}{2}})$) of the
nonlinear term is positive automatically in the former one while the sign of
the corresponding coefficient is uncertain for the later one. Second, to
guarantee the MBP-preserving property, we add the stabilization term as an
extra artificial term, that is, add and subtract a linear term in the scheme
instead of a quadratic term in the energy functional; they are equivalent
mutually for the classic stabilization or convex splitting method, but not for
the SAV approach. Moreover, we find that the MBP preservation and the energy
dissipation of the sESAV schemes can be established in parallel and
independently, unlike the purely stabilized semi-implicit scheme discussed in
TaYa16  where the MBP is needed first to bound the nonlinear term in the proof
of the stability with respect to the original energy. Since the schemes we
studied are all one-step methods, adaptive time-stepping strategies (such as
QiaoZhTa11 ) can be inherently adopted to accelerate the computation.

Some generalizations can be carried out by replacing the Laplace operator in
(1.1) by some analogues, for instance, the nonlocal diffusion DuGuLeZh12  and
the fractional Laplace operators SmakoKiMa93  which also satisfy the semigroup
property with their discretizations satisfying the analogues of Lemma 1.
Furthermore, the proposed stabilizing approaches in this paper can be
naturally extended to many other type of gradient flow problems, which can be
handled by the existing SAV schemes. For example, the fourth-order
Cahn–Hilliard equation is the $H^{-1}$ gradient flow of the energy functional
(1.2), and satisfies the same energy dissipation law as the Allen–Cahn
equation. The MBP is not valid anymore, but the solution is still $L^{\infty}$
stable. In the similar spirit of this paper, it is interesting to develop the
sESAV schemes for the Cahn–Hilliard equation, and the discrete $L^{\infty}$
stability of the sESAV solution can be established by combining the high-order
consistency analysis and stability estimate, as done in GuWaWi14 ; LiQiWa21 ,
which will also be one of our future works.

## References

  * (1) S. M. Allen and J. W. Cahn, A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta Metall., 27 (1979), 1085–1095.
  * (2) G. Akrivis, B. Y. Li, and D. F. Li, Energy-decaying extrapolated RK-SAV methods for the Allen–Cahn and Cahn–Hilliard equations, SIAM J. Sci. Comput., 41 (2019), A3703–A3727.
  * (3) P. W. Bates, On some nonlocal evolution equations arising in materials science, Fields Inst. Commun., 48 (2006), 13–52.
  * (4) C. J. Chen and X. F. Yang, Fast, provably unconditionally energy stable, and second-order accurate algorithms for the anisotropic Cahn–Hilliard model, Comput. Methods Appl. Mech. Engrg., 351 (2019), 35–59.
  * (5) Q. Cheng, C. Liu, and J. Shen, A new Lagrange multiplier approach for gradient flows, Comput. Methods Appl. Mech. Engrg., 367 (2020), 113070.
  * (6) Q. Cheng, C. Liu, and J. Shen, Generalized SAV approaches for gradient systems, J. Comput. Appl. Math., 394 (2021), 113532.
  * (7) Q. Du, M. Gunzburger, R. B. Lehoucq, and K. Zhou, Analysis and approximation of nonlocal diffusion problems with volume constraints, SIAM Rev., 54 (2012), 667–696.
  * (8) Q. Du, L. Ju, X. Li, and Z. H. Qiao, Maximum principle preserving exponential time differencing schemes for the nonlocal Allen–Cahn equation, SIAM J. Numer. Anal., 57 (2019), 875–898.
  * (9) Q. Du, L. Ju, X. Li, and Z. H. Qiao, Maximum bound principles for a class of semilinear parabolic equations and exponential time-differencing schemes, SIAM Rev., 63 (2021), 317–359.
  * (10) Q. Du, J. Yang, and Z. Zhou, Time-fractional Allen–Cahn equations: analysis and numerical methods, J. Sci. Comput., 85 (2020), 42.
  * (11) X. L. Feng, T. Tang, and J. Yang, Stabilized Crank–Nicolson/Adams–Bashforth schemes for phase field models, East Asian J. Appl. Math., 3 (2013), 59–80.
  * (12) Z. Guan, C. Wang, and S. M. Wise, A convergent convex splitting scheme for the periodic nonlocal Cahn–Hilliard equation, Numer. Math., 128 (2014), 377–406.
  * (13) C. F. Gui and M. F. Zhao, Traveling wave solutions of Allen–Cahn equation with a fractional Laplacian, Ann. Inst. H. Poincaré-An., 32 (2015), 785–812.
  * (14) D. M. Hou, M. Azaiez, and C. J. Xu, A variant of scalar auxiliary variable approaches for gradient flows, J. Comput. Phys., 395 (2019), 307–332.
  * (15) T. L. Hou and H. T. Leng, Numerical analysis of a stabilized Crank–Nicolson/Adams–Bashforth finite difference scheme for Allen–Cahn equations, Appl. Math. Lett., 102 (2020), 106150.
  * (16) T. L. Hou, T. Tang, and J. Yang, Numerical analysis of fully discretized Crank–Nicolson scheme for fractional-in-space Allen–Cahn equations, J. Sci. Comput., 72 (2017), 1214–1231.
  * (17) F. K. Huang, J. Shen, and Z. G. Yang, A highly efficient and accurate new scalar auxiliary variable approach for gradient flows, SIAM J. Sci. Comput., 42 (2020), A2514–A2536.
  * (18) L. Ju, X. Li, Z. H. Qiao, and J. Yang, Maximum bound principle preserving integrating factor Runge–Kutta methods for semilinear parabolic equations, J. Comput. Phys., 439 (2021), 110405.
  * (19) L. Ju, X. Li, Z. H. Qiao, and H. Zhang, Energy stability and error estimates of exponential time differencing schemes for the epitaxial growth model without slope selection, Math. Comp., 87 (2018), 1859–1885.
  * (20) L. Ju, J. Zhang, and Q. Du, Fast and accurate algorithms for simulating coarsening dynamics of Cahn–Hilliard equations, Comput. Mater. Sci., 108 (2015), 272–282.
  * (21) J. W. Li, L. Ju, Y. Y. Cai, and X. L. Feng, Unconditionally maximum bound principle preserving linear schemes for the conservative Allen–Cahn equation with nonlocal constraint, J. Sci. Comput., 87 (2021), 98.
  * (22) J. W. Li, X. Li, L. Ju, and X. L. Feng, Stabilized integrating factor Runge–Kutta method and unconditional preservation of maximum bound principle, SIAM J. Sci. Comput., 43 (2021), A1780–A1802.
  * (23) X. Li, Z. H. Qiao, and C. Wang, Convergence analysis for a stabilized linear semi-implicit numerical scheme for the nonlocal Cahn–Hilliard equation, Math. Comp., 90 (2021), 171–188.
  * (24) H. L. Liao, T. Tang, and T. Zhou, On energy stable, maximum-principle preserving, second-order BDF scheme with variable steps for the Allen–Cahn equation, SIAM J. Numer. Anal., 58 (2020), 2294–2314.
  * (25) Z. G. Liu and X. L. Li, The exponential scalar auxiliary variable (E-SAV) approach for phase field models and its explicit computing, SIAM J. Sci. Comput., 42 (2020), B630–B655.
  * (26) Z. H. Qiao, Z. R. Zhang, and T. Tang, An adaptive time-stepping strategy for the molecular beam epitaxy models, SIAM J. Sci. Comput., 33 (2011), 1395–1414.
  * (27) S. G. Samko, A. A. Kilbas, and O. I. Marichev, Fractional Integrals and Derivatives, Gordon and Breach, Yverdon, 1993.
  * (28) J. Shen, T. Tang, and J. Yang, On the maximum principle preserving schemes for the generalized Allen–Cahn equation, Commun. Math. Sci., 14 (2016), 1517–1534.
  * (29) J. Shen, C. Wang, X. M. Wang, and S. M. Wise, Second-order convex splitting schemes for gradient flows with Ehrlich–Schwoebel type energy: application to thin film epitaxy, SIAM J. Numer. Anal., 50 (2012), 105–125.
  * (30) J. Shen and J. Xu, Convergence and error analysis for the scalar auxiliary variable (SAV) schemes to gradient flows, SIAM J. Numer. Anal., 56 (2018), 2895–2912.
  * (31) J. Shen, J. Xu, and J. Yang, The scalar auxiliary variable (SAV) approach for gradient flows, J. Comput. Phys., 353 (2018), 407–416.
  * (32) J. Shen, J. Xu, and J. Yang, A new class of efficient and robust energy stable schemes for gradient flows, SIAM Rev., 61 (2019), 474–506.
  * (33) J. Shen and X. F. Yang, Numerical approximations of Allen–Cahn and Cahn–Hilliard equations, Discrete Contin. Dyn. Syst., 28 (2010), 1669–1691.
  * (34) T. Tang and J. Yang, Implicit-explicit scheme for the Allen–Cahn equation preserves the maximum principle, J. Comput. Math., 34 (2016), 471–481.
  * (35) S. M. Wise, C. Wang, and J. S. Lowengrub, An energy stable and convergent finite difference scheme for the phase field crystal equation, SIAM J. Numer. Anal., 47 (2009), 2269–2288.
  * (36) X. F. Xiao, X. L. Feng, and J. Y. Yuan, The stabilized semi-implicit finite element method for the surface Allen–Cahn equation, Discrete Contin. Dyn. Syst. Ser. B, 22 (2017), 2857–2877.
  * (37) C. J. Xu and T. Tang, Stability analysis of large time-stepping methods for epitaxial growth models, SIAM J. Numer. Anal., 44 (2006), 1759–1779.
  * (38) Z. Xu, X. F. Yang, H. Zhang, and Z. Q. Xie, Efficient and linear schemes for anisotropic Cahn–Hilliard model using the stabilized-invariant energy quadratization (S-IEQ) approach, Comput. Phys. Commun., 238 (2019), 36–49.
  * (39) X. F. Yang, Linear, first and second-order, unconditionally energy stable numerical schemes for the phase field model of homopolymer blends, J. Comput. Phys., 327 (2016), 294–316.
  * (40) X. F. Yang and G. D. Zhang, Convergence analysis for the invariant energy quadratization (IEQ) schemes for solving the Cahn–Hilliard and Allen–Cahn equations with general nonlinear potential, J. Sci. Comput., 82 (2020), 55.