# Global-in-time energy stability: a powerful analysis tool for the gradient
flow problem without maximum principle or Lipschitz assumption

Jingwei Sun  Department of Mathematics, National University of Defense
Technology, Changsha 410073, Hunan Province, China.  Haifeng Wang  Department
of Mathematics, National University of Defense Technology, Changsha 410073,
Hunan Province, China.  Hong Zhang Corresponding author.
E-mail address<EMAIL_ADDRESS>Department of Mathematics, National
University of Defense Technology, Changsha 410073, Hunan Province, China.  Xu
Qian  College of Mathematics Science, University of Wuhan, Wuhan 430072, Hubei
Province, China.  Songhe Song  Department of Mathematics, National University
of Defense Technology, Changsha 410073, Hunan Province, China.

###### Abstract

Before proving (unconditional) energy stability for gradient flows, most
existing studies either require a strong Lipschitz condition regarding the
non-linearity or certain $L^{\infty}$ bounds on the numerical solutions (the
maximum principle). However, proving energy stability without such premises is
a very challenging task. In this paper, we aim to develop a novel analytical
tool, namely global-in-time energy stability, to demonstrate energy
dissipation without assuming any strong Lipschitz condition or $L^{\infty}$
boundedness. The fourth-order-in-space Swift-Hohenberg equation is used to
elucidate the theoretical results in detail. We also propose a temporal
second-order accurate scheme for efficiently solving such a strongly stiff
equation. Furthermore, we present the corresponding optimal $L^{2}$ error
estimate and provide several numerical simulations to demonstrate the
dynamics.

Mathematics Subject Classification: Primary 35K30, 35K55, 65L06, 65T40
Key words: Global-in-time energy stability, gradient flow, Lispchitz
assumption, maximum principle, original energy, exponential Runge–Kutta,
Swift–Hohenberg equation.

##  1 Introduction

Many physical problems can be modeled by PDEs that take the form of gradient
flows, which are often derived from the second law of thermodynamics. It is
well-known that a gradient flow is determined by not only the driving free
energy but also the dissipation mechanism. Given a free energy functional
$E(u)$, denote its variational derivative as $\mu=\delta E/\delta u$. The
general form of the gradient flow can be written as

$\frac{\partial u}{\partial t}=\mathcal{G}\mu,$ (1.1)

equipped with suitable boundary conditions. In the above, a non-positive
symmetric operator $\mathcal{G}$ gives the dissipation mechanism, thus the
free energy is non-increasing:

$\frac{\mathrm{d}E}{\mathrm{d}t}=\langle\mu,\mathcal{G}\mu\rangle\leq 0,$
(1.2)

where
$\langle\phi,\psi\rangle=\int_{\Omega}\phi\cdot\psi\mathrm{d}\mathbf{x}$.
Familiar dissipative operators $\mathcal{G}$ include but not limited to
$-I,\leavevmode\nobreak\ \Delta,\leavevmode\nobreak\
-(-\Delta)^{\alpha}\leavevmode\nobreak\ (0<\alpha<1)$, corresponding to
$L^{2},\leavevmode\nobreak\ H^{-1}$ and nonlocal $H^{-\alpha}$ gradient flows,
respectively. In general, the free energy functional contains linear and
nonlinear terms, which we write explicitly as

$E(u)=\frac{1}{2}\langle u,\mathcal{L}u\rangle+\langle F(u),1\rangle,$ (1.3)

where $\mathcal{L}$ is a symmetric non-negative linear operator, and $F(u)$
represents the nonlinear potential. In particular, if the energy functional
associated with (1.1) is

$E_{\varepsilon}(u)=\int_{\Omega}\left(\frac{\varepsilon^{2}}{2}|\nabla
u|^{2}+F(u)\right)\,\mathrm{d}\mathbf{x},$ (1.4)

where $\varepsilon$ is a certain positive parameter, the well-known Allen–Cahn
(AC, $L^{2}$ gradient flow) and Cahn–Hilliard (CH, $H^{-1}$ gradient flow)
equations could be derived with different nonlinear potentials. Moreover, it
follows from (1.2) that

$E_{\varepsilon}(u(t))\leq E_{\varepsilon}(u(s)),\quad\forall t\leq s,$

which gives a prior control of $H^{1}$-norm of the solution.

The energy dissipation property (1.2) is always viewed as a key criterion for
designing efficient and long-time stable numerical schemes, such as [2, 4, 18,
15, 19, 13]. Among them, Shen and Yang [28] proved energy stability for the AC
and CH equations by a nonlinear term truncation technique. More precisely, it
is assumed that

$\max\limits_{u\in\mathbb{R}}|\tilde{f}{{}^{\prime}}(u)|\leq\beta,\leavevmode\nobreak\
\tilde{f}{{}^{\prime}}(u)\leavevmode\nobreak\ \text{is a suitable modification
of the original function derivative}\leavevmode\nobreak\ f{{}^{\prime}}(u),$

which is what we referred to as the Lipschitz assumption on the non-linearity
in the abstract. Also, see the same strategy for example [25, 26, 27, 29] and
the references therein. As Li and Qiao [34] mentioned, the main drawback of
the aforementioned numerical developments is that to obtain energy stability,
one either makes a Lipschitz assumption on the non-linearity such as [28], or
one assumes some additional $L^{\infty}$ bounds on the numerical solution,
which is automatically satisfied by AC type equations with the standard
maximum principle. It is very desirable to remove these technical obstacles
and establish a more reasonable energy stability theory. Therefore, the
purpose of this work is to prove energy stability on applying an efficient
scheme to the gradient flow without maximum principle or any Lipschitz
assumption.

In fact, very few work has been devoted to such analysis. This is partly due
to some technical difficulties related to the $H^{m}$-norm $(m\geq 1)$ prior
control and $L^{\infty}$ bounds of the solutions, making it difficult to
handle nonlinear terms optimally. Notably, by utilizing a log-type
interpolation inequality, Li et al. [34, 35] provided rigorous analyses for
the CH equation several years ago, employing first-order stabilized semi-
implicit (SSI1) and second-order backward differential formula (BDF2) time-
stepping schemes. In this paper, however, we shall adopt a different approach
by considering another gradient flow, the Swift–Hohenberg (SH, $L^{2}$
gradient flow) equation, using a second-order accurate, energy-stable
exponential-type scheme. Our exponential-type scheme effectively addresses the
severe stiffness introduced by the biharmonic linear term and achieves second-
order accuracy in a single step. Furthermore, unlike the previous studies [34,
35], our analytical framework leverages the specific energy structure and
employs specialized operator estimations to achieve the desired theoretical
results. This framework is also applicable to the phase field crystal (PFC,
$H^{-1}$ gradient flow) equation [30], which we leave the discussion in a
forthcoming paper.

The SH equation [20] is widely used in the study of phenomena such as
Rayleigh-Bénard convection and more elaborate density functional theories of
liquid interfaces [10, 5, 17, 11, 12]. It differs from the classical AC and CH
systems in that the stable phase is periodic, and is built with the following
Lyapunov energy functional (cf. [20, 6, 14]):

$E(u)=\int_{\Omega}\frac{1}{2}u\underbrace{(\Delta+1)^{2}}_{\mathcal{L}}u+\underbrace{\frac{1}{4}u^{4}-\frac{\varepsilon}{2}u^{2}}_{F(u)}\,\mathrm{d}\mathbf{x},$
(1.5)

which leads to the following fourth-order-in-space system:

$\frac{\partial u}{\partial
t}=-(\Delta+1)^{2}u-f(u),\quad(\mathbf{x},t)\in\Omega\times(0,T],$ (1.6)

where $f(u)=u^{3}-\varepsilon u$, $\Omega\subset\mathbb{R}^{2}$,
$u:\Omega\to\mathbb{R}$ stands for the order parameter, and $\varepsilon$ is a
positive constant with certain physical significance. Similar to the PFC
equation, the SH equation can describe many basic properties of
polycrystalline materials that arise during non-equilibrium processing. The
equation of motion governing these non-equilibrium phenomena is a nonlinear
PDE that generally cannot be solved analytically for random initial
conditions. Therefore, numerous efforts have been made to design appropriate
numerical methods to help researchers understand and characterize non-
equilibrium phenomena, see, e.g., [13, 23, 24, 22, 16, 15] and the references
therein.

A significant challenge is that obtaining energy dissipation using explicit
schemes is very difficult due to severe time-step restrictions. To ensure the
decay of the total energy while employing a moderately large time-step size, a
feasible choice is to use implicit-explicit (IMEX) schemes, where the linear
part is treated implicitly and the nonlinear part is evaluated explicitly.
Motivated by this idea and Duhamel’s formula with respect to (1.1):

$u(t)=\mathrm{e}^{tG}u(t_{0})+\int_{t_{0}}^{t}\mathrm{e}^{G(t-s)}f(u(s))\,\mathrm{d}s,$
(1.7)

where $\mathrm{e}^{tG}$ is the standard semi-group (kernel) with $G$ as the
spatial discretization of $\mathcal{G}$, researchers have widely used energy-
stable exponential-type methods to deal with various gradient flows. For
instance, Ju et al. [31, 32, 33] used first- and second-order exponential time
differencing Runge-Kutta (ETDRK) methods to present numerical analyses for AC
and CH equations, while Chen et al. [40, 41] conducted extensive research on
the application of stabilization ETD multi-step (sETDMS) schemes for thin film
models. Building on these remarkable studies, we shall propose an improved
second-order exponential Runge-Kutta method, namely ERK(2,2) (2.12).
Furthermore, original energy stability of the SH equation (1.6) is rigorously
proved by our proposed method, while multi-step methods only achieve modified
energy stability, which includes a few correction terms that cannot be
removed.

We now state our main results by three steps as follows.
Step 1: by assuming $l^{\infty}$ bounds of the numerical solutions, we are
able to derive original energy stability of the SH equation:

######  Theorem 1.1 (Original energy stability).

The ERK(2,2) scheme (2.12) unconditionally preserves the original energy of
the SH equation; that is,

$E_{N}(u^{n+1})\leq E_{N}(u^{n}),\quad\forall\tau>0,\leavevmode\nobreak\ 0\leq
n\leq N_{t}-1,$

provided that

$\kappa\geq\max\limits_{|\xi|\leq\beta}\frac{|3\xi^{2}-\varepsilon|}{2}\leavevmode\nobreak\
\mbox{with}\leavevmode\nobreak\ \beta:=\max_{0\leq n\leq
N_{t},\leavevmode\nobreak\ i=0,1,2}\|u_{n,i}\|_{\ell^{\infty}}.$ (1.8)

To finish the complete proof, it is necessary to recover the $l^{\infty}$
boundedness of $u^{n+1}$.
Step 2: We suppose energy stable in the first $n$ steps to obtain an upper
bound on the discrete energy, i.e. $E_{N}(u^{n})\leq
E_{N}(u^{n-1})\leq\ldots\leq E_{N}(u^{0})=:C_{e}$ with $C_{e}$ a constant.
Subsequently, with the help of many state-of-art theoretical analysis
techniques, such as the discrete Sobolev inequality, elliptic regularity, as
well as repeated eigenvalue estimates for various Fourier-space operators at
the RK stages, we are able to establish $\ell^{2}$ and $H_{h}^{2}$ bounds of
the numerical solutions at every RK stages.

Step 3: In turn, $\|u^{n+1}\|_{l^{\infty}}$ becomes a direct consequence of an
application of discrete Sobolev embedding. Such an $\ell^{\infty}$ bound
enables us to derive the following desirable result under a certain time-step
size constraint:

######  Theorem 1.2 (Global-in-time energy stability).

With the chosen stabilization parameter in (1.8), which only depends on the
parameter $\varepsilon$, initial energy $E_{0}$ and domain $\Omega$, we select
a time-step size that satisfies the following $\mathcal{O}(1)$ constraint:

$\tau\leq\min\left\\{\frac{1}{256},\frac{1}{64}\tilde{C}_{1}^{-4},(64\kappa)^{-\frac{1}{2}},\frac{1}{4}\tilde{C}_{0}^{-2}\kappa^{-\frac{1}{2}}\right\\},$

where both $\tilde{C}_{1}$ and $\tilde{C}_{0}$ are global-in-time constants.
Then the numerical solution $\\{u^{n}\\}_{0\leq n\leq N_{t}}$ produced by the
ERK(2,2) scheme (2.12) always satisfies $E_{N}(u^{n+1})\leq E_{N}(u^{n})$.

It is observed that with only one condition concerning the time step, the
proof is established without additional assumptions. Moreover, the above two
theorems complement each other, and the constants $\tilde{C}_{1}$ and
$\tilde{C}_{0}$ do not depend on exponential factors involving the time
parameter. Therefore, even if the total time is large, these factors will not
cause instability. This is precisely why we refer to it as “global-in-time
energy stability”.

The remainder of the paper is structured as follows. In Section 2, a fully
discrete numerical scheme is devised, employing the Fourier spectral
collocation spatial discretization and two-stage, second-order exponential
Runge–Kutta temporal integration. The detailed proof of main results is
provided in Section 3, followed by an optimal rate convergence estimate in
Section 4. In Section 5, we present some numerical results to illustrate the
temporal accuracy and long-time dynamic performance of the proposed scheme.
Moreover, some concluding remarks are made in Section 6. The proofs of two
essential propositions in Section 3 are placed in Appendices A and B.

##  2 Fully discretization

To simplify the presentation, a two-dimensional (2-D) domain is assumed, and
an extension to the three-dimensional case could be similarly handled without
an essential difficulty. Furthermore, we assume throughout the paper that
periodic boundary conditions are chosen such that all boundary terms will
vanish when integration by parts is performed; of course, an extension to the
homogeneous Neumann boundary condition case is straightforward.

###  2.1 Review of the Fourier pseudo-spectral approximation

We assume that the domain is given by $\Omega=[0,L]^{2}$, with a uniform mesh
size: $N_{x}=N_{y}=N$, $N\cdot h=L$. The number of grid points (in each
direction) is set as $N=2K+1$, and the case for an even $N$ could be similarly
treated. All the spatial variables are evaluated on the $N\times N$ uniform
mesh $\Omega_{N}$, in which $x_{p}=ph,\leavevmode\nobreak\ y_{q}=qh$, $0\leq
p,q\leq N-1$, $h=\frac{L}{N}$. We also denote a uniform time step size
$\tau=\frac{T}{N_{t}}$, where $N_{t}$ is a positive integer, and $t_{n}=n\tau$
for $0\leq n\leq N_{t}$. Let $\mathcal{M}_{N}$ denote the set of 2-D periodic
grid functions defined on $\Omega_{N}$. For any gird functions
${f},{g}\in\mathcal{M}_{N}$, the discrete $\ell^{2}$ norm and inner product,
and discrete $\ell^{\infty}$ norm are introduced as

$\|{f}\|_{2}:=\sqrt{\langle
f,f\rangle}\quad\text{with}\quad\langle{f},{g}\rangle:=h^{2}\sum_{p,q=0}^{N-1}f_{p,q}g_{p,q};\quad\|{f}\|_{\infty}:=\max\limits_{p,q}|f_{p,q}|.$

For ${f}\in\mathcal{M}_{N}$, we set the discrete Fourier expansion as

$f_{p,q}=\sum_{\ell,m=-K}^{K}\hat{f}_{\ell,m}\mathrm{e}^{2\pi\mathrm{i}(\ell
x_{p}+my_{q})/L},\quad\text{with}\leavevmode\nobreak\ \mathrm{i}=\sqrt{-1},$
(2.1)

where the coefficients $\hat{f}_{\ell,m}$ are obtained by the discrete Fourier
transform.

Although no aliasing error needs to be considered in the numerical analysis,
due to the lack of spatial derivative terms in the nonlinear parts of the SH
equation (1.6), we have to introduce a periodic extension of a grid function
and a Fourier collocation interpolation operator to facilitate the later
analysis.

######  Definition 2.1.

For any $f\in\mathcal{M}_{N}$, we denote its periodic continuous extension
into $\mathcal{B}^{K}$ (the space of trigonometric polynomials of degree at
most $K$) as $f_{S}$, given by

$f_{S}(x,y)=S_{N}(f)=\sum_{\ell,m=-K}^{K}\hat{f}_{\ell,m}\mathrm{e}^{2\pi\mathrm{i}(\ell
x+my)/L}.$

We call $S_{N}:\mathcal{M}_{N}\to\mathcal{B}^{K}$ the spectral interpolation
operator.

######  Definition 2.2.

The discrete differentiation operators in the $x$-direction are defined as

$\displaystyle(\mathcal{D}_{x}f)_{p,q}:=\sum_{\ell,m=-K}^{K}\frac{2\ell\pi\mathrm{i}}{L}\hat{f}_{\ell,m}\mathrm{e}^{2\pi\mathrm{i}(\ell
x_{p}+my_{q})/L},$
$\displaystyle(\mathcal{D}_{x}^{2}f)_{p,q}:=\sum_{\ell,m=-K}^{K}\frac{-4\ell^{2}\pi^{2}}{L^{2}}\hat{f}_{\ell,m}\mathrm{e}^{2\pi\mathrm{i}(\ell
x_{p}+my_{q})/L}.$

The differentiation operators in the $y$-direction, $\mathcal{D}_{y}$ and
$\mathcal{D}_{y}^{2}$, can be introduced in the same fashion. In turn, the
discrete Laplacian, gradient and divergence operators are given by

$\Delta_{N}{f}:=\mathcal{D}_{x}^{2}{f}+\mathcal{D}_{y}^{2}{f},\quad\nabla_{N}{f}:=\begin{pmatrix}\mathcal{D}_{x}{f},\\\
\mathcal{D}_{y}{f}\end{pmatrix},\quad\nabla_{N}\cdot\begin{pmatrix}f_{1}\\\
f_{2}\end{pmatrix}:=\mathcal{D}_{x}f_{1}+\mathcal{D}_{y}f_{2},$

at the point-wise level.

Detailed calculations reveal that the following summation-by-parts formulas
(cf. [2, 3, 7, 8]) are valid .

######  Proposition 2.1.

For any periodic grid functions ${f},{g}\in\mathcal{M}_{N}$, we have

$\langle{f},\Delta_{N}{g}\rangle=-\langle\nabla_{N}{f},\nabla_{N}{g}\rangle=\langle\Delta_{N}{f},{g}\rangle,\quad\langle{f},\nabla_{N}\cdot{g}\rangle=-\langle\nabla_{N}{f},{g}\rangle.$

###  2.2 Time-stepping integrator

The space-discrete problem of equation (1.6) turns out to be: find
${u}:[0,+\infty)\to\mathcal{M}_{N}$ which satisfies

$\frac{\mathrm{d}u}{\mathrm{d}t}=-(\Delta_{N}+I)^{2}u-f(u).$ (2.2)

Meanwhile, the discrete energy is given by

$E_{N}(u)=\frac{1}{2}\|(\Delta_{N}+I)u\|_{2}^{2}+\langle\frac{1}{4}u^{4}-\frac{\varepsilon}{2}u^{2},1\rangle.$

Adding and subtracting the stabilization term $\kappa{u}$ ($\kappa>0$) on the
right-hand side (RHS) of equation (2.2), we obtain the following ODE system:

$\left\\{\begin{aligned}
&\frac{\mathrm{d}{u}}{\mathrm{d}t}=-L_{\kappa}{u}+N_{\kappa}(u),\quad
t\in(0,T],\\\ &{u}(0)=u_{0},\qquad\qquad\quad\leavevmode\nobreak\
\leavevmode\nobreak\ t=0,\end{aligned}\right.$

where $u_{0}\in\mathcal{M}_{N}$ is given by the initial data,
$L_{\kappa}=(\Delta_{N}+I)^{2}+\kappa I$, and $N_{\kappa}(u)=\kappa u-f(u)$.
Let $\\{\mathrm{e}^{-tL_{\kappa}}\\}_{t\geq 0}$ denote the semigroup on
$\Omega_{N}$ with the generator $(-L_{\kappa})$. Based on the integrating
factor, an update of the exact solution from time instant $t_{n}$ to the next
time step $t_{n+1}$ could be expressed as

${u}(t_{n+1})=\mathrm{e}^{-\tau
L_{\kappa}}{u}(t_{n})+\int_{0}^{\tau}\mathrm{e}^{-(\tau-r)L_{\kappa}}N_{\kappa}({u}(t_{n}+r))\,\mathrm{d}r.$
(2.3)

Denoting by $u^{n}$ the fully discrete numerical solution at the time step
$t_{n}$, a family of second-order ERK approach derived by Hochbruck and
Ostermann [9] for solving equation (2.3) is formulated as below: for
$n=0,1,\ldots,N_{t}-1$,

$u_{n,i}=\varphi_{0}(\tau L_{\kappa})u^{n}+\tau\sum_{j=0}^{i-1}a_{i,j}(\tau
L_{\kappa})N_{\kappa}(u_{n,j}),\quad i=1,2,$ (2.4)

in which the coefficients $a_{i,j}(z)$ that satisfy the second-order
conditions in Section 5.1 of [9] are linear combinations of the
$\varphi_{k}(z)$ functions defined by

$\varphi_{0}(z)=\mathrm{e}^{-z},\quad\varphi_{1}(z)=\frac{1-\mathrm{e}^{-z}}{z},\quad\varphi_{2}(z)=\frac{\mathrm{e}^{-z}-1+z}{z^{2}},\leavevmode\nobreak\
\text{with}\leavevmode\nobreak\ \varphi_{k}(0)=\frac{1}{k!}.$

It is observed that $a_{i,j}(z)$ are specified in the following one-parameter
Butcher-like tableau (cf. equation (5.8) in [9]):

$\displaystyle\begin{array}[]{>{\centering\arraybackslash$} p{0.4cm} <{$} | >{\centering\arraybackslash$} p{1.2cm} <{$} >{\centering\arraybackslash$} p{0.8cm} <{$}}c_{0}$\@add@centering&0$\@add@centering&$\@add@centering\\\ c_{1}$\@add@centering&a_{1,0}(z)$\@add@centering&0$\@add@centering\cr\hline\cr c_{2}$\@add@centering&a_{2,0}(z)$\@add@centering&a_{2,1}(z)$\@add@centering\end{array}\leavevmode\nobreak\ \leavevmode\nobreak\ =\leavevmode\nobreak\ \leavevmode\nobreak\ \begin{array}[]{>{\centering\arraybackslash$} p{0.4cm} <{$} | >{\centering\arraybackslash$} p{1.8cm} <{$} >{\centering\arraybackslash$} p{0.8cm} <{$}}0$\@add@centering&0$\@add@centering&$\@add@centering\\\ c_{1}$\@add@centering&c_{1}\varphi_{1,1}$\@add@centering&0$\@add@centering\cr\hline\cr 1$\@add@centering&(1-\frac{1}{2c_{1}})\varphi_{1}$\@add@centering&\frac{1}{2c_{1}}\varphi_{1}$\@add@centering\end{array},$ (2.11)

where $\varphi_{1,1}(z):=\varphi_{1}(c_{1}z)$. Notice that the choice
$c_{1}=\frac{1}{2}$ yields $a_{2,0}=0$. We apply this particular instance to
derive a two-stage, second-order fully discrete scheme, denoted as ERK(2,2):

$\displaystyle{u}^{n+1}$ $\displaystyle=\varphi_{0}(\tau
L_{\kappa}){u}^{n}+\tau\varphi_{1}(\tau L_{\kappa})N_{\kappa}(u_{n,1})$ (2.12)
$\displaystyle=P^{-1}\left[\varphi_{0}(\tau\widehat{L}_{\kappa})Pu^{n}+\tau\varphi_{1}(\tau\widehat{L}_{\kappa})PN_{\kappa}(u_{n,1})\right],$
$\displaystyle u_{n,1}$ $\displaystyle=\varphi_{0}(\frac{1}{2}\tau
L_{\kappa}){u}^{n}+\frac{\tau}{2}\varphi_{1}(\frac{1}{2}\tau
L_{\kappa})N_{\kappa}({u}^{n})$
$\displaystyle=P^{-1}\left[\varphi_{0}(\frac{1}{2}\tau\widehat{L}_{\kappa})Pu^{n}+\frac{\tau}{2}\varphi_{1}(\frac{1}{2}\tau\widehat{L}_{\kappa})PN_{\kappa}({u}^{n})\right],$

where the operators $P$ and $P^{-1}$ can be implemented by the 2-D fast
Fourier transform (FFT) and the corresponding inverse transform, respectively.
Therefore, the overall computational complexity is
$\mathcal{O}(N^{2}\log_{2}N)$ per time step. In addition,
$\widehat{L}_{\kappa}=PL_{\kappa}P^{-1}$.

######  Remark 2.1.

Note that another second-order typical ERK scheme, ETDRK2 (C.7), shares the
same stage as ERK(2,2) but possesses a more complex structure. In a previous
work [18], the authors conducted a numerical investigation revealing that
ETDRK2 achieves better accuracy than the scheme (2.11) with $c_{1}=1$.
Nonetheless, when compared with ETDRK2, ERK(2,2) demonstrates both enhanced
accuracy and faster computation (cf. Figure 1 in Section 5), representing a
significant computational advantage.

######  Remark 2.2.

In fact, there is one more difference between the two methods; that is, ETDRK2
is of stiff order two, while ERK(2,2) has stiff convergence order two [39].
Existing numerical evidence suggests that the latter method demonstrates
superior error accuracy due to its second-order global error, whereas the
former exhibits enhanced stability when addressing stiff problems. Remaining
this observation in actual numerical implementation remains to be studied
further.

##  3 Proof of Main results

Now we turn to the energy stability analysis of the ERK(2,2) scheme. For any
${u}\in\mathcal{M}_{N}$, the discrete energy functional could be rewritten as
$E_{N}(u)=E_{N,c}(u)+E_{N,e}(u)$, with

$E_{N,c}({u})=\frac{1}{2}\langle{u},(\Delta_{N}+I)^{2}{u}\rangle,\quad
E_{N,e}({u})=\langle\frac{1}{4}{u}^{4}-\frac{\varepsilon}{2}{u}^{2},{1}\rangle.$

###  3.1 Proof of 1.1

We first prove that the ERK(2,2) scheme (2.12), with a sufficiently large
$\kappa$, is original energy stable.

###### Proof.

A difference of the two energy functional $E_{N,c}({v})$ and $E_{N,c}({w})$
yields

$\displaystyle
E_{N,c}({v})-E_{N,c}({w})=\frac{1}{2}\langle{v},(\Delta_{N}+I)^{2}{v}\rangle-\frac{1}{2}\langle{w},(\Delta_{N}+I)^{2}{w}\rangle$
(3.1)
$\displaystyle=\langle{v}-{w},(\Delta_{N}+I)^{2}{v}\rangle-\frac{1}{2}\langle{v}-{w},(\Delta_{N}+I)^{2}({v}-{w})\rangle$
$\displaystyle\leq\langle{v}-{w},L_{\kappa}{v}\rangle-\kappa\langle{v}-{w},{v}\rangle.$

For $E_{N,e}({v})$ and $E_{N,e}({w})$, a careful application of Taylor’s
expansion indicates that

$\displaystyle
E_{N,e}({v})-E_{N,e}({w})=\langle{v}-{w},f({w})\rangle+\frac{1}{2}\langle{v}-{w},f{{}^{\prime}}({\xi})({v}-{w})\rangle$
(3.2)
$\displaystyle=-\langle{v}-{w},N_{\kappa}({w})\rangle+\kappa\langle{v}-{w},{w}\rangle+\frac{1}{2}\langle{v}-{w},f{{}^{\prime}}({\xi})({v}-{w})\rangle,$

in which the variable ${\xi}$ is between ${v}$ and ${w}$, at a point-wise
level.

A combination of (3.1) and (3.2) leads to

$\displaystyle
E_{N}({v})-E_{N}({w})\leq\langle{v}-{w},L_{\kappa}{v}-\kappa{v}\rangle+\langle{v}-{w},-N_{\kappa}({w})+\kappa{w}\rangle+\frac{1}{2}\langle{v}-{w},f{{}^{\prime}}({\xi})({v}-{w})\rangle$
$\displaystyle=\langle{v}-{w},L_{\kappa}{v}-N_{\kappa}({w})\rangle+\langle{v}-{w},[\frac{1}{2}f{{}^{\prime}}({\xi})-\kappa
I]({v}-{w})\rangle.$

In turn, under the condition that
$\kappa\geq\max\limits_{|\xi|\leq\beta}\frac{|3\xi^{2}-\varepsilon|}{2}\geq\frac{1}{2}f^{\prime}(\xi)$,
the following inequality is available:

$E_{N}({v})-E_{N}({w})\leq\langle{v}-{w},L_{\kappa}{v}-N_{\kappa}({w})\rangle.$

Then we arrive at

$\displaystyle E_{N}(u_{n,1})-E_{N}({u}^{n})\leq\langle
u_{n,1}-{u}^{n},L_{\kappa}u_{n,1}-N_{\kappa}({u}^{n})\rangle$ (3.3)
$\displaystyle=\langle
u_{n,1}-{u}^{n},L_{\kappa}u_{n,1}-[c_{1}\tau\varphi_{1}(c_{1}\tau
L_{\kappa})]^{-1}[u_{n,1}-u^{n}+c_{1}\tau L_{\kappa}\varphi_{1}(c_{1}\tau
L_{\kappa}){u}^{n}]\rangle$ $\displaystyle=\langle
u_{n,1}-{u}^{n},\\{L_{\kappa}-\left[c_{1}\tau\varphi_{1}\left(c_{1}\tau
L_{\kappa}\right)\right]^{-1}\\}(u_{n,1}-{u}^{n})\rangle$
$\displaystyle=:\langle u_{n,1}-{u}^{n},\Delta_{1}(u_{n,1}-{u}^{n})\rangle,$

where $\Delta_{1}=2\tau^{-1}h_{1}(c_{1}\tau L_{\kappa})$, with
$h_{1}(z)=z-[\varphi_{1}(z)]^{-1}={z\mathrm{e}^{-z}}/{(\mathrm{e}^{-z}-1)}$.
Notice that $L_{\kappa}$ is symmetric positive definite, thus $h_{1}(z)<0$ for
any $z\neq 0$. Consequently, the operator $\Delta_{1}$ is negative semi-
definite and $E_{N}(u_{n,1})-E_{N}({u}^{n})\leq 0$.

As for the second step, a difference between $E_{N}({u}^{n+1})$ and
$E_{N}(u_{n,1})$ reveals that

$\displaystyle
E_{N}({u}^{n+1})-E_{N}(u_{n,1})\leq\langle{u}^{n+1}-u_{n,1},L_{\kappa}{u}^{n+1}-N_{\kappa}(u_{n,1})\rangle$
(3.4)
$\displaystyle=\langle{u}^{n+1}-u_{n,1},L_{\kappa}{u}^{n+1}-\left[\tau\varphi_{1}(\tau
L_{\kappa})\right]^{-1}[{u}^{n+1}-{u}^{n}+\tau L_{\kappa}\varphi_{1}(\tau
L_{\kappa}){u}^{n}]\rangle$
$\displaystyle=\langle{u}^{n+1}-u_{n,1},\\{L_{\kappa}-[\tau\varphi_{1}(\tau
L_{\kappa})]^{-1}\\}({u}^{n+1}-{u}^{n})\rangle$
$\displaystyle=:\langle{u}^{n+1}-u_{n,1},\Delta_{2}({u}^{n+1}-{u}^{n})\rangle,$

with $\Delta_{2}=\tau^{-1}h_{1}(\tau L_{\kappa})$ a negative semi-definite
operator. In turn, we obtain $E_{N}({u}^{n+1})-E_{N}(u_{n,1})\leq 0$.
Subsequently, a summation of inequalities (3.3) and (3.4) yields

$\displaystyle E_{N}({u}^{n+1})-E_{N}({u}^{n})\leq\langle
u_{n,1}-{u}^{n},\Delta_{1}(u_{n,1}-{u}^{n})\rangle+\langle{u}^{n+1}-u_{n,1},\Delta_{2}({u}^{n+1}-{u}^{n})\rangle$
(3.5) $\displaystyle=\underbrace{\langle
u_{n,1}-{u}^{n},\Delta_{1}(u_{n,1}-{u}^{n})\rangle}_{A}+\underbrace{\langle{u}^{n+1}-u_{n,1},\Delta_{2}({u}^{n+1}-u_{n,1})\rangle}_{B}$
$\displaystyle+\underbrace{\langle{u}^{n+1}-u_{n,1},\Delta_{2}(u_{n,1}-{u}^{n})\rangle}_{C},$
$A=\underbrace{\langle
u_{n,1}-{u}^{n},(\Delta_{1}-\frac{1}{2}\Delta_{2})(u_{n,1}-{u}^{n})\rangle}_{A_{1}}+\underbrace{\frac{1}{2}\langle
u_{n,1}-{u}^{n},\Delta_{2}(u_{n,1}-{u}^{n})\rangle}_{A_{2}},$

where $\Delta_{1}-\frac{1}{2}\Delta_{2}=2\tau^{-1}h_{2}(\tau L_{\kappa})$,
with
$h_{2}(z)=h_{1}(\frac{z}{2})-\frac{1}{4}h_{1}(z)=\frac{z\mathrm{e}^{-z/2}}{2\left(\mathrm{e}^{-z/2}-1\right)}-\frac{z\mathrm{e}^{-z}}{4\left(\mathrm{e}^{-z}-1\right)}$.
It can be verified that $h_{2}\leq 0$ for any $z\neq 0$. Therefore,
$\Delta_{1}-\frac{1}{2}\Delta_{2}$ is symmetric negative semi-definite and
$A_{1}\leq 0$. We also notice that

$\displaystyle
A_{2}+\frac{1}{2}C=\frac{1}{2}\langle{u}^{n+1}-{u}^{n},\Delta_{2}(u_{n,1}-{u}^{n})\rangle,\quad\frac{1}{2}(B+C)=\frac{1}{2}\langle{u}^{n+1}-u_{n,1},\Delta_{2}({u}^{n+1}-{u}^{n})\rangle,$

which in turn leads to

$A_{2}+\frac{1}{2}B+C=\frac{1}{2}\langle{u}^{n+1}-{u}^{n},\Delta_{2}({u}^{n+1}-{u}^{n})\rangle\leq
0.$

As a result, inequality (3.5) turns out to be

$E_{N}({u}^{n+1})-E_{N}({u}^{n})\leq
A_{1}+(A_{2}+\frac{1}{2}B+C)+\frac{1}{2}B\leq 0.$

This completes the proof. ∎

###  3.2 Proof of 1.2

1.1 indicates that, the choice of stabilization parameter in (1.8) implicitly
uses an a-priori $\ell^{\infty}$ bound assumption on $u_{n,i}$ in order to
make $\kappa$ a controllable constant. Of course, it is desirable to remove
this technical restriction and establish a more reasonable energy stability
theory.

In this subsection, we perform a direct analysis for the numerical solution of
the ERK(2,2) scheme (2.12), so that uniform-in-time $\ell^{2}$ and $H_{h}^{2}$
estimates become available for the numerical solutions at all the stages. With
the help of the discrete Sobolev embedding, we are able to recover the global-
in-time values for $\beta$ and $\kappa$ in (1.8), which allows us to derive a
global-in-time energy stability estimate for the ERK(2,2) scheme.

####  3.2.1 A few preliminary notations and results

The verification of the following Calculus-style analysis is straightforward.

######  Lemma 3.1.

For the ERK(2,2) scheme (2.12), the functions $\varphi_{i}(c_{1}z)$ are
decreasing and $\varphi_{i}(c_{1}z)\leq 1$, $\forall z>0$, where $i=0,1$.

We introduce the following linear operators to facilitate the energy stability
analysis:

$G_{i}=c_{i}\varphi_{1}(c_{i}\tau L_{\kappa})=(\tau
L_{\kappa})^{-1}(I-\mathrm{e}^{-c_{i}\tau L_{\kappa}}),\quad i=1,2.$

In more details, for ${f}\in\mathcal{M}_{N}$ with a discrete Fourier expansion
as (2.1), an application of of $G_{i}$ becomes

$\left(G_{i}f\right)_{p,q}=\sum_{\ell,m=-K}^{K}c_{i}\varphi_{1}(c_{i}\tau\Lambda_{\ell,m})\hat{f}_{\ell,m}\mathrm{e}^{2\pi\mathrm{i}(\ell
x_{p}+my_{q})/L},\leavevmode\nobreak\ \text{with}\leavevmode\nobreak\
\Lambda_{\ell,m}=(-1+\lambda_{\ell,m})^{2}+\kappa,$

where $\lambda_{\ell,m}=\frac{4\pi^{2}}{L^{2}}(\ell^{2}+m^{2})$. Since all the
eigenvalues of $G_{i}$ are non-negative, we define $G_{i}^{\frac{1}{2}}$ and
$G_{1,2}^{\frac{1}{2}}=G_{1}^{\frac{1}{2}}G_{2}^{\frac{1}{2}}$ as

$\displaystyle(G_{i}^{\frac{1}{2}}f)_{p,q}=\sum_{\ell,m=-K}^{K}\left(c_{i}\varphi_{1}(c_{i}\tau\Lambda_{\ell,m})\right)^{\frac{1}{2}}\hat{f}_{\ell,m}\mathrm{e}^{2\pi\mathrm{i}(\ell
x_{p}+my_{q})/L};$
$\displaystyle(G_{1,2}^{\frac{1}{2}}f)_{p,q}=\sum_{\ell,m=-K}^{K}\left(c_{1}c_{2}\varphi_{1}(c_{1}\tau\Lambda_{\ell,m})\varphi_{1}(c_{2}\tau\Lambda_{\ell,m})\right)^{\frac{1}{2}}\hat{f}_{\ell,m}\mathrm{e}^{2\pi\mathrm{i}(\ell
x_{p}+my_{q})/L}.$

It is obvious that the operators $G_{i}^{\frac{1}{2}}$ and
$G_{1,2}^{\frac{1}{2}}$ are commutative with any differential operators in the
Fourier pseudo-spectral space, and the summation by parts formulas are
available:

$\langle{f},G_{i}{g}\rangle=\langle
G_{i}^{\frac{1}{2}}{f},G_{i}^{\frac{1}{2}}{g}\rangle,\quad\langle
G_{1}{f},G_{2}{g}\rangle=\langle
G_{1,2}^{\frac{1}{2}}{f},G_{1,2}^{\frac{1}{2}}{g}\rangle.$

In addition, the operators $G_{i}^{*}$, $G_{i}^{**}$, $G^{*}$, and $G^{**}$
are introduced to facilitate the analysis for the diffusion part:

$\displaystyle\left(G_{i}^{*}f\right)_{p,q}=\sum_{\ell,m=-K}^{K}\left(c_{i}\Lambda_{\ell,m}\varphi_{1}(c_{i}\tau\Lambda_{\ell,m})\right)^{\frac{1}{2}}\hat{f}_{\ell,m}\mathrm{e}^{2\pi\mathrm{i}(\ell
x_{p}+my_{q})/L};$
$\displaystyle\left(G_{i}^{**}f\right)_{p,q}=\sum_{\ell,m=-K}^{K}\left(c_{i}\Lambda_{\ell,m}\varphi_{1}(c_{i}\tau\Lambda_{\ell,m})\right)^{\frac{1}{2}}\lambda_{\ell,m}\hat{f}_{\ell,m}\mathrm{e}^{2\pi\mathrm{i}(\ell
x_{p}+my_{q})/L};$
$\displaystyle\left(G^{*}f\right)_{p,q}=\sum_{\ell,m=-K}^{K}\left(c_{1}c_{2}\Lambda_{\ell,m}^{2}\varphi_{1}(c_{1}\tau\Lambda_{\ell,m})\varphi_{1}(c_{2}\tau\Lambda_{\ell,m})\right)^{\frac{1}{2}}\hat{f}_{\ell,m}\mathrm{e}^{2\pi\mathrm{i}(\ell
x_{p}+my_{q})/L};$
$\displaystyle\left(G^{**}f\right)_{p,q}=\sum_{\ell,m=-K}^{K}\left(c_{1}c_{2}\Lambda_{\ell,m}^{2}\varphi_{1}(c_{1}\tau\Lambda_{\ell,m})\varphi_{1}(c_{2}\tau\Lambda_{\ell,m})\right)^{\frac{1}{2}}\lambda_{\ell,m}\hat{f}_{\ell,m}\mathrm{e}^{2\pi\mathrm{i}(\ell
x_{p}+my_{q})/L}.$

The following identities could be verified in a straightforward way:

$\langle
G_{i}L_{\kappa}{f},(I+\Delta_{N}^{2}){f}\rangle=\|G_{i}^{*}{f}\|_{2}^{2}+\|G_{i}^{**}f\|_{2}^{2};\leavevmode\nobreak\
\langle
G_{1}L_{\kappa}{f},G_{2}(I+\Delta_{N}^{2}){f}\rangle=\|G^{*}{f}\|_{2}^{2}+\|G^{**}f\|_{2}^{2}.$

By the fact that $\varphi_{1}(c_{1}z)\leq 1$ (by 3.1), it is easy to derive
3.1. We also introduce other two propositions whose proofs are postponed to
Appendices A and B.

######  Proposition 3.1.

For any ${f}\in\mathcal{M}_{N}$, the following two estimates are valid:

$\|G_{i}^{\frac{1}{2}}f\|_{2}\leq\|f\|_{2},\qquad\|G_{1,2}^{\frac{1}{2}}f\|_{2}\leq\|f\|_{2}.$
(3.6)

######  Proposition 3.2.

For $\kappa\geq 1$ and $f\in\mathcal{M}_{N}$, we have

$\displaystyle\|G_{i}^{*}{f}\|_{2}^{2}+\|G_{i}^{**}{f}\|_{2}^{2}\geq\frac{1}{4}(\|G_{i}^{\frac{1}{2}}\Delta_{N}{f}\|_{2}^{2}+\|G_{i}^{\frac{1}{2}}\Delta_{N}^{2}{f}\|_{2}^{2})$
(3.7)
$\displaystyle+(\kappa-1)(\|G_{i}^{\frac{1}{2}}{f}\|_{2}^{2}+\|G_{i}^{\frac{1}{2}}\Delta_{N}{f}\|_{2}^{2})+\frac{2}{3}(\|G_{i}^{\frac{1}{2}}{f}\|_{2}^{2}+\|G_{i}^{\frac{1}{2}}\Delta_{N}{f}\|_{2}^{2}).$
$\displaystyle\|G^{*}{f}\|_{2}^{2}+\|G^{**}{f}\|_{2}^{2}\geq\frac{1}{4}(\|G_{1,2}^{\frac{1}{2}}\Delta_{N}{f}\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}^{2}{f}\|_{2}^{2})$
(3.8)
$\displaystyle+(\kappa-1)(\|G_{1,2}^{\frac{1}{2}}{f}\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}{f}\|_{2}^{2})+\frac{2}{3}(\|G_{1,2}^{\frac{1}{2}}f\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}f\|_{2}^{2}).$

######  Proposition 3.3.

For $f,g\in\mathcal{M}_{N}$, we have

$\displaystyle\tau\langle G_{i}L_{\kappa}{f},\varphi_{0}(c_{i}\tau
L_{\kappa})f\rangle+\|{g}-\varphi_{0}(c_{i}\tau
L_{\kappa})f\|_{2}^{2}\geq\tau\|G_{i}^{*}{g}\|_{2}^{2};$ (3.9)
$\displaystyle\tau\langle
G_{i}L_{\kappa}{f},\Delta_{N}^{2}\varphi_{0}(c_{i}\tau
L_{\kappa}){f}\rangle+\|\Delta_{N}({g}-\varphi_{0}(c_{i}\tau
L_{\kappa}){f})\|_{2}^{2}\geq\tau\|G_{i}^{**}{g}\|_{2}^{2};$ (3.10)
$\displaystyle\tau\langle G_{1}L_{\kappa}{f},G_{2}\varphi_{0}(c_{i}\tau
L_{\kappa}){f}\rangle+\|G_{2}^{\frac{1}{2}}(g-\varphi_{0}(c_{i}\tau
L_{\kappa})f)\|_{2}^{2}\geq\tau\|G^{*}{g}\|_{2}^{2};$ (3.11)
$\displaystyle\tau\langle
G_{1}L_{\kappa}{f},G_{2}\Delta_{N}^{2}\varphi_{0}(c_{i}\tau
L_{\kappa}){f}\rangle+\|G_{2}^{\frac{1}{2}}\Delta_{N}({g}-\varphi_{0}(c_{i}\tau
L_{\kappa}){f})\|_{2}^{2}\geq\tau\|G^{**}{g}\|_{2}^{2}.$ (3.12)

######  Lemma 3.2.

For $f\in\mathcal{M}_{N}$, it follows that

$\|{f}\|_{\infty}\leq\hat{C}\left(\|{f}\|_{2}+\|\Delta_{N}{f}\|_{2}\right),$

where the constant $\hat{C}$ is independent of ${f}$, $h$, and $\kappa$.

###### Proof.

For any periodic function $f\in\mathcal{M}_{N}$, we recall its continuous
extension, $f_{S}=S_{N}(f)\in\mathcal{B}^{K}$, as introduced in 2.1. Since $f$
is the point-wise interpolation of $f_{S}$, we see that
$\|f\|_{\infty}\leq\|f_{S}\|_{L^{\infty}}$. For any smooth function $f_{S}$,
applying the 2-D Sobolev inequality associated the embedding
$H^{2}\hookrightarrow L^{\infty}$ and the elliptic regularity, it holds that

$\|f_{S}\|_{L^{\infty}}\leq\hat{C}(\|f_{S}\|_{L^{2}}+\|\Delta
f_{S}\|_{L^{2}})=\hat{C}(\|f\|_{2}+\|\Delta_{N}f\|_{2}).$

∎

To proceed the energy analysis, we make an a-priori assumption at the previous
time step:

$E_{N}({u}^{n})\leq E_{N}({u}^{0})=:C_{e}.$ (3.13)

Such an assumption will be recovered at the next time step. Afterwards, the
$\ell^{2}$ and $H_{h}^{2}$ bounds for the numerical solutions $u^{n}$ could be
derived (cf. Lemma 3.7 of [21]),

$\|{u}^{n}\|_{2},\leavevmode\nobreak\ \|\Delta_{N}{u}^{n}\|_{2}\leq
C_{0}:=2\sqrt{C_{e}+|\Omega|}.$ (3.14)

Based on an application of 3.2, we have a discrete $\ell^{\infty}$ bound at
time step $t^{n}$:

$\|u^{n}\|_{\infty}\leq\hat{C}(\|u^{n}\|_{2}+\|\Delta_{N}u^{n}\|_{2})\leq
2\hat{C}C_{0}=:\tilde{C}_{0}.$ (3.15)

####  3.2.2 Preliminary $\ell^{2}$ and $H_{h}^{2}$ estimates of $u_{n,1}$

To obtain a rough $\ell^{\infty}$ estimate for the intermediate-stage
numerical solution, we have to derive the $\ell^{2}$ and $H_{h}^{2}$ estimates
for $u_{n,1}$, given by (2.12). Meanwhile, an intuitional interaction between
the linear and nonlinear terms is not clearly presented in the current
numerical formulation (2.12). To remedy this issue, by introducing
$u_{n,i}^{*}=\varphi_{0}\left(c_{i}\tau L_{\kappa}\right){u}^{n}$, the
algorithm can be recast as a two-substage system, facilitating a more
convenient theoretical analysis:

$\displaystyle\frac{{u}_{n,i}^{*}-{u}^{n}}{\tau}$
$\displaystyle=-c_{i}L_{\kappa}\varphi_{1}(c_{i}\tau L_{\kappa}){u}^{n},$
(3.16) $\displaystyle\frac{u_{n,i}-{u}_{n,i}^{*}}{\tau}$
$\displaystyle=c_{i}\varphi_{1}(c_{i}\tau L_{\kappa})N_{\kappa}(u_{n,j}),\quad
i=1,2,\leavevmode\nobreak\ j=i-1.$ (3.17)

Taking a discrete $\ell^{2}$ inner product with equality (3.16) by
$(I+\Delta_{N}^{2})(u_{n,i}^{*}+u^{n})$, combined with the summation by parts
formula, yields

$(I+\Delta_{N}^{2})(\|{u}_{n,i}^{*}\|_{2}^{2}-\|{u}^{n}\|_{2}^{2})+\tau(\|G_{i}^{*}{u}^{n}\|_{2}^{2}+\|G_{i}^{**}{u}^{n}\|_{2}^{2}+\langle
G_{i}L_{\kappa}{u}^{n},(I+\Delta_{N}^{2}){u}_{n,i}^{*}\rangle)=0.$ (3.18)

Taking a discrete $\ell^{2}$ inner product with equality (3.17) by
$2(I+\Delta_{N}^{2})u_{n,i}$ leads to

$\langle u_{n,i}-{u}_{n,i}^{*},2(I+\Delta_{N}^{2})u_{n,i}\rangle=2\tau\langle
G_{i}N_{\kappa}(u_{n,j}),(I+\Delta_{N}^{2})u_{n,i}\rangle.$ (3.19)

The term on the left-hand side (LHS) of equality (3.19) can be rewritten as:

$\displaystyle\langle
u_{n,i}-{u}_{n,i}^{*},2(I+\Delta_{N}^{2})u_{n,i}\rangle=\|u_{n,i}\|_{2}^{2}-\|u_{n,i}^{*}\|_{2}^{2}+\|u_{n,i}-{u}_{n,i}^{*}\|_{2}^{2}+\|\Delta_{N}u_{n,i}\|_{2}^{2}$
(3.20)
$\displaystyle-\|\Delta_{N}u_{n,i}^{*}\|_{2}^{2}+\|\Delta_{N}(u_{n,i}-{u}_{n,i}^{*})\|_{2}^{2},$

where the identity $a^{2}-b^{2}=2a(a-b)-(a-b)^{2}$ has been employed. In turn,
a combination of equalities (3.18)-(3.20) leads to

$\displaystyle\|u_{n,i}\|_{2}^{2}-\|u^{n}\|_{2}^{2}+\|u_{n,i}-{u}_{n,i}^{*}\|_{2}^{2}+\|\Delta_{N}u_{n,i}\|_{2}^{2}-\|\Delta_{N}u_{n,i}^{*}\|_{2}^{2}+\|\Delta_{N}(u_{n,i}-{u}_{n,i}^{*})\|_{2}^{2}$
(3.21)
$\displaystyle+\tau(\|G_{i}^{*}{u}^{n}\|_{2}^{2}+\|G_{i}^{**}{u}^{n}\|_{2}^{2}+\langle
G_{i}L_{\kappa}{u}^{n},(I+\Delta_{N}^{2}){u}_{n,i}^{*}\rangle)=2\tau\langle
G_{i}N_{\kappa}(u_{n,j}),(I+\Delta_{N}^{2})u_{n,i}\rangle.$

Meanwhile, by inequalities (3.9)-(3.10) in 3.3, we see that

$\displaystyle\|u_{n,i}-u_{n,i}^{*}\|_{2}^{2}+\|\Delta_{N}(u_{n,i}-u_{n,i}^{*})\|_{2}^{2}+\tau\langle
G_{i}L_{\kappa}{u}^{n},(I+\Delta_{N}^{2})u_{n,i}^{*}\rangle$
$\displaystyle\geq\tau(\|G_{i}^{*}u_{n,i}\|_{2}^{2}+\|G_{i}^{**}u_{n,i}\|_{2}^{2}).$

Going back to equality (3.21), we arrive at

$\displaystyle\|u_{n,i}\|_{2}^{2}-\|u^{n}\|_{2}^{2}+\|\Delta_{N}u_{n,i}\|_{2}^{2}-\|\Delta_{N}u^{n}\|_{2}^{2}+\tau(\|G_{i}^{*}u_{n,i}\|_{2}^{2}+\|G_{i}^{*}{u}^{n}\|_{2}^{2}$
(3.22)
$\displaystyle+\|G_{i}^{**}u_{n,i}\|_{2}^{2}+\|G_{i}^{**}{u}^{n}\|_{2}^{2})\leq
2\tau\langle G_{i}N_{\kappa}(u_{n,j}),(I+\Delta_{N}^{2})u_{n,i}\rangle.$

The RHS of inequality (3.22) contains two parts:

$\displaystyle 2\langle
G_{i}N_{\kappa}(u_{n,j}),(I+\Delta_{N}^{2})u_{n,i}\rangle=\langle-2G_{i}u^{3}_{n,j}+2(\kappa+\varepsilon)G_{i}u_{n,j},(I+\Delta_{N}^{2})u_{n,i}\rangle.$
(3.23)

For $i=1$, the first term could be analyzed as follows:

$\displaystyle-2\langle
G_{1}(u^{n})^{3},(I+\Delta_{N}^{2})u_{n,1}\rangle=-2\langle
G_{1}^{\frac{1}{2}}(u^{n})^{3},G_{1}^{\frac{1}{2}}(I+\Delta_{N}^{2})u_{n,1}\rangle$
(3.24) $\displaystyle\leq
2\|({u}^{n})^{3}\|_{2}\cdot\|G_{1}^{\frac{1}{2}}(I+\Delta_{N}^{2})u_{n,1}\|_{2}\leq
16\|(u^{n})^{3}\|_{2}^{2}+\frac{1}{8}(\|G_{1}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}+\|G_{1}^{\frac{1}{2}}\Delta_{N}^{2}u_{n,1}\|_{2}^{2}),$

in which the inequality $(a+b)^{2}\leq 2(a^{2}+b^{2})$, summations by parts
formula, discrete Cauchy–Schwartz inequality, as well as 3.1, have been
applied in the analysis. The second term of equality (3.23) could be
decomposed into two parts:

$\displaystyle 2(\kappa-1)\langle
G_{1}u^{n},(I+\Delta_{N}^{2})u_{n,1}\rangle\leq
2(\kappa-1)[\|G_{1}^{\frac{1}{2}}{u}^{n}\|_{2}\cdot(\|G_{1}^{\frac{1}{2}}\Delta_{N}^{2}u_{n,1}\|_{2}+\|G_{1}^{\frac{1}{2}}u_{n,1}\|_{2})]$
(3.25)
$\displaystyle\leq(\kappa-1)(\|G_{1}^{\frac{1}{2}}{u}^{n}\|_{2}^{2}+\|G_{1}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}+\|G_{1}^{\frac{1}{2}}\Delta_{N}{u}^{n}\|_{2}^{2}+\|G_{1}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}^{2});$
$\displaystyle 2(1+\varepsilon)\langle
G_{1}u^{n},(I+\Delta_{N}^{2})u_{n,1}\rangle\leq
64\|{u}^{n}\|_{2}^{2}+\frac{1}{16}\|G_{1}^{\frac{1}{2}}(I+\Delta_{N}^{2})u_{n,1}\|_{2}^{2}$
(3.26) $\displaystyle\leq
64\|{u}^{n}\|_{2}^{2}+\frac{1}{8}(\|G_{1}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}+\|G_{1}^{\frac{1}{2}}\Delta_{N}^{2}u_{n,1}\|_{2}^{2}).$

Meanwhile, we observe the following inequality, as given by inequalitiy (3.7)
in 3.2:

$\displaystyle\|G_{1}^{*}u_{n,1}\|_{2}^{2}+\|G_{1}^{*}u^{n}\|_{2}^{2}+\|G_{1}^{**}u_{n,1}\|_{2}^{2}+\|G_{1}^{**}{u}^{n}\|_{2}^{2}\geq\frac{2}{3}(\|G_{1}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}+\|G_{1}^{\frac{1}{2}}u^{n}\|_{2}^{2})$
(3.27)
$\displaystyle+\frac{1}{4}(\|G_{1}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}^{2}+\|G_{1}^{\frac{1}{2}}\Delta_{N}u^{n}\|_{2}^{2}+\|G_{1}^{\frac{1}{2}}\Delta_{N}^{2}u_{n,1}\|_{2}^{2}+\|G_{1}^{\frac{1}{2}}\Delta_{N}^{2}{u}^{n}\|_{2}^{2})$
$\displaystyle+(\kappa-1)(\|G_{1}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}+\|G_{1}^{\frac{1}{2}}u^{n}\|_{2}^{2}+\|G_{1}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}^{2}+\|G_{1}^{\frac{1}{2}}\Delta_{N}{u}^{n}\|_{2}^{2}).$

Subsequently, a substitution of inequalities (3.24)-(3.27) into (3.22) leads
to

$\displaystyle\|u_{n,1}\|_{2}^{2}-\|{u}^{n}\|_{2}^{2}+\|\Delta_{N}u_{n,1}\|_{2}^{2}-\|\Delta_{N}{u}^{n}\|_{2}^{2}+\frac{\tau}{4}(\|G_{1}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}^{2}+\|G_{1}^{\frac{1}{2}}\Delta_{N}{u}^{n}\|_{2}^{2}$
$\displaystyle+\|G_{1}^{\frac{1}{2}}\Delta_{N}^{2}u^{n}\|_{2}^{2})+\frac{5}{12}\tau\|G_{1}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}+\frac{2}{3}\tau\|G_{1}^{\frac{1}{2}}u^{n}\|_{2}^{2}\leq
64\tau\|{u}^{n}\|_{2}^{2}+16\tau\|({u}^{n})^{3}\|_{2}^{2},$

or in an equivalent manner,

$\|u_{n,1}\|_{2}^{2}+\|\Delta_{N}u_{n,1}\|_{2}^{2}\leq(1+64\tau)\|{u}^{n}\|_{2}^{2}+16\tau\|({u}^{n})^{3}\|_{2}^{2}+\|\Delta_{N}u^{n}\|_{2}^{2}\leq(2+64\tau)C_{0}^{2}+16\tau\tilde{C}_{0}^{4}C_{0}^{2},$

where the following inequality has been employed:

$\|(u^{n})^{3}\|_{2}^{2}=\|u^{n}\|_{\infty}^{4}\cdot\|u^{n}\|_{2}^{2}\leq\tilde{C}_{0}^{4}C_{0}^{2}.$
(3.28)

Under an $\mathcal{O}(1)$ constraint for the time-step size

$\tau\leq\min\left\\{\frac{1}{64},\frac{1}{16}\tilde{C}_{0}^{-4}\right\\},$
(3.29)

we see that

$\|u_{n,1}\|_{2}^{2}+\|\Delta_{N}u_{n,1}\|_{2}^{2}\leq
4C_{0}^{2},\quad\text{so that}\leavevmode\nobreak\
\|u_{n,1}\|_{2}+\|\Delta_{N}u_{n,1}\|_{2}\leq 2\sqrt{2}C_{0}.$

Therefore, we obtain a rough $\ell^{2}$, $H_{h}^{2}$, and $\ell^{\infty}$
estimates of $u_{n,1}$:

$\|u_{n,1}\|_{2},\leavevmode\nobreak\ \|\Delta_{N}u_{n,1}\|_{2}\leq
2C_{0};\quad\|u_{n,1}\|_{\infty}\leq\hat{C}(\|u_{n,1}\|_{2}+\|\Delta_{N}u_{n,1}\|_{2})\leq
2\sqrt{2}\hat{C}C_{0}=:\tilde{C}_{1}.$ (3.30)

####  3.2.3 Preliminary $\ell^{2}$ and $H_{h}^{2}$ estimates of
$G_{2}^{\frac{1}{2}}u_{n,1}$

In addition, preliminary estimates of $\|G_{2}^{\frac{1}{2}}u_{n,1}\|_{2}$ and
$\|G_{2}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}$ are needed to obtain a refined
bound of the numerical solution at the next time step. Again, the reformulated
numerical system is used in this estimate. Therefore, we take $i=1$ in the
system (3.16)-(3.17).

Taking a discrete $\ell^{2}$ inner product with equality (3.16) by
$G_{2}(I+\Delta_{N}^{2})(u_{n,1}^{*}+u^{n})$ results in

$\displaystyle\|G_{2}^{\frac{1}{2}}u_{n,1}^{*}\|_{2}^{2}-\|G_{2}^{\frac{1}{2}}u^{n}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}u_{n,1}^{*}\|_{2}^{2}-\|G_{2}^{\frac{1}{2}}\Delta_{N}u^{n}\|_{2}^{2}+\tau(\|G^{*}u^{n}\|_{2}^{2}+\|G^{**}u^{n}\|_{2}^{2}$
(3.31) $\displaystyle+\langle
G_{1}L_{\kappa}u^{n},G_{2}(I+\Delta_{N}^{2})u_{n,1}^{*}\rangle)=0.$

Taking a discrete $\ell^{2}$ inner product with equality (3.17) by
$2G_{2}(I+\Delta_{N}^{2})u_{n,1}$ gives

$\displaystyle\|G_{2}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}-\|G_{2}^{\frac{1}{2}}u_{n,1}^{*}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}(u_{n,1}-u_{n,1}^{*})\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}^{2}-\|G_{2}^{\frac{1}{2}}\Delta_{N}u_{n,1}^{*}\|_{2}^{2}$
(3.32)
$\displaystyle+\|G_{2}^{\frac{1}{2}}\Delta_{N}(u_{n,1}-u_{n,1}^{*})\|_{2}^{2}=2\tau\langle
G_{1}N_{\kappa}(u^{n}),G_{2}(I+\Delta_{N}^{2})u_{n,1}\rangle.$

A combination of equalities (3.31) and (3.32) yields

$\displaystyle\|G_{2}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}-\|G_{2}^{\frac{1}{2}}u^{n}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}(u_{n,1}-u_{n,1}^{*})\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}^{2}-\|G_{2}^{\frac{1}{2}}\Delta_{N}u^{n}\|_{2}^{2}$
(3.33)
$\displaystyle+\|G_{2}^{\frac{1}{2}}\Delta_{N}(u_{n,1}-u_{n,1}^{*})\|_{2}^{2}+\tau(\|G^{*}u^{n}\|_{2}^{2}+\|G^{**}u^{n}\|_{2}^{2}+\langle
G_{1}L_{\kappa}u^{n},G_{2}(I+\Delta_{N}^{2})u_{n,i}^{*}\rangle)$
$\displaystyle=2\tau\langle
G_{1}N_{\kappa}(u^{n}),G_{2}(I+\Delta_{N}^{2})u_{n,1}\rangle.$

Meanwhile, an application of inequalities (3.11)-(3.12) in 3.3 indicates that

$\displaystyle\tau\langle
G_{1}L_{\kappa}u^{n},G_{2}(I+\Delta_{N}^{2})u_{n,1}^{*}\rangle+\|G_{2}^{\frac{1}{2}}(u_{n,1}-u_{n,1}^{*})\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}(u_{n,1}-u_{n,1}^{*})\|_{2}^{2}$
$\displaystyle\geq\tau(\|G^{*}u_{n,1}\|_{2}^{2}+\|G^{**}u_{n,1}\|_{2}^{2}).$

Going back (3.33), we arrive at

$\displaystyle\|G_{2}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}-\|G_{2}^{\frac{1}{2}}u^{n}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}^{2}-\|G_{2}^{\frac{1}{2}}\Delta_{N}u^{n}\|_{2}^{2}+\tau(\|G^{*}u^{n}\|_{2}^{2}+\|G^{*}u_{n,1}\|_{2}^{2}$
(3.34) $\displaystyle+\|G^{**}u^{n}\|_{2}^{2}+\|G^{**}u_{n,1}\|_{2}^{2})\leq
2\tau\langle G_{1}N_{\kappa}(u^{n}),G_{2}(I+\Delta_{N}^{2})u_{n,1}\rangle.$

The RHS of (3.34) contains two parts:

$\displaystyle 2\langle
G_{1}N_{\kappa}(u^{n}),G_{2}(I+\Delta_{N}^{2})u_{n,1}\rangle=-2\langle
G_{1}(u^{n})^{3},G_{2}(I+\Delta_{N}^{2})u_{n,1}\rangle$ (3.35)
$\displaystyle+2(\kappa+\varepsilon)\langle
G_{1}u^{n},G_{2}(I+\Delta_{N}^{2})u_{n,1}\rangle.$

The first term could be analyzed as follows:

$\displaystyle-2\langle
G_{1}(u^{n})^{3},G_{2}(I+\Delta_{N}^{2})u_{n,1}\rangle\leq
2\|(u^{n})^{3}\|_{2}\cdot\|G_{1,2}^{\frac{1}{2}}(I+\Delta_{N}^{2})u_{n,1}\|_{2}$
(3.36) $\displaystyle\leq
16\tilde{C}_{0}^{4}C_{0}^{2}+\frac{1}{16}\|G_{1,2}^{\frac{1}{2}}(I+\Delta_{N}^{2})u_{n,1}\|_{2}^{2}\leq
16\tilde{C}_{0}^{4}C_{0}^{2}+\frac{1}{8}(\|G_{1,2}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}^{2}u_{n,1}\|_{2}^{2}).$

The second term of (3.35) could be decomposed as two parts:

$\displaystyle 2(1+\varepsilon)\langle
G_{1}u^{n},G_{2}(I+\Delta_{N}^{2})u_{n,1}\rangle\leq
2(1+\varepsilon)\|u^{n}\|_{2}\cdot\|G_{1,2}^{\frac{1}{2}}(I+\Delta_{N}^{2})u_{n,1}\|_{2}$
(3.37) $\displaystyle\leq
64C_{0}+\frac{1}{16}\|G_{1,2}^{\frac{1}{2}}(I+\Delta_{N}^{2})u_{n,1}\|_{2}^{2}\leq
64C_{0}+\frac{1}{8}(\|G_{1,2}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}^{2}u_{n,1}\|_{2}^{2});$
$\displaystyle 2(\kappa-1)\langle
G_{1}u^{n},G_{2}(I+\Delta_{N}^{2})u_{n,1}\rangle\leq(\kappa-1)(\|G_{1,2}^{\frac{1}{2}}u^{n}\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}u^{n}\|_{2}^{2}$
(3.38) $\displaystyle+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}^{2}).$

For the two positive $\kappa$-independent terms in (3.38), we observe the
following inequality, as given by (3.8) in 3.2:

$\displaystyle\|G^{*}u^{n}\|_{2}^{2}+\|G^{*}u_{n,1}\|_{2}^{2}+\|G^{**}u^{n}\|_{2}^{2}+\|G^{**}u_{n,1}\|_{2}^{2}\geq\frac{2}{3}(\|G_{1,2}^{\frac{1}{2}}u^{n}\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}$
(3.39)
$\displaystyle+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}u^{n}\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}^{2})+\frac{1}{4}(\|G_{1,2}^{\frac{1}{2}}\Delta_{N}u^{n}\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}^{2}u^{n}\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}^{2}$
$\displaystyle+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}^{2}u_{n,1}\|_{2}^{2})+(\kappa-1)(\|G_{1,2}^{\frac{1}{2}}u^{n}\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}u^{n}\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}^{2}).$

As a result, a substitution of (3.36)-(3.39) into (3.34) yields

$\displaystyle\|G_{2}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}-\|G_{2}^{\frac{1}{2}}u^{n}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}^{2}-\|G_{2}^{\frac{1}{2}}\Delta_{N}u^{n}\|_{2}^{2}+\frac{2\tau}{3}\|G_{1,2}^{\frac{1}{2}}u^{n}\|_{2}^{2}+\frac{5\tau}{12}\|G_{1,2}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}$
$\displaystyle+\frac{\tau}{4}(\|G_{1,2}^{\frac{1}{2}}\Delta_{N}u^{n}\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}^{2}u^{n}\|_{2}^{2}+\|G_{1,2}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}^{2})\leq\tau(64C_{0}+16\tilde{C}_{0}^{4}C_{0}^{2}).$

Consequently, the following combination of $\ell^{2}$ and $H_{h}^{2}$ bounds
become available for $G_{2}^{\frac{1}{2}}u_{n,1}$:

$\|G_{2}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}^{2}\leq\|G_{2}^{\frac{1}{2}}u^{n}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}u^{n}\|_{2}^{2}+\tau(64C_{0}+16\tilde{C}_{0}^{4}C_{0}^{2}).$
(3.40)

####  3.2.4 Preliminary $\ell^{\infty}$ estimate of $u^{n+1}$

We aim to derive a bound for $\|u^{n+1}\|_{\infty}$. By taking $i=2$, the
first term of equality (3.23) can be analyzed as follows:

$2\langle G_{2}u_{n,1}^{3},(I+\Delta_{N}^{2})u^{n+1}\rangle\leq
16\|u_{n,1}^{3}\|_{2}^{2}+\frac{1}{8}(\|G_{2}^{\frac{1}{2}}u^{n+1}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}^{2}u^{n+1}\|_{2}^{2}).$
(3.41)

The second term is also decomposed as two parts:

$\displaystyle 2(\kappa-1)\langle
G_{2}u_{n,1},(I+\Delta_{N}^{2})u^{n+1}\rangle$ (3.42)
$\displaystyle\leq(\kappa-1)(\|G_{2}^{\frac{1}{2}}u_{n,1}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}{u}^{n+1}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}u_{n,1}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}u^{n+1}\|_{2}^{2})$
$\displaystyle\leq(\kappa-1)[(64C_{0}+16\tilde{C}_{0}^{4}C_{0}^{2})\tau+\|G_{2}^{\frac{1}{2}}u^{n}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}u^{n}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}u^{n+1}\|_{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}u^{n+1}\|_{2}^{2}];$
$\displaystyle 2(1+\varepsilon)\langle
G_{2}u_{n,1},(I+\Delta_{N}^{2}){u}^{n+1}\rangle\leq
64\|u_{n,1}\|_{2}^{2}+\frac{1}{8}(\|G_{2}^{\frac{1}{2}}{u}^{n+1}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}^{2}u^{n+1}\|_{2}^{2}).$

where the preliminary estimate (3.40) for the intermediate-stage solution has
been applied in the last step. Following inequality (3.7) in 3.2, we have

$\displaystyle\|G_{2}^{*}{u}^{n+1}\|_{2}^{2}+\|G_{2}^{*}{u}^{n}\|_{2}^{2}+\|G_{2}^{**}{u}^{n+1}\|_{2}^{2}+\|G_{2}^{**}{u}^{n}\|_{2}^{2}\geq\frac{1}{4}(\|G_{2}^{\frac{1}{2}}\Delta_{N}{u}^{n+1}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}{u}^{n}\|_{2}^{2}$
(3.43)
$\displaystyle+\|G_{2}^{\frac{1}{2}}\Delta_{N}^{2}{u}^{n+1}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}^{2}{u}^{n}\|_{2}^{2})+(\kappa-1)(\|G_{2}^{\frac{1}{2}}{u}^{n+1}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}{u}^{n}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}{u}^{n+1}\|_{2}^{2}$
$\displaystyle+\|G_{2}^{\frac{1}{2}}\Delta_{N}{u}^{n}\|_{2}^{2})+\frac{2}{3}(\|G_{2}^{\frac{1}{2}}u^{n+1}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}u^{n+1}\|_{2}^{2}).$

Of course, a substitution of inequalities (3.41)-(3.43) into (3.22) leads to

$\displaystyle\|u^{n+1}\|_{2}^{2}-\|u^{n}\|_{2}^{2}+\|\Delta_{N}u^{n+1}\|_{2}^{2}-\|\Delta_{N}u^{n}\|_{2}^{2}+\frac{\tau}{4}(\|G_{2}^{\frac{1}{2}}\Delta_{N}^{2}u^{n}\|_{2}^{2}+\|G_{2}^{\frac{1}{2}}\Delta_{N}{u}^{n+1}\|_{2}^{2}$
$\displaystyle+\|G_{2}^{\frac{1}{2}}\Delta_{N}{u}^{n}\|_{2}^{2})+\frac{5}{12}\tau\|G_{2}^{\frac{1}{2}}u^{n+1}\|_{2}^{2}\leq
64\tau\|u_{n,1}\|_{2}^{2}+16\tau\|u_{n,1}^{3}\|_{2}^{2}+(\kappa-1)(64C_{0}+16\tilde{C}_{0}^{4}C_{0}^{2})\tau^{2},$

or, equivalently,

$\displaystyle\|{u}^{n+1}\|_{2}^{2}+\|\Delta_{N}u^{n+1}\|_{2}^{2}\leq\|{u}^{n}\|_{2}^{2}+\|\Delta_{N}u^{n}\|_{2}^{2}+64\tau\|u_{n,1}\|_{2}^{2}+16\tau\|u_{n,1}^{3}\|_{2}^{2}$
$\displaystyle+(\kappa-1)(64C_{0}+16\tilde{C}_{0}^{4}C_{0}^{2})\tau^{2}\leq
2C_{0}^{2}+256C_{0}^{2}\tau+64\tau\tilde{C}_{1}^{4}C_{0}^{2}+(\kappa-1)(64C_{0}+16\tilde{C}_{0}^{4}C_{0}^{2})\tau^{2},$

in which the following inequality has been used:

$\|u_{n,1}^{3}\|_{2}^{2}=\|u_{n,1}\|_{\infty}^{4}\|u_{n,1}\|_{2}^{2}\leq
4\tilde{C}_{1}^{4}C_{0}^{2}.$

Under an $\mathcal{O}(1)$ constraint for the time-step size (more constraint
than (3.29))

$\tau\leq\min\left\\{\frac{1}{256},\frac{1}{64}\tilde{C}_{1}^{-4},(64\kappa)^{-\frac{1}{2}},\frac{1}{4}\tilde{C}_{0}^{-2}\kappa^{-\frac{1}{2}}\right\\},$
(3.44)

we get

$\|{u}^{n+1}\|_{2}^{2}+\|\Delta_{N}u^{n+1}\|_{2}^{2}\leq
6C_{0}^{2},\quad\text{so that}\leavevmode\nobreak\
\|{u}^{n+1}\|_{2}+\|\Delta_{N}u^{n+1}\|_{2}\leq 2\sqrt{3}C_{0}.$ (3.45)

Also note that $C_{0}$ is $\kappa$-independent and time-independent. Again, an
application of 3.2 implies the following $\|\cdot\|_{\infty}$ bound at the
next time step

$\|{u}^{n+1}\|_{\infty}\leq\hat{C}\left(\|{u}^{n+1}\|_{2}+\|\Delta_{N}{u}^{n+1}\|_{2}\right)\leq
2\sqrt{3}\hat{C}C_{0}=:\tilde{C}_{2}.$ (3.46)

Notice that $\tilde{C}_{2}$ is also a $\kappa$-independent and global-in-time
constant.

####  3.2.5 Justification of the stabilization parameter $\kappa$ and a-priori
assumption (3.13)

On the other hand, by making comparison between the $\ell^{\infty}$ bounds for
$u_{n,i}\leavevmode\nobreak\ (i=0,1,2)$, given by (3.15), (3.30), and (3.46),
respectively, it is clear that
$\tilde{C}_{2}\geq\tilde{C}_{1}\geq\tilde{C}_{0}$ for the ERK(2,2) scheme
(2.12). In turn, we could take

$\kappa=\max\left\\{\frac{|3\tilde{C}_{2}^{2}-\varepsilon|}{2},1\right\\},$
(3.47)

in which $\tilde{C}_{2}\geq\beta=\max_{i=0,1,2}\|u_{n,i}\|_{\infty}$. Notice
that $\kappa$ is an $\mathcal{O}(1)$ constant, and contains no singular
dependence on any physical parameter. With this choice of $\kappa$, a fixed
constant, we could take the time-step size $\tau$ satisfying (3.44) for
ERK(2,2), such that original energy stability becomes available at the next
time step by 1.1:

$E_{N}({u}^{n+1})\leq E_{N}({u}^{n})\leq E_{N}(u^{0})=C_{e}.$ (3.48)

This in turn recovers the a-priori assumption (3.13) at the next time step, so
that an induction argument can be effectively applied. Therefore, we have
proved the main theorem 1.2.

######  Remark 3.1.

Obviously, the above $\ell^{2}$, $H_{h}^{2}$, and $\ell^{\infty}$ estimates of
$u^{n+1}$, namely (3.45) and (3.46), turn out to be too rough, since we did
not make use of the variational energy structure in the analysis. In fact, to
obtain an energy dissipation at the theoretical level, an $\ell^{\infty}$
bound of the numerical solution at the time step $t^{n+1}$ has to be derived,
due to the nonlinear term involved. On the other hand, with such a rough bound
at hand, we are able to justify the artificial parameter value in (3.47), so
that energy stability becomes theoretically available at the next time step.
With a theoretical justification of the energy stability analysis, we are able
to obtain much sharper $\ell^{2}$, $H_{h}^{2}$, and $\ell^{\infty}$ bounds for
the numerical solution $u^{n+1}$.

In more details, with the energy stability result (3.48), we apply a similar
analysis in (3.14) and obtain

$\|u^{n+1}\|_{2},\leavevmode\nobreak\ \|\Delta_{N}u^{n+1}\|_{2}\leq
C_{0}:=2\sqrt{C_{e}+|\Omega|},$

which is a global-in-time constant. In turn, a much sharper maximum-norm bound
for $u^{n+1}$ also becomes available, with the help of 3.2:

$\|u^{n+1}\|_{\infty}\leq\hat{C}(\|u^{n+1}\|_{2}+\|\Delta_{N}u^{n+1}\|_{2})\leq
2\hat{C}C_{0}=:\tilde{C}_{3}.$

In other words, the $H_{h}^{2}$ bound $C_{0}$ and the $\ell^{\infty}$ bound
$\tilde{C}_{3}$ turns out to be global-in-time constants.

##  4 Optimal rate convergence analysis

We denote by $u_{e}$ the exact solution to equation (1.6), and assume it
satisfies the following regularity:

$u_{e}\in\mathcal{R}:=H^{3}(0,T;C^{0})\cap H^{2}(0,T;H^{m_{0}})\cap
L^{\infty}(0,T;H^{m_{0}+4}).$

A rigorous error estimate for the ERK(2,2) scheme (2.12) will be derived under
this regularity. To this end, the following lemma is needed.

######  Lemma 4.1.

For any ${u},{v}\in\mathcal{M}_{N}$ satisfying
$\|{u}\|_{\infty},\leavevmode\nobreak\ \|{v}\|_{\infty}\leq\beta$ that is
introduced in 1.1, we have

$\|N_{\kappa}({u})-N_{\kappa}({v})\|_{2}\leq 3\kappa\|{u}-{v}\|_{2}.$

###### Proof.

Since $N_{\kappa}(u)-N_{\kappa}(v)=(u-v)(\varepsilon+\kappa-u^{2}-v^{2}-uv)$,
using $\|u\|_{\infty}\leq\beta$, $\|v\|_{\infty}\leq\beta$, and
$\kappa\geq\max\limits_{|\xi|\leq\beta}\frac{|3\xi^{2}-\varepsilon|}{2}$, we
obtain the result. ∎

Meanwhile, we denote $U^{n}$ as the interpolation values of the projection
solution $U_{N}$ at the grid points at time instant
$t_{n}:\leavevmode\nobreak\ U_{p,q}^{n}:=U_{N}(x_{i},y_{j},t_{n})$. The
initial data is given by

$u_{p,q}^{0}=U_{p,q}^{0}:=U_{N}(x_{p},y_{q},t=0).$

The error grid function is defined as

$e^{n}:=U^{n}-u^{n},\quad 0\leq n\leq N_{t}-1.$

For the proposed ERK(2,2) scheme (2.12), the convergence result is stated
below.

######  Theorem 4.1.

Given an initial data with sufficient regularity, suppose the unique solution
for the SH equation (1.6) is of regularity class $\mathcal{R}$. Provided that
$\tau$ and $h$ are sufficiently small, then, for the ERK(2,2) scheme (2.12),
the following $\ell^{\infty}(0,T;\ell^{2})$ convergence estimate is valid for
any $\kappa$ satisfying (3.47):

$\|u_{e}(t^{n})-u^{n}\|_{2}\leq C(\tau^{2}+h^{m_{0}}),\quad\forall n\leq
N_{t},$

where $C>0$ is dependent of $\kappa$ and $\Omega$, but independent of $\tau$
and $h$.

###### Proof.

For the exact solution $u_{e}$ and its interpolation $U$, a careful
consistency analysis implies that

$\displaystyle U_{n,1}$ $\displaystyle=\varphi_{0}(c_{1}\tau
L_{\kappa})U^{n}+c_{1}\tau\varphi_{1}(c_{1}\tau L_{\kappa})N_{\kappa}(U^{n}),$
(4.1) $\displaystyle U^{n+1}$ $\displaystyle=\varphi_{0}(\tau
L_{\kappa})U^{n}+\tau\varphi_{1}(\tau
L_{\kappa})N_{\kappa}(U_{n,1})+\tau\zeta^{n},$

with $\|\zeta^{n}\|_{2}\leq\bar{C}(\tau^{2}+h^{m_{0}})$. Notice that the
profile $U_{n,1}$ is constructed, based on the projection solution $U^{n}$. In
turn, subtracting the numerical solution (2.12) from the consistency estimate
(4.1) yields

$\displaystyle e_{n,1}$ $\displaystyle=\varphi_{0}(c_{1}\tau
L_{\kappa})e^{n}+c_{1}\tau\varphi_{1}(c_{1}\tau
L_{\kappa})\tilde{N}_{\kappa}(U^{n},{u}^{n}),$ (4.2) $\displaystyle e^{n+1}$
$\displaystyle=\varphi_{0}(\tau L_{\kappa})e^{n}+\tau\varphi_{1}(\tau
L_{\kappa})\tilde{N}_{\kappa}(U_{n,1},u_{n,1})+\tau\zeta^{n},$

with $\tilde{N}_{\kappa}(a,b)=N_{\kappa}(a)-N_{\kappa}(b)$. To carry out the
error analysis in a more convenient way, we denote
$e_{n,i}^{*}=\varphi_{0}\left(c_{i}\tau L_{\kappa}\right)e^{n}$ ($i=1,2$), so
that the evolutionary equation (4.2) could be rewritten as the following two-
substage system:

$\displaystyle\frac{e_{n,i}^{*}-e^{n}}{\tau}$
$\displaystyle=-c_{i}L_{\kappa}\varphi_{1}(c_{i}\tau L_{\kappa})e^{n},$ (4.3)
$\displaystyle\frac{e_{n,i}-e_{n,i}^{*}}{\tau}$
$\displaystyle=c_{i}\varphi_{1}(c_{i}\tau
L_{\kappa})\tilde{N}_{\kappa}(U_{n,j},u_{n,j})+(i-1)\zeta^{n}.$ (4.4)

Taking a discrete $\ell^{2}$ inner product with (4.3) by $e_{n,i}^{*}+e^{n}$
gives

$\|e_{n,i}^{*}\|_{2}^{2}-\|e^{n}\|_{2}^{2}+\tau(\|G_{i}^{*}e^{n}\|_{2}^{2}+\langle
G_{i}L_{\kappa}e^{n},e_{n,i}^{*}\rangle)=0.$ (4.5)

Taking a discrete $\ell^{2}$ inner product with (4.4) by $2e_{n,i}$ yields

$\|e_{n,i}\|_{2}^{2}-\|e_{n,i}^{*}\|_{2}^{2}+\|e_{n,i}-e_{n,i}^{*}\|_{2}^{2}=2\tau\langle
G_{i}\tilde{N}_{\kappa}(U_{n,j},u_{n,j}),e_{n,i}\rangle+2(i-1)\tau\langle\zeta^{n},e_{n,i}\rangle.$
(4.6)

As a result, a combination of equalities (4.5) and (4.6) leads to

$\displaystyle\|e_{n,i}\|_{2}^{2}-\|e^{n}\|_{2}^{2}+\|e_{n,i}-e_{n,i}^{*}\|_{2}^{2}+\tau(\|G_{i}^{*}e^{n}\|_{2}^{2}+\langle
G_{i}L_{\kappa}e^{n},e_{n,i}^{*}\rangle)$ $\displaystyle=2\tau\langle
G_{i}\tilde{N}_{\kappa}(U_{n,j},u_{n,j}),e_{n,i}\rangle+2(i-1)\tau\langle\zeta^{n},e_{n,i}\rangle.$

Meanwhile, an application of (3.7) in 3.2 and (3.9) in 3.3 results in

$\displaystyle\|G_{i}^{*}e_{n,i}\|_{2}^{2}\geq\frac{1}{4}\|G_{i}^{\frac{1}{2}}\Delta_{N}e_{n,i}\|_{2}^{2}+(\kappa-1)\|G_{i}^{\frac{1}{2}}e_{n,i}\|_{2}^{2}+\frac{2}{3}\|G_{i}^{\frac{1}{2}}e_{n,i}\|_{2}^{2};$
$\displaystyle\tau\langle
G_{i}L_{\kappa}e^{n},e_{n,i}^{*}\rangle+\|e_{n,i}-e_{n,i}^{*}\|_{2}^{2}\geq\tau\|G_{i}^{*}e_{n,i}\|_{2}^{2}.$

Then we obtain

$\displaystyle\|e_{n,i}\|_{2}^{2}-\|e^{n}\|_{2}^{2}+\tau(\|G_{i}^{*}e^{n}\|_{2}^{2}+\frac{1}{4}\|G_{i}^{\frac{1}{2}}\Delta_{N}e_{n,i}\|_{2}^{2}+(\kappa-1)\|G_{i}^{\frac{1}{2}}e_{n,i}\|_{2}^{2}+\frac{2}{3}\|G_{i}^{\frac{1}{2}}e_{n,i}\|_{2}^{2})$
(4.7) $\displaystyle\leq 2\tau\langle
G_{i}\tilde{N}_{\kappa}(U_{n,j},u_{n,j}),e_{n,i}\rangle+2(i-1)\tau\langle\zeta^{n},e_{n,i}\rangle.$

Moreover, an application of 4.1 implies that

$\|\tilde{N}_{\kappa}(U_{n,j},u_{n,j})\|_{2}\leq 3\kappa\|e_{n,j}\|_{2},$

which would be used for the derivation of the following nonlinear inner
product term estimate ($i=1$):

$2\langle G_{1}\tilde{N}_{\kappa}(U^{n},{u}^{n}),e_{n,1}\rangle\leq
36\kappa^{2}\|e^{n}\|_{2}^{2}+\frac{1}{4}\|G_{1}^{\frac{1}{2}}e_{n,1}\|_{2}^{2}.$
(4.8)

Its substitution into inequality (4.7) yields

$\|e_{n,1}\|_{2}^{2}-\|e^{n}\|_{2}^{2}\leq 36\kappa^{2}\tau\|e^{n}\|_{2}^{2}.$

Provided that $\tau\leq(6\kappa)^{-2}$, a preliminary $\ell^{2}$ error
estimate for the intermediate-stage error solution $u_{n,1}$ is obtained

$\|e_{n,1}\|_{2}\leq\left[1+36\kappa^{2}\tau\right]^{\frac{1}{2}}\|e^{n}\|_{2}\leq\sqrt{2}\|e^{n}\|_{2}.$

For $i=2$, the analysis for the nonlinear error inner product term of
inequality (4.8) could be similarly established:

$2\langle G_{2}\tilde{N}_{\kappa}(U_{n,1},u_{n,1}),e^{n+1}\rangle\leq
72\kappa^{2}\|e^{n}\|_{2}^{2}+\frac{1}{4}\|G_{2}^{\frac{1}{2}}e^{n+1}\|_{2}^{2},$
(4.9)

in which the estimate $\|e_{n,1}\|_{2}\leq\sqrt{2}\|e^{n}\|_{2}$ has been
applied in the last step. A bound for the truncation error inner product term
is more straightforward:

$2\langle\zeta^{n},e^{n+1}\rangle\leq\|\zeta^{n}\|_{2}^{2}+\|e^{n+1}\|_{2}^{2}.$
(4.10)

Afterwards, a substitution of inequalities (4.9) and (4.10) into (4.7) leads
to

$\|e^{n+1}\|_{2}^{2}-\|e^{n}\|_{2}^{2}\leq
72\kappa^{2}\tau\|e^{n}\|_{2}^{2}+\tau\|e^{n+1}\|_{2}^{2}+\tau\|\zeta^{n}\|_{2}^{2}.$

In turn, an application of the Grönwall’s inequality results in the desired
convergence estimate:

$\|e^{n+1}\|_{2}\leq C(\tau^{2}+h^{m_{0}}),$

due to the fact that $\|\zeta^{n}\|_{2}\leq\bar{C}(\tau^{2}+h^{m_{0}})$. This
validates the convergence estimate. ∎

##  5 Numerical experiments: Comparison with other energy-stable methods

In this section, we present a few 2-D numerical results for the SH equation
(1.6)), to demonstrate the efficiency, accuracy, and long-time performance of
the ERK(2,2) scheme (2.12). To preserve energy stability, the condition
$\kappa=\max\left\\{\frac{|3\tilde{C}_{2}^{2}-\varepsilon|}{2},1\right\\}$ is
theoretically required in 1.2. Although there is no maximum principle for the
SH equation (1.6), in practice we observe that the numerical solution is
always bounded in $[-1,1]$ during the whole simulation process. Therefore, it
suffices to set $\kappa=2$ in the computation.

###  5.1 Convergence test

The computational domain is taken as $\Omega=[0,32]^{2}$, and the following
smooth initial data is enforced:

$u_{0}(x,y)=0.01\times[\cos(\pi x)+\cos(\pi y)+\cos(0.25\pi x)+\cos(0.25\pi
y)],$

on the uniform mesh $N=N_{x}=N_{y}$ with $N=256$. The final time is set as
$T=5$. To obtain the numerical errors, we set the numerical solution obtained
by the ERK(2,2) scheme with $\tau=0.1\times 2^{-9}$ as the reference one.
Afterwards, we perform the numerical simulations of the first-order ETD1,
IMEX1, and second-order ERK(2,2), ETDRK2, IMEX-RK(2,2) schemes using time-step
sizes $\tau=2^{-k}(k=1,\ldots,9$), with two different parameters,
$\varepsilon=0.25$ and $0.025$. The results are displayed in Figure 1, and
through observation we observe that the exponential-type schemes consistently
outperform the non-exponential-type ones, in terms of computational accuracy
and efficiency. Further, ERK(2,2) does the best.

Figure 1: Convergence of the first- and second-order Fourier pseudo-spectral
schemes in time with fixed $\tau$ (left) and $\tau^{2}$ (right) for the 2-D SH
equation. Top: $\varepsilon=0.25$. Bottom: $\varepsilon=0.025$. It is seen
that the numerical error magnitude and computational cost of the ERK(2,2)
scheme are smaller than those of both the ETDRK2 and IMEX-RK(2,2) schemes,
although they share the same convergence order.

###  5.2 Energy stable test

We simulate the pattern formation and evolution of the SH equation (1.6),
which arises in the Rayleigh–Bénard convection. We conduct the simulation on a
rectangular domain $\Omega=[0,100]^{2}$ from $T=0$ to $100$, subject to the
following initial condition:

$u_{0}(x,y)=0.1+0.02\times\cos\left(\frac{\pi
x}{100}\right)\sin\left(\frac{\pi y}{100}\right)+0.05\times\sin\left(\frac{\pi
x}{20}\right)\cos\left(\frac{\pi y}{20}\right).$

Setting the spatial mesh $256\times 256$, we compare the original energy
evolution by different time step sizes and energy-stable methods, which are
depicted in Figure 2. In actual implementation, we discover that ERK(2,2)
allows a rather mild time-step restriction, which will be useful in future
work on large-scale scientific computing.

(a)

(b)

Figure 2: The evolution of the original energy using ERK(2,2) with different
time-step sizes (left) and a fixed size $\tau=0.1$ with various energy-stable
methods (right) is shown. It is observed that, while only the red dashed line
($\tau=0.1$) approximately matches the reference line, the other two lines
also exhibit the same energy-decreasing trend and eventually converge to the
same steady state. From the right subplot, it can be seen that all methods
maintain discrete energy stability over extended periods, but ERK(2,2) reaches
the steady state more rapidly

###  5.3 Polycrystal growth in a supercooled liquid

In the existing studies [20, 22, 24], the polycrystal growth in a supercooled
liquid was considered as an important 2-D benchmark test. Here, we look at the
growth of three crystal nucleuses with the following initial data:

$u_{0}(x,y)=0.287+\alpha\times\mathrm{rand}(x,y),$

in which $\mathrm{rand}(x,y)$ is the random number uniformly distributed
between $-1$ and $1$, and $\alpha$ takes the values of
$0.1,\leavevmode\nobreak\ 0.2$, and $0.4$ for three crystal nucleuses locating
at $(375,125)$, $(375,375)$, and $(125,250)$, respectively. The length of each
nucleus is $10$. We also set the computational domain, spatial resolution, and
time-step size as $(0,500)^{2}$, $512\times 512$, and $0.5$, respectively. In
this test, the parameter $\varepsilon$ is chosen to be $0.25$. Figure 3
displays snapshots of the crystal micro-structure at several time instants.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Evolution of 2-D polycrystal growth in a supercooled liquid at
$T=16,\leavevmode\nobreak\ 40,\leavevmode\nobreak\ 72,\leavevmode\nobreak\
96,\leavevmode\nobreak\ 120$, and $160$ computed by ERK(2,2). It can be seen
that three different nuclei grains grow and eventually become sufficiently
large to form grain boundaries.

##  6 Concluding Remarks

Taking the analysis of the SH equation as a demonstration, we have proposed a
new strategy for proving discrete energy stability without any preconditions.
Further, to numerically overcome the difficulty caused by strong stiffness, we
have devised a stabilization exponential Runge–Kutta (ERK) scheme, which is
equipped with an appropriate linear stabilization term, preserving the
dissipation property of the original energy. We determined the $\ell^{\infty}$
bounds of the numerical solution at all ERK stages so that global-in-time
energy stability becomes available. Such an $\ell^{\infty}$ analysis was
accomplished by the $\ell^{2}$ and $H_{h}^{2}$ estimates of the numerical
solution at each ERK stage, with the help of summation-by-parts formulas in
the Fourier pseudo-spectral space, the discrete Sobolev inequality and
elliptic regularity, as well as the eigenvalue estimates in the Fourier space.
The global-in-time stability analysis for the original energy is the first
such result for SH-type equations. Furthermore, we have provided an optimal
rate convergence analysis for the proposed scheme, under a sufficient
regularity assumption for the exact solution. A few numerical results have
also been presented in this paper. In the convergence test, it was discovered
that ERK(2,2) outperforms ETDRK2 and IMEX-RK(2,2), in terms of efficiency and
accuracy. Besides, the long-time numerical simulation of pattern formation and
further evolution have also validated the robustness of the proposed scheme.

## Acknowledgements

This research was supported by the National Key R&D Program of China [Grant
No. 2020YFA0709803], the National Natural Science Foundation of China [Grant
Nos. 12071481 and 12271523], the Defense Science Foundation of China [Grant
No. 2021-JCJQ-JJ-0538], and the Science and Technology Innovation Program of
Hunan Province [Grant Nos. 2021RC3082 and 2022RC1192].

## References

  * [1]
  * [2] W. Chen, S. Conde, C. Wang, X. Wang, and S. M. Wise, A linear energy stable scheme for a thin film model without slope selection, Journal of Scientific Computing, 52 (2012), pp. 546–562.
  * [3] W. Chen, C. Wang, X. Wang, and S. M. Wise, A linear iteration algorithm for a second-order energy stable scheme for a thin film model without slope selection, Journal of Scientific Computing, 59 (2014), pp. 574–601.
  * [4] M. Cheng and J. A. Warren, An efficient algorithm for solving the phase field crystal model, Journal of Computational Physics, 227 (2008), pp. 6241–6248.
  * [5] M. C. Cross and P. C. Hohenberg, Pattern formation outside of equilibrium, Reviews of modern physics, 65 (1993), p. 851.
  * [6] M. Dehghan and V. Mohammadi, The numerical simulation of the phase field crystal (PFC) and modified phase field crystal (MPFC) models via global and local meshless methods, Computer Methods in Applied Mechanics and Engineering, 298 (2016), pp. 453–484.
  * [7] S. Gottlieb, F. Tone, C. Wang, X. Wang, and D. Wirosoetisno, Long time stability of a classical efficient scheme for two-dimensional Navier–Stokes equations, SIAM Journal on Numerical Analysis, 50 (2012), pp. 126–150.
  * [8] S. Gottlieb and C. Wang, Stability and convergence analysis of fully discrete Fourier collocation spectral method for 3-D viscous Burgers’ equation, Journal of Scientific Computing, 53 (2012), pp. 102–128.
  * [9] M. Hochbruck and A. Ostermann, Explicit exponential Runge–Kutta methods for semilinear parabolic problems, SIAM Journal on Numerical Analysis, 43 (2005), pp. 1069–1090.
  * [10] P. Hohenberg and J. Swift, Effects of additive noise at the onset of Rayleigh-Bénard convection, Physical Review A, 46 (1992), p. 4773.
  * [11] A. Hutt and F. M. Atay, Analysis of nonlocal neural fields for both general and gamma-distributed connectivities, Physica D: Nonlinear Phenomena, 203 (2005), pp. 30–54.
  * [12] A. Hutt, A. Longtin, and L. Schimansky-Geier, Additive noise-induced Turing transitions in spatial systems with application to neural fields and the SH equation, Physica D: Nonlinear Phenomena, 237 (2008), pp. 755–773.
  * [13] H. G. Lee, An energy stable method for the SH equation with quadratic–cubic non-linearity, Computer Methods in Applied Mechanics and Engineering, 343 (2019), pp. 40–51.
  * [14] H. G. Lee and J. Kim, A simple and efficient finite difference method for the phase-field crystal equation on curved surfaces, Computer Methods in Applied Mechanics and Engineering, 307 (2016), pp. 32–43.
  * [15] X. Li and Z. Qiao, A second-order, linear, $L^{\infty}$-convergent, and energy stable scheme for the phase field crystal equation, SIAM Journal on Scientific Computing, 46 (2024), pp. A429–A451.
  * [16] Z. Liu and C. Chen, On efficient semi-implicit auxiliary variable methods for the six-order SH model, Journal of Computational and Applied Mathematics, 419 (2023), p. 114730.
  * [17] R. Rosa, J. Pontès, C. Christov, F. M. Ramos, C. R. Neto, E. L. Rempel, and D. Walgraef, Gradient pattern analysis of SH dynamics: phase disorder characterization, Physica A: Statistical Mechanics and its Applications, 283 (2000), pp. 156–159.
  * [18] J. Sun, H. Zhang, X. Qian, and S. Song, A family of Structure-preserving Exponential Time Differencing Runge–Kutta Schemes for the Viscous Cahn–Hilliard Equation, Journal of Computational Physics, (2023), p. 112414.
  * [19] Z. Fu and J. Yang, Energy-decreasing exponential time differencing Runge–Kutta methods for phase-field models, Journal of Computational Physics, (2022), p. 110943.
  * [20] J. Swift and P. C. Hohenberg, Hydrodynamic fluctuations at the convective instability, Physical Review A, 15 (1977), p. 319.
  * [21] S. M. Wise, C. Wang, and J. S. Lowengrub, An energy-stable and convergent finite-difference scheme for the phase field crystal equation, SIAM Journal on Numerical Analysis, 47 (2009), pp. 2269–2288.
  * [22] J. Yang and J. Kim, Linear and energy stable schemes for the SH equation with quadratic-cubic non-linearity based on a modified scalar auxiliary variable approach, Journal of Engineering Mathematics, 128 (2021), p. 21.
  * [23] J. Yang, Z, Tan, and J. Kim, High-order time-accurate, efficient, and structure-preserving numerical methods for the conservative SH model, Computers & Mathematics with Applications, 102 (2021), pp. 160–174.
  * [24] J. Su, W. Fang, Q. Yu, and Y. Li, Numerical simulation of SH equation by the fourth-order compact scheme, Computational and Applied Mathematics, 38 (2019), p. 54.
  * [25] X. Feng, T. Tang, and J. Yang, Stabilized Crank-Nicolson/Adams-Bashforth schemes for phase field models, East Asain Journal on Applied Mathematics, 3 (2013), pp. 59–80.
  * [26] N. Condette, C. Melcher, and E. Suli, Spectral approximation of pattern-forming nonlinear evolution equations with double-well potentials of quadratic growth, Mathematics of Computation, 80 (2011), pp. 205–223.
  * [27] S. Pei, Y. Hou, and B. You, A linearly second-order energy stable scheme for the phase field crystal model, Applied Numerical Mathematics, 140 (2019), pp. 134–164.
  * [28] J. Shen and J. Yang, Numerical approximations of Allen–Cahn and Cahn–Hilliard equations, Discrete and Continuous Dynamical Systems, 28 (2010), pp. 1669–1691.
  * [29] M. Cui, Y. Niu, and Z. Xu, A second order exponential time differencing multi-step energy stable scheme for SH equation with quadratic-cubic nonlinear term, Journal of Scientific Computing, 99 (2024), p. 26.
  * [30] K.R. Elder and M. Grant, Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals, Physical Review E, 70 (2004), p. 051605.
  * [31] X. Li, L. Ju, and X. Meng, Convergence Analysis of Exponential Time Differencing Schemes for the Cahn–Hilliard equation, Communications in Computational Physics, 26 (2019), p. 5.
  * [32] L. Zhu, L. Ju, and W. Zhao, Fast high-order compact exponential time differencing Runge–Kutta methods for second-order semilinear parabolic equations, Journal of Scientific Computing, 67 (2016), pp. 1043–1065.
  * [33] J. Li, L. Ju, Y. Cai, and X. Feng, Unconditionally maximum bound principle preserving linear schemes for the conservative Allen–Cahn equation with nonlocal constraint, Journal of Scientific Computing, 87 (2021), pp. 1–32.
  * [34] D. Li, and Z. Qiao, On second order semi-implicit Fourier spectral methods for 2D Cahn–Hilliard equation, Journal of Scientific Computing, 70 (2017), pp. 301–341.
  * [35] D. Li, Z. Qiao, and T. Tang, Characterizing the stabilization size for semi-implicit Fourier-spectral method to phase field equations, SIAM Journal on Numerical Analysis, 54 (2016), pp. 1653–1681.
  * [36] M. Caliari, F. Cassini, L. Einkemmer, and A. Ostermann, Accelerating Exponential Integrators to Efficiently Solve Semilinear Advection-Diffusion-Reaction Equations, SIAM Journal on Scientific Computing, 46 (2024), pp. A906–A928.
  * [37] U.M. Ascher, S.J. Ruuth, and R.J. Spiteri, Implicit-explicit Runge–Kutta methods for time-dependent partial differential equations, Applied Numerical Mathematics, 25 (1997), pp. 151–167.
  * [38] Z. Fu, T Tang, and J. Yang, Energy diminishing implicit-explicit Runge–Kutta methods for gradient flows, Mathematics of Computation, (2024).
  * [39] S. Maset and M. Zennaro, Unconditionally stability of explicit exponential Runge–Kutta methods for semi-linear ordinary differential equations, Mathematics of Computation, 78 (2009), pp. 957–967.
  * [40] W. Chen, W. Li, Z. Luo, C. Wang, and X. Wang, A stabilized second order exponential time differencing multistep method for thin film growth model without slope selection, ESAIM: Mathematical Modeling and Numerical Analysis, 54 (2020), pp. 727–750.
  * [41] W. Chen, W. Li, C. Wang, S. Wang, and X. Wang, Energy stable higher-order linear ETD multi-step methods for gradient flows: application to thin film epitaxy, Research in the Mathematical Sciences, 7 (2020), pp. 1–27.

##  Appendix A Proof of 3.2

It could be easily verified that

$\displaystyle\left(-1+\lambda_{\ell,m}\right)^{2}+\kappa=2-2\lambda_{\ell,m}+\lambda_{\ell,m}^{2}+(\kappa-1)=\frac{2}{3}+\frac{4}{3}-2\lambda_{\ell,m}+\frac{3}{4}\lambda_{\ell,m}^{2}+\frac{1}{4}\lambda_{\ell,m}^{2}+(\kappa-1)$
$\displaystyle=\frac{2}{3}+(\frac{2}{\sqrt{3}}-\frac{\sqrt{3}}{2}\lambda_{\ell,m})^{2}+\frac{1}{4}\lambda_{\ell,m}^{2}+(\kappa-1),$

so that

$\Lambda_{\ell,m}=\left(-1+\lambda_{\ell,m}\right)^{2}+\kappa\geq\frac{1}{4}\lambda_{\ell,m}^{2}+\frac{2}{3}+(\kappa-1).$

This in turn leads to the following inequality:

$\|G_{i}^{*}{f}\|_{2}^{2}\geq
L^{2}\sum_{\ell,m=-K}^{K}\frac{1-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}}}{\tau\Lambda_{\ell,m}}[\frac{1}{4}\lambda_{\ell,m}^{2}+(\kappa-1)+\frac{2}{3}]|\hat{f}_{\ell,m}|^{2}.$
(A.1)

On the other hand, an application of Parseval’s equality to the discrete
Fourier expansion of $G_{i}^{\frac{1}{2}}\Delta_{N}{f}$ and
$G_{i}^{\frac{1}{2}}{f}$ reveals that

$\displaystyle\|G_{i}^{\frac{1}{2}}\Delta_{N}{f}\|_{2}^{2}$
$\displaystyle=L^{2}\sum_{\ell,m=-K}^{K}\frac{1-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}}}{\tau\Lambda_{\ell,m}}\cdot\lambda_{\ell,m}^{2}\cdot|\hat{f}_{\ell,m}|^{2},$
(A.2) $\displaystyle\|G_{i}^{\frac{1}{2}}{f}\|_{2}^{2}$
$\displaystyle=L^{2}\sum_{\ell,m=-K}^{K}\frac{1-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}}}{\tau\Lambda_{\ell,m}}\cdot|\hat{f}_{\ell,m}|^{2}.$
(A.3)

Making a comparison between inequality (A.1) and equations (A.2)-(A.3), we
conclude that inequality (3.7) has been established. Inequality (3.8) is a
direct application of the above derivation process, and the details are
skipped for the sake of brevity. This finishes the proof of 3.2.

##  Appendix B Proof of 3.3

The following expansion is assumed for the grid function $g$:

$g_{p,q}=\sum_{\ell,m=-K}^{K}\hat{g}_{\ell,m}\mathrm{e}^{2\pi\mathrm{i}(\ell
x_{p}+my_{q})/L}.$

In turn, the discrete Fourier expansion of $g-\varphi_{0}(c_{1}\tau
L_{\kappa})f$ becomes

$(g-\varphi_{0}(c_{i}\tau
L_{\kappa})f)_{p,q}=\sum_{\ell,m=-K}^{K}(\hat{g}_{\ell,m}-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}}\hat{f}_{\ell,m})\mathrm{e}^{2\pi\mathrm{i}(\ell
x_{p}+my_{q})/L},$ (B.1)

so that its discrete $\ell^{2}$ norm turns out to be

$\|g-\varphi_{0}(c_{i}\tau
L_{\kappa})f\|_{2}^{2}=L^{2}\sum_{\ell,m=-K}^{K}|\hat{g}_{\ell,m}-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}}\hat{f}_{\ell,m}|^{2}.$
(B.2)

Subsequently, a combination of equations (B.1) and (B.2) produces

$\displaystyle\tau\langle G_{i}L_{\kappa}f,\varphi_{0}(c_{i}\tau
L_{\kappa})f\rangle+\|g-\varphi_{0}(c_{i}\tau L_{\kappa})f\|_{2}^{2}$
$\displaystyle=L^{2}\sum_{\ell,m=-K}^{K}\left((1-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}})\cdot\mathrm{e}^{c_{i}\tau\Lambda_{\ell,m}}\cdot|\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}}\hat{f}_{\ell,m}|^{2}+|\hat{g}_{\ell,m}-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}}\hat{f}_{\ell,m}|^{2}\right)$
$\displaystyle=L^{2}\sum_{\ell,m=-K}^{K}(1-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}})\left({\mathrm{e}^{c_{i}\tau\Lambda_{\ell,m}}|\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}}\hat{f}_{\ell,m}|^{2}}\right.$
$\displaystyle\left.{+(1-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}})^{-1}|\hat{g}_{\ell,m}-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}}\hat{f}_{\ell,m}|^{2}}\right).$

On the other hand, for each fixed mode frequency $(\ell,m)$, the following
lower bound is clearly observed:

$\displaystyle\mathrm{e}^{c_{i}\tau\Lambda_{\ell,m}}a^{2}+(1-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}})^{-1}b^{2}=a^{2}+b^{2}+(\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}}-1)a^{2}+[(1-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}})^{-1}-1]b^{2}$
$\displaystyle\geq a^{2}+b^{2}+2ab=(a+b)^{2},$

for any $a,b\geq 0$, in which the Cauchy inequality has been applied in the
second step. Then we get

$\displaystyle\mathrm{e}^{c_{i}\tau\Lambda_{\ell,m}}|\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}}\hat{f}_{\ell,m}|^{2}+(1-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}})^{-1}|\hat{g}_{\ell,m}-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}}\hat{f}_{\ell,m}|^{2}$
$\displaystyle\geq(|\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}}\hat{f}_{\ell,m}|+|\hat{g}_{\ell,m}-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}}\hat{f}_{\ell,m}|)^{2}\geq|\hat{g}_{\ell,m}|^{2},$

so that

$\displaystyle\tau\|G_{i}L_{\kappa}f,\varphi_{0}(c_{i}\tau
L_{\kappa})f\|_{2}+\|g-\varphi_{0}(c_{i}\tau L_{\kappa})f\|_{2}^{2}\geq
L^{2}\sum_{\ell,m=-K}^{K}(1-\mathrm{e}^{-c_{i}\tau\Lambda_{\ell,m}})|\hat{g}_{\ell,m}|^{2}=\tau\|G_{i}^{*}g\|_{2}^{2}.$

Therefore, the proof of inequality (3.9) has been finished. The proofs of
inequalities (3.10), (3.11), and (3.12) could be similarly derived as that of
(3.9), and the details are skipped to save space.

##  Appendix C Butcher-like tableaux of energy-stable methods

We give the coefficients of the energy-stable methods that we used in the
numerical experiments.

  * •

ETD1, also referred to as exponential Euler method:

$\displaystyle\begin{array}[]{>{\centering\arraybackslash$} p{0.4cm} <{$} | >{\centering\arraybackslash$} p{0.6cm} <{$} }0$\@add@centering&0$\@add@centering\cr\hline\cr 1$\@add@centering&\varphi_{1}$\<EMAIL_ADDRESS>(C.3)

  * •

ETDRK2 reads as:

$\displaystyle\begin{array}[]{>{\centering\arraybackslash$} p{0.4cm} <{$} | >{\centering\arraybackslash$} p{1.4cm} <{$} >{\centering\arraybackslash$} p{0.6cm} <{$}}0$\@add@centering&0$\@add@centering&$\@add@centering\\\ 1$\@add@centering&\varphi_{1}$\@add@centering&0$\@add@centering\cr\hline\cr 1$\@add@centering&\varphi_{1}-\varphi_{2}$\@add@centering&\varphi_{2}$\<EMAIL_ADDRESS>(C.7)

  * •

IMEX1, also referred to as first-order semi-implicit method:

$\displaystyle\begin{array}[]{>{\centering\arraybackslash$} p{0.4cm} <{$} | >{\centering\arraybackslash$} p{0.5cm} <{$} >{\centering\arraybackslash$} p{0.5cm} <{$}}0$\@add@centering&0$\@add@centering&0$\@add@centering\cr 1$\@add@centering&0$\@add@centering&1$\@add@centering\cr\hline\cr 1$\@add@centering&0$\@add@centering&1$\@add@centering\end{array},\qquad\begin{array}[]{>{\centering\arraybackslash$} p{0.4cm} <{$} | >{\centering\arraybackslash$} p{0.5cm} <{$} >{\centering\arraybackslash$} p{0.5cm} <{$}}0$\@add@centering&0$\@add@centering&0$\@add@centering\cr 1$\@add@centering&1$\@add@centering&0$\@add@centering\cr\hline\cr 1$\@add@centering&1$\@add@centering&0$\<EMAIL_ADDRESS>(C.14)

  * •

IMEX-RK(2,2) reads as (cf. [37]):

$\displaystyle\begin{array}[]{>{\centering\arraybackslash$} p{0.4cm} <{$} | >{\centering\arraybackslash$} p{0.5cm} <{$} >{\centering\arraybackslash$} p{1.2cm} <{$} >{\centering\arraybackslash$} p{0.5cm} <{$}}0$\@add@centering&0$\@add@centering&0$\@add@centering&0$\@add@centering\cr\gamma$\@add@centering&0$\@add@centering&\gamma$\@add@centering&0$\@add@centering\cr 1$\@add@centering&0$\@add@centering&1-\gamma$\@add@centering&\gamma$\@add@centering\cr\hline\cr 1$\@add@centering&0$\@add@centering&1-\gamma$\@add@centering&\gamma$\@add@centering\end{array},\qquad\begin{array}[]{>{\centering\arraybackslash$} p{0.4cm} <{$} | >{\centering\arraybackslash$} p{0.5cm} <{$} >{\centering\arraybackslash$} p{1.2cm} <{$} >{\centering\arraybackslash$} p{0.5cm} <{$}}0$\@add@centering&0$\@add@centering&0$\@add@centering&0$\@add@centering\cr\gamma$\@add@centering&\gamma$\@add@centering&0$\@add@centering&0$\@add@centering\cr 1$\@add@centering&\delta$\@add@centering&1-\delta$\@add@centering&0$\@add@centering\cr\hline\cr 1$\@add@centering&\delta$\@add@centering&1-\delta$\@add@centering&0$\@add@centering\end{array},$ (C.23)

where $\gamma=\frac{2+\sqrt{2}}{2}$ and $\delta=\frac{2\gamma-1}{2\gamma}$.

All the above schemes can preserve the original energy dissipation property
for gradient flows including SH equation, see, e.g. [15, 19, 38].