# Unconditionally stable exponential integrator schemes for the 2D Cahn-
Hilliard equation

Xinyu Cheng  X. Cheng, Research Institute of Intelligent Complex Systems,
Fudan University, Shanghai, P.R. China<EMAIL_ADDRESS>
###### Abstract.

Phase field models are gradient flows with their energy naturally dissipating
in time. In order to preserve this property, many numerical schemes have been
well-studied. In this paper we consider a well-known method, namely the
exponential integrator method (EI). In the literature a few works studied
several EI schemes for various phase field models and proved the energy
dissipation by either requiring a strong Lipschitz condition on the nonlinear
source term or certain $L^{\infty}$ bounds on the numerical solutions (maximum
principle). However for phase field models such as the (non-local) Cahn-
Hilliard equation, the maximum principle no longer exists. As a result,
solving such models via EI schemes remains open for a long time. In this paper
we aim to give a systematic approach on applying EI-type schemes to such
models by solving the Cahn-Hilliard equation with a first order EI scheme and
showing the energy dissipation. In fact second order EI schemes can be handled
similarly and we leave the discussion in a subsequent paper. To our best
knowledge, this is the first work to handle phase field models without
assuming any strong Lipschitz condition or $L^{\infty}$ boundedness.
Furthermore, we will analyze the $L^{2}$ error and present some numerical
simulations to demonstrate the dynamics.

##  1\. Introduction

###  1.1. Introduction to the models and historical review

Many mathematical models of physical, biological phenomena can be decribed via
partial differential equations (PDEs). Among PDEs, phase field equations are
models of essential importance in the area of material sciences. In this work
we consider two classic phase field models: Allen-Cahn (AC) and Cahn-Hilliard
(CH) equations. The (AC) model was developed in [2] by Allen and Cahn to study
the competition of crystal grain orientations in an annealing process
separation of different metals in a binary alloy; while the (CH) was
introduced in [7] by Cahn and Hilliard to describe the process of phase
separation of different metals in a binary alloy. These equations are
presented as:

$\left\\{\begin{aligned} &\partial_{t}u=\nu\Delta u-f(u),\ \ \
(x,t)\in\Omega\times(0,\infty),\\\ &u(x,0)=u_{0},\end{aligned}\right.\ $ (AC)

and

$\left\\{\begin{aligned} &\partial_{t}u=\Delta(-\nu\Delta u+f(u)),\ \ \
(x,t)\in\Omega\times(0,\infty),\\\ &u(x,0)=u_{0},\end{aligned}\right.\ $ (CH)

where $u(x,t)$ is a real valued function and values of $u$ in $(-1,1)$
represent a mixture of the two phases, with $-1$ representing the pure state
of one phase and $+1$ representing the pure state of the other phase. Vector
position $x$ is in the spatial domain $\Omega$, which is oftentimes taken to
be two or three dimensional periodic domain and $t$ is the time variable. Here
$\nu$ is a small parameter, occasionally we denote $\varepsilon=\sqrt{\nu}$ to
represent an average distance over which phases mix. The energy term $f(u)$ is
often chosen to be

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

It is well known that, as $\varepsilon\to 0$, the limiting problem of (AC) is
driven under a (mean) curvature flow while the limiting problem of (CH)
follows from the Mullins-Sekerka flow; we refer to [16] for AC and [22], [1]
for CH and a recent work for matrix-valued AC [13]. Both asymptotic and
rigorous analysis are well-studided. This current work will present an idea
about how to approach these limiting behaviors in an appropriate way
numerically.

To be more specific, throughout this paper we choose the spatial domain
$\Omega$ to be the two dimensional $2\pi$-periodic torus
$\mathbb{T}^{2}=(\mathbb{R}/2\pi\mathbb{Z})^{2}$. In fact, our analysis can be
applied to more general settings such as Dirichlet and Neumann boundary
conditions in a bounded domain. However, considering the periodic domain
allows us to use the efficient and accurate Fourier-spectral numerical
algorithms. Beyond this, periodic domain is often appropriate for practical
questions, which involve the formation of micro-structure away from physical
boundaries. As is well-known that both the AC and CH behave as gradient flows,
therefore there energy dissipates in time. Here the associated energy
functional is given by

$\displaystyle E(u)=\int_{\Omega}\left(\frac{1}{2}\nu|\nabla
u|^{2}+F(u)\right)\ dx.$ (1.1)

Assume that $u(x,t)$ is a smooth solution, it is clear that the energy decays
in time: $\frac{d}{dt}E(u(t))\leq 0$. Such property thus provides an $a\
priori$ $H^{1}$-norm bound and since the scaling-critical space for (CH) is
$L^{2}$ in 2D (and $H^{\frac{1}{2}}$ in 3D), the global well-posedness follow
from standard energy estimates. Therefore from the analysis point of view, the
energy dissipation property is an important index for whether a numerical
scheme is “stable” or not.

There are various approaches existing in the literature that are developed to
study numerical behaviors of the Cahn-Hilliard and other related phase field
models, cf. [10, 15, 11, 14, 8, 4, 12, 3, 17]. Among which different time
stepping approaches are applied including the fully explicit (forward Euler)
scheme, fully implicit (backward Euler) scheme, finite element scheme and
convex splitting scheme; moreover various schemes are adopted for the spatial
discretization including the Fourier-spectral method and finite difference
schemes. To guarantee the accuracy and stability, numerical approximations
usually need to obey certain qualitative behaviors and features. In this paper
we focus on one of the key features, namely the energy dissipation or
conservation as mentioned earlier.

Based on the discussion above, we briefly go through the related results on
the energy stability in the literature. To start with, Feng and Prohl [14]
introduced a semi-discrete in time and fully spatially discrete finite element
method for CH where they obtained an error bound of size of powers of $1/\nu$.
However, explicit time-stepping schemes usually require strict time-step
restrictions and do not obey energy decay in general (one can compare with
Remark 1.4). To guarantee the energy dissipation with bigger time steps, a
good alternative is to use semi-implicit schemes in which the linear term is
implicit (such as backward time differentiation) and the nonlinear term is
treated explicitly. Having only a linear implicit at every time step has
computational advantages, as suggested in [8], where a semi-implicit Fourier-
spectral scheme for CH was considered. On the other hand, semi-implicit
schemes can lose stability for large time steps and thus smaller time steps
are needed in practice. To resolve this problem, semi-implicit methods with
better stability have been introduced, cf. [19, 20, 9]. Motivated by the work
on semi-implicit schemes, the so-called “exponential integrator schemes” (EI
schemes) start to draw attention in the study of PDEs, cf. [18, 23]. In a
short word, the EI schemes discretize the semi-group presentation of the
solutions to nonlinear PDEs. In more details, take the Allen-Cahn equation as
an example. It is well-known that the smooth solution to (AC) is given by the
semi-group mild form (or Duhamel’s Formula):

$u(t)=e^{\nu t\Delta}u(t_{0})+\int_{t_{0}}^{t}e^{\nu(t-s)\Delta}f(u(s))\ ds,$
(1.2)

where $e^{t\Delta}$ is the standard heat semi-group (heat kernel). Then the
idea of the classic EI is to discretize the semi-group presentation in each
time interval $(t_{n},t_{n+1})$ by approximating the integrand function by a
constant (cf.[18]):

$\displaystyle u(t_{n+1})=$ $\displaystyle
e^{\nu\tau\Delta}u(t_{n})+\int_{t_{n}}^{t_{n+1}}e^{\nu(t_{n+1}-s)\Delta}f(u(s))\
ds$ $\displaystyle\approx$ $\displaystyle
e^{\nu\tau\Delta}u(t_{n})+\int_{t_{n}}^{t_{n+1}}e^{\nu(t_{n+1}-s)\Delta}f(u(t_{n}))\
ds$ $\displaystyle\approx$ $\displaystyle
e^{\nu\tau\Delta}u(t_{n})+\frac{1-e^{\nu(t_{n+1}-t_{n})\Delta}}{\nu\Delta}f(u(t_{n})),$

or equivalently writing $u(t_{n})=u^{n},u(t_{n+1})=u^{n+1}$ and
$\tau=t_{n+1}-t_{n}$ we have

$\displaystyle
u^{n+1}=e^{\nu\tau\Delta}u^{n}+\frac{1-e^{\nu\tau\Delta}}{\nu\Delta}f(u^{n}).$

However, those works either require a strong Lipschitz condition on the
nonlinear source term, or require certain $L^{\infty}$ bounds on the numerical
solutions. Such requirement is satisfied by the Allen-Cahn equation by the
standard maximum principle (we refer the readers to [19, 9] for example).
While for other phase field models such as the (non-local) Cahn-Hilliard
equation, the maximum principle no longer exists. As a result, solving phase
field models via EI schemes where strong Lipschitz conditions or maximum
principle are no longer applicable remains open for a long time. One of our
modest goals of this paper is to give a systematic approach on applying EI-
type schemes to such models. In particular, we consider one of the classic
phase models, the Cahn-Hilliard equation with a first order EI scheme. In fact
second order EI schemes can be handled similarly and we leave the discussion
in a forthcoming paper.

In this work we consider the following stabilized EI scheme for (CH):

$\left\\{\begin{aligned}
&u^{n+1}=e^{-\nu\Delta^{2}\tau}u^{n}+\frac{1-e^{-\nu\Delta^{2}\tau}}{\nu\Delta}\Pi_{N}[f(u^{n})]-S\tau\Delta^{2}(u^{n+1}-u^{n}),\\\
&u^{0}=\Pi_{N}u_{0},\end{aligned}\right.$ (1.3)

where $\tau$ is the time step and $S>0$ is the coefficient for the $O(\tau)$
regularization term. In fact (1.3) above can be reorganized and rewritten into
the following form:

$\left\\{\begin{aligned}
&\frac{u^{n+1}-u^{n}}{\tau}+S\Delta^{2}(u^{n+1}-u^{n})=\frac{e^{-\nu\Delta^{2}\tau}-1}{\tau}u^{n}+\frac{1-e^{-\nu\Delta^{2}\tau}}{\tau\nu\Delta}\Pi_{N}[f(u^{n})],\\\
&u^{0}=\Pi_{N}u_{0}.\end{aligned}\right.$ (1.4)

######  Remark 1.1.

It is worth mentioning that from (1.4) one can observe that our exponential
integrator scheme behaves very similarly to an explicit forward Euler method
since

$\displaystyle\frac{e^{-\nu\Delta^{2}\tau}-1}{\tau}=-\nu\Delta^{2}+\frac{1}{2}\nu^{2}\Delta^{4}\tau+O(\tau^{2}),$
$\displaystyle\frac{1-e^{-\nu\Delta^{2}\tau}}{\tau\nu\Delta}=\Delta-\frac{1}{2}\nu\Delta^{3}\tau+O(\tau^{2}).$

Therefore informally speaking, (1.4) can be understood as the usual explicit
forward Euler scheme for (CH) with a stabilization term:

$\displaystyle\frac{u^{n+1}-u^{n}}{\tau}+S\Delta^{2}(u^{n+1}-u^{n})=-\nu\Delta^{2}u^{n}+\Delta\Pi_{N}[f(u^{n})]+O(\nu\tau,\tau^{2}).$
(1.5)

######  Remark 1.2.

The stabilizer $S\Delta^{2}(u^{n+1}-u^{n})$ is not unique. In fact our
analysis works for many other stabilizers such as
$Se^{-\Delta}(u^{n+1}-u^{n})$ and so on. Another particular choice is
$\frac{\nu}{2}\Delta^{2}(u^{n+1}-u^{n})-S\Delta(u^{n+1}-u^{n})$. Such
particular choice results in the analysis of the energy dissipation. More
precisely speaking, we need $\frac{\nu}{2}\Delta^{2}(u^{n+1}-u^{n})$ to
balance the drawback caused by the “explicit forward” diffusion; the real
stabilizer is indeed $-S\Delta(u^{n+1}-u^{n})$. We choose
$S\Delta^{2}(u^{n+1}-u^{n})$ just for the sake of simplicity. Even though such
choices of stabilizers are not unique, they are necessary. As will be
presented in Section 5 with numerical evidence supporting, we will see
examples without stabilizers have growing energy.

######  Remark 1.3.

The exact explicit forward Euler method seems challenging to obtain the energy
stability because it seems unlikely to obtain a uniform Sobolev $H^{s}$-bound
estimate. To our best knowledge, we believe such problem arises from the
essence of the forward Euler method since the regularity at step $n$ is higher
than the regularity at step $n+1$ due to the “explicit forward” diffusion
$-\Delta^{2}u^{n}$. This seems to be the crucial technical difference between
our exponential integrator scheme (1.4) and the exact explicit forward Euler
scheme (1.5) without the stabilizer $S\Delta^{2}(u^{n+1}-u^{n})$ and higher
order $O(\nu\tau,\tau^{2})$ perturbation.

######  Remark 1.4.

As mentioned earlier in Remark 1.1, our scheme (1.4) can be understood as an
explicit forward Euler scheme. It is well known that after spatial
discretization on a grid of size $h$, any explicit scheme will have a time
step $\tau$ restriction of size $\tau<O(h^{4})$ for the fourth order CH
equation. However as will be explained later (see Theorem 1.7), in order to
show the convergence of our numerical scheme, we only require
$\tau\ll(\log(N))^{-1}$, which is a very mild restriction. Moreover the energy
stability holds for any size of $\tau$, see Theorem 1.5.

###  1.2. Main results

Our main results state below.

######  Theorem 1.5 (Unconditional energy stability for CH).

Consider (1.3) or (1.4) with $\nu>0$ and assume $u_{0}\in
H^{2}(\mathbb{T}^{2})$. Then there exists a constant $\beta_{0}$ depending
only on the initial energy $E_{0}=E(u_{0})$ such that if

$\begin{aligned}
S\geq\beta\cdot\left(\|u_{0}\|^{2}_{H^{2}}+\nu^{-1}|\log\nu|+\nu^{-1}\log
N+\nu\right)\end{aligned},\quad\beta\geq\beta_{0}.$

Then $E(u^{n+1})\leq E(u^{n})$, $\forall n\geq 0$ and for any choice of the
time step $\tau$, where the energy $E$ is defined in (1.1).

######  Remark 1.6.

The choice of $S$ is not optimal here; in fact we refer the readers to Section
5 for numerical evidence with much smaller $S$. Determining the optimal bound
of the stabilizers can be a very interesting question, however it is not the
focus of the presenting paper; we leave this question to the readers.

######  Theorem 1.7 ($L^{2}$ error estimate).

Fix $\nu>0$ and let $u_{0}\in H^{s}$, $s\geq 10$. Let $0<\tau\leq M$ for some
$M>0$. Let $u(t)$ be the continuous solution to the 2D Cahn-Hilliard equation
with initial data $u_{0}$. Let $u^{m}$, $m\geq 1$ be defined in (1.4) with
initial data $u^{0}$. Assume $S$ satisfies the same condition in Theorem 1.5.
Define $t_{0}=0$ and $t_{m}=m\tau$ for $m\geq 1$. Then for any $m\geq 1$,

$\|u^{m}-u(t_{m})\|_{2}\leq(1+S)\cdot C_{2}\cdot
e^{C_{1}t_{m}}\left(N^{-s}+\tau+\tau\cdot N^{-s+4}\right),$

where $C_{1}>0$ is a constant depending on $\nu\ ,u_{0}$; $C_{2}>0$ is a
constant depending on $s\ ,\nu$ and $u_{0}$.

######  Remark 1.8.

It is worth mentioning that the requirement $S\geq\beta\nu^{-1}\log N$ in
Theorem 1.5 is troublesome in the convergence of the numerical scheme. However
as stated in Theorem 1.7 above, $S\cdot N^{-s}\to 0$ as $N\to\infty$ therefore
such choice of $S$ will not ruin the convergence. On the other hand, we do
need $\tau\cdot S\ll 1$, which means $\tau\ll(\log N)^{-1}$. This is still a
very mild restriction as discussed in Remark 1.4 because if one uses spatial
discretization, the restriction on the size of the time step $\tau$ is usually
much smaller (of polynomial order of $\frac{1}{N}$) to be numerically stable.

######  Remark 1.9.

Here we need to require that $\tau$ is not arbitrarily large. However, in
practice it is not a big issue as we always use small (or at least not
arbitrarily large) time steps.

###  1.3. Organization of the presenting paper

The presenting paper is organized as follows. In Section 2 we list the
notation and preliminaries including several useful lemmas. The energy
stability of the EI scheme of the 2D Cahn-Hilliard will be shown in Section 3
while the error estimate is given in Section 4. Numerical experiments will be
presented in Section 5. Further discussion and concluding remarks will be in
Section 6.

##  2\. Notation and preliminaries

Throughout this paper, for any two (non-negative in particular) quantities $X$
and $Y$, we denote $X\lesssim Y$ if $X\leq CY$ for some constant $C>0$.
Similarly $X\gtrsim Y$ if $X\geq CY$ for some $C>0$. We denote $X\sim Y$ if
$X\lesssim Y$ and $Y\lesssim X$. The dependence of the constant $C$ on other
parameters or constants are usually clear from the context and we will often
suppress this dependence. We shall denote
$X\lesssim_{Z_{1},Z_{2},\cdots,Z_{k}}Y$ if $X\leq CY$ and the constant $C$
depends on the quantities $Z_{1},\cdots,Z_{k}$.

For any two quantities $X$ and $Y$, we shall denote $X\ll Y$ if $X\leq cY$ for
some sufficiently small constant $c$. The smallness of the constant $c$ is
usually clear from the context. The notation $X\gg Y$ is similarly defined.
Note that our use of $\ll$ and $\gg$ here is _different_ from the usual
Vinogradov notation in number theory or asymptotic analysis.

For a real-valued function $u:\Omega\to\mathbb{R}$ we denote its usual
Lebesgue $L^{p}$-norm by

$\displaystyle\|u\|_{p}=\|u\|_{L^{p}(\Omega)}=\begin{cases}&\left(\int_{\Omega}|u|^{p}\
dx\right)^{\frac{1}{p}},\quad 1\leq p<\infty;\\\
&\operatorname{esssup}_{x\in\Omega}|u(x)|,\quad p=\infty.\end{cases}$ (2.1)

Similarly, we use the weak derivative in the following sense: For $u$, $v\in
L^{1}_{loc}(\Omega)$, (i.e they are locally integrable); $\forall\phi\in
C^{\infty}_{0}(\Omega)$, i.e $\phi$ is infinitely differentiable (smooth) and
compactly supported; and

$\int_{\Omega}u(x)\ \partial^{\alpha}\phi(x)\
dx=(-1)^{\alpha_{1}+\cdots+\alpha_{n}}\int_{\Omega}v(x)\ \phi(x)\ dx,$

then $v$ is defined to be the weak partial derivative of $u$, denoted by
$\partial^{\alpha}u$. Suppose $u\in L^{p}(\Omega)$ and all weak derivatives
$\partial^{\alpha}u$ exist for $|\alpha|=\alpha_{1}+\cdots+\alpha_{n}\leq k$ ,
such that $\partial^{\alpha}u\in L^{p}(\Omega)$ for $|\alpha|\leq k$, then we
denote $u\in W^{k,p}(\Omega)$ to be the standard Sobolev space. The
corresponding norm of $W^{k,p}(\Omega)$ is :

$\|u\|_{W^{k,p}(\Omega)}=\left(\sum_{|\alpha|\leq
k}\int_{\Omega}|\partial^{\alpha}u|^{p}\ dx\right)^{\frac{1}{p}}\ .$

For $p=2$ case, we use the convention $H^{k}(\Omega)$ to denote the space
$W^{k,2}(\Omega)$. We often use $D^{m}u$ to denote any differential operator
$D^{\alpha}u$ for any $|\alpha|=m$: $D^{2}$ denotes
$\partial^{2}_{x_{i}x_{j}}u$ for $1\leq i,j\leq d$, as an example.

In this paper we use the following convention for Fourier expansion on
$\mathbb{T}^{d}$:

$f(x)=\frac{1}{(2\pi)^{d}}\sum_{k\in\mathbb{Z}^{d}}\widehat{f}(k)e^{ik\cdot
x}\ ,\ \widehat{f}(k)=\int_{\Omega}f(x)e^{-ik\cdot x}\ dx\ .$

Taking advantage of the Fourier expansion, we use the well-known equivalent
$H^{s}$-norm and $\dot{H}^{s}$-semi-norm of function $f$ by

$\|f\|_{H^{s}}=\frac{1}{(2\pi)^{d/2}}\left(\sum_{k\in\mathbb{Z}^{d}}(1+|k|^{2s})|\widehat{f}(k)|^{2}\right)^{\frac{1}{2}}\
,\
\|f\|_{\dot{H}^{s}}=\frac{1}{(2\pi)^{d/2}}\left(\sum_{k\in\mathbb{Z}^{d}}|k|^{2s}|\widehat{f}(k)|^{2}\right)^{\frac{1}{2}}.$

In addition, for $N\geq 2$, we define

$X_{N}=\text{span}\left\\{\cos(k\cdot x)\ ,\sin(k\cdot x):\ \
k=(k_{1},k_{2})\in\mathbb{Z}^{2}\ ,\
|k|_{\infty}=\max\\{|k_{1}|,|k_{2}|\\}\leq N\right\\}\ .$

######  Lemma 2.1 (Sobolev inequality on $\mathbb{T}^{d}$).

Let $0<s<d$ and $f\in L^{q}(\mathbb{T}^{d})$ for any $\frac{d}{d-s}<p<\infty$,
then

$\displaystyle\|\left<\nabla\right>^{-s}f\|_{L^{p}(\mathbb{T}^{d})}\lesssim_{s,p,d}\|f\|_{L^{q}(\mathbb{T}^{d})}\
,\ \mbox{where}\ \frac{1}{q}=\frac{1}{p}+\frac{s}{d}\ ,$

where $\left<\nabla\right>^{-s}$ denotes $(1-\Delta)^{-\frac{s}{2}}$ and
$A\lesssim_{s,p,d}B$ is defined as $A\leq C_{s,p,d}\ B$ where $C_{s,p,d}$ is a
constant dependent on $s,p$ and $d$.

######  Remark 2.2.

Note that the this Sobolev inequality is a variety of the standard version.
Note that on the Fourier side the symbol of $\left<\nabla\right>^{-s}$ is
given by $(1+|k|^{2})^{-\frac{s}{2}}$. In particular,
$\|f\|_{\infty(\mathbb{T}^{d})}\lesssim\|f\|_{H^{2}(\mathbb{T}^{d})}$, known
as Morrey’s inequality.

######  Remark 2.3.

If one further requires that $f\in H^{1}(\mathbb{T}^{2})$ has zero mean, i.e.
$\widehat{f}(0)=0$ we have

$\displaystyle\|f\|_{4}\leq\|f\|_{2}^{\frac{1}{2}}\|\nabla
f\|_{2}^{\frac{1}{2}},$ (2.2)

and

$\displaystyle\|f\|_{6}\leq\|f\|_{2}^{\frac{1}{3}}\|\nabla
f\|_{2}^{\frac{2}{3}}.$ (2.3)

The proof of (2.2) and (2.3) follows from Lemma 2.1 and standard
interpolation.

######  Lemma 2.4.

Suppose $f\in H^{1}(\mathbb{T}^{2})$ and $f$ has zero mean, i.e.
$\widehat{f}(0)=0$. Then

$\displaystyle\|f\|_{2}\leq\||\nabla|^{-1}f\|_{2}^{\frac{1}{2}}\|\nabla
f\|_{2}^{\frac{1}{2}},$ (2.4) $\displaystyle\|f\|_{2}\leq\|\nabla f\|_{2}.$
(2.5)

###### Proof.

The proof follows from Parseval’s identity directly:

$\displaystyle\int_{\mathbb{T}^{2}}f\
dx=\int_{\mathbb{T}^{2}}|\nabla|^{-1}f\cdot|\nabla|f\ dx.$

Then (2.4) follows from the classic Cauchy-Schwarz inequality by noting
$\||\nabla|f\|_{2}=\|\nabla f\|_{2}$. (2.5) follows from the standard Poincaré
inequality (or directly from the Fourier side). ∎

######  Lemma 2.5 (Log-type interpolation).

For all $f\in H^{s}(\mathbb{T}^{2})\ ,\ s>1$ and suppose $f$ has zero mean,
i.e. $\widehat{f}(0)=0$, then

$\|f\|_{\infty}\leq
C_{s}\cdot\left(\|f\|_{\dot{H}^{1}}\sqrt{\log(\|f\|_{\dot{H}^{s}}+3)}+1\right)\
.$

Here $C_{s}$ is a constant which only depends on s.

###### Proof.

This lemma is a special case of Lemma 3.1 in [9] and we only sketch the proof
here. Note that $f$ is mean-zero, we then consider the Fourier series
$f(x)=\frac{1}{(2\pi)^{2}}\sum_{|k|\geq 1}\widehat{f}(k)\ e^{ik\cdot x}$. It
then follows that

$\displaystyle\|f\|_{\infty}$
$\displaystyle\leq\frac{1}{(2\pi)^{2}}\sum_{|k|\geq 1}|\widehat{f}(k)|$
$\displaystyle\leq\frac{1}{(2\pi)^{2}}\left(\sum_{1\leq|k|\leq
N}|\widehat{f}(k)|+\sum_{|k|>N}|\widehat{f}(k)|\right)$
$\displaystyle\lesssim\sum_{1\leq|k|\leq
N}(|\widehat{f}(k)\|k|\cdot|k|^{-1})+\sum_{|k|>N}(|\widehat{f}(k)\|k|^{s}\cdot|k|^{-s})$
$\displaystyle\lesssim\left(\sum_{1\leq|k|\leq
N}|\widehat{f}(k)|^{2}|k|^{2}\right)^{\frac{1}{2}}\cdot\left(\sum_{1\leq|k|\leq
N}|k|^{-2}\right)^{\frac{1}{2}}+\left(\sum_{|k|>N}|\widehat{f}(k)|^{2}|k|^{2s}\right)^{\frac{1}{2}}\cdot(\sum_{|k|>N}|k|^{-2s})^{\frac{1}{2}}$
$\displaystyle\lesssim\frac{1}{N^{s-1}}\left(\sum_{|k|>N}|\widehat{f}(k)|^{2}|k|^{2s}\right)^{\frac{1}{2}}+\left(\sum_{1\leq|k|\leq
N}|\widehat{f}(k)|^{2}|k|^{2}\right)^{\frac{1}{2}}\cdot\sqrt{\log(N+3)}$
$\displaystyle\lesssim\frac{1}{N^{s-1}}\|f\|_{\dot{H}^{s}}+\sqrt{\log(N+3)}\|f\|_{\dot{H}^{1}}\
.$

If $\|f\|_{\dot{H}^{s}}\leq 3$, we can simply take $N=1$; otherwise take
$N^{s-1}$ close to $\|f\|_{\dot{H}^{s}}$ . ∎

######  Lemma 2.6 (Discrete Grönwall’s inequality).

Let $\tau>0$ and $y_{n}\geq 0$, $\alpha_{n}\geq 0$, $\beta_{n}\geq 0$ for
$n=1,2,3\cdots$. Suppose

$\frac{y_{n+1}-y_{n}}{\tau}\leq\alpha_{n}y_{n}+\beta_{n}\ ,\ \forall\ n\geq 0\
.$

Then for any $m\geq 1$, we have

$\displaystyle
y_{m}\leq\exp\left(\tau\sum_{n=0}^{m-1}\alpha_{n}\right)\left(y_{0}+\sum_{k=0}^{m-1}\beta_{k}\right)\
.$

###### Proof.

We only sketch the proof here for the sake of completeness. By the assumption,
it follows that for $n\geq 0$,

$y_{n+1}\leq(1+\alpha_{n}\tau)y_{n}+\tau\beta_{n}\leq
e^{\tau\alpha_{n}}y_{n}+\tau\beta_{n};$

therefore one can derive that

$\exp\left(-\tau\sum_{j=0}^{n}\alpha_{j}\right)y_{n+1}\leq\exp\left(-\tau\sum_{j=0}^{n-1}\alpha_{j}\right)y_{n}+\exp\left(-\tau\sum_{j=0}^{n}\alpha_{j}\right)\beta_{n}.$

We thus obtain the desired result by performing a telescoping summation. ∎

######  Lemma 2.7 ($H^{k}$ boundedness of the exact CH solution).

Assume $u(x,t)$ is a smooth solution to the Cahn-Hilliard equation in
$\mathbb{T}^{2}$ and the initial data $u_{0}\in H^{k}(\mathbb{T}^{2})$ for
$k\geq 2$. Then,

$\sup_{t\geq 0}\|u(t)\|_{H^{k}(\mathbb{T}^{2})}\lesssim_{k}1$ (2.6)

where we omit the dependence on $\nu$ and $u_{0}$.

###### Proof.

The proof is standard and we postpone the proof to the Appendix. ∎

##  3\. Energy stability of the scheme (1.4) on the 2D Cahn-Hilliard equation

Recall that the Cahn-Hilliard equation (CH) takes the following form:

$\displaystyle\begin{cases}\partial_{t}u=-\nu\Delta^{2}u+\Delta f(u),\\\
u(x,0)=u_{0}.\end{cases}$

Here $f(u)=u^{3}-u$, and the spatial domain $\Omega$ is taken to be the two
dimensional $2\pi-$periodic torus $\mathbb{T}^{2}$. The corresponding energy
is defined by $E(u)=\int_{\Omega}(\frac{\nu}{2}|\nabla u|^{2}+F(u))\ dx$ ,
where $F(u)=\frac{1}{4}(u^{2}-1)^{2}$, the anti-derivative of $f(u)$. Recall
that we consider the stabilized exponential integrator scheme (1.4):

$\left\\{\begin{aligned}
&\frac{u^{n+1}-u^{n}}{\tau}+S\Delta^{2}(u^{n+1}-u^{n})=\frac{e^{-\nu\Delta^{2}\tau}-1}{\tau}u^{n}+\frac{1-e^{-\nu\Delta^{2}\tau}}{\tau\nu\Delta}\Pi_{N}f(u^{n}),\\\
&u^{0}=\Pi_{N}u_{0}.\end{aligned}\right.$ (3.1)

Here $\Pi_{N}$ is truncation of Fourier modes of $L^{2}$ functions to
$|k|_{\infty}\leq N$. Note that here $E^{0}=E(\Pi_{N}u_{0})$ while
$E_{0}=E(u_{0})$ and in general $E_{0}\neq E^{0}$. In particular the following
statement holds.

######  Lemma 3.1.

Suppose $E^{0}=E(\Pi_{N}u_{0})$ and $E_{0}=E(u_{0})$ as defined above, the
following inequality holds:

$\sup_{N}\ E(\Pi_{N}u_{0})\lesssim 1+E_{0},\mbox{where}\ u_{0}\in
H^{1}({\mathbb{T}^{2}})\ .$

###### Proof.

We rewrite $\Pi_{N}u_{0}$ as $\frac{1}{(2\pi)^{2}}\sum_{|k|\leq
N}\widehat{u_{0}}(k)e^{ik\cdot x}$ , namely the Dirichlet partial sum of
$u_{0}$.

$\begin{aligned}
\|\nabla\left(\Pi_{N}u_{0}\right)\|_{L^{2}(\mathbb{T}^{2})}^{2}=\frac{1}{(2\pi)^{2}}\sum_{|k|\leq
N}|k|^{2}|\widehat{u_{0}}(k)|^{2}\leq\frac{1}{(2\pi)^{2}}\sum_{|k|\in\mathbb{Z}^{2}}|k|^{2}|\widehat{u_{0}}(k)|^{2}=\|\nabla\left(u_{0}\right)\|_{L^{2}(\mathbb{T}^{2})}^{2}\end{aligned}\
.$

On the potential energy part, by the Sobolev inequality Lemma 2.1,
$\|u_{0}\|_{L^{4}(\mathbb{T}^{2})}\lesssim\|u_{0}\|_{H^{1}(\mathbb{T}^{2})}$,
this shows $u_{0}\in L^{4}(\mathbb{T}^{2})$ and hence the Dirichlet partial
sum $\Pi_{N}u_{0}$ converges to $u_{0}$ in $L^{4}(\mathbb{T}^{2})$. Then
$\lim_{N}\|\Pi_{N}u_{0}\|_{L^{4}(\mathbb{T}^{2})}=\|u_{0}\|_{L^{4}(\mathbb{T}^{2})}$,
which leads to $\sup_{N}\|\Pi_{N}u_{0}\|_{L^{4}(\mathbb{T}^{2})}<\infty$. By
the Uniform Boundedness Principle, we derive $\sup_{N}\|\Pi_{N}\|<\infty$,
i.e. $\sup_{N}\|\Pi_{N}u_{0}\|_{L^{4}(\mathbb{T}^{2})}\leq
c\|u_{0}\|_{L^{4}(\mathbb{T}^{2})}$ for an absolute constant $c$. Combining
the two estimates above we finish the proof. ∎

We will prove Theorem 1.5 by induction. To start with, let us multiply the
equation by $\frac{\nu\Delta\tau}{e^{-\nu\tau\Delta^{2}}-1}(u^{n+1}-u^{n})$
and integrate over $\mathbb{T}^{2}$. Then the LHS of (3.1) is

$\displaystyle\mbox{LHS}=$
$\displaystyle\frac{1}{\tau}\left(u^{n+1}-u^{n},\frac{\nu\Delta\tau}{e^{-\nu\tau\Delta^{2}}-1}(u^{n+1}-u^{n})\right)+S\left(\Delta^{2}(u^{n+1}-u^{n}),\frac{\nu\Delta\tau}{e^{-\nu\tau\Delta^{2}}-1}(u^{n+1}-u^{n})\right)$
(3.2) $\displaystyle=$
$\displaystyle\left(u^{n+1}-u^{n},\frac{-\nu\Delta}{1-e^{-\nu\tau\Delta^{2}}}(u^{n+1}-u^{n})\right)+S\left(u^{n+1}-u^{n},\frac{-\nu\Delta^{3}\tau}{1-e^{-\nu\tau\Delta^{2}}}(u^{n+1}-u^{n})\right),$
(3.3)

where $(\cdot,\cdot)$ is the usual $L^{2}$ inner product. We observe that both
the two operators $A=\frac{-\nu\Delta}{1-e^{-\nu\tau\Delta^{2}}}$ and
$B=\frac{-\nu\Delta^{3}\tau}{1-e^{-\nu\tau\Delta^{2}}}$ are positive semi-
definite and self-adjoint, therefore we have

$\displaystyle\mbox{LHS}=\|\sqrt{A}(u^{n+1}-u^{n})\|_{2}^{2}+S\|\sqrt{B}(u^{n+1}-u^{n})\|_{2}^{2},$
(3.4)

where the operator $\sqrt{A}$ and $\sqrt{B}$ correspond to the following
Fourier symbols:

$\displaystyle\widehat{\sqrt{A}u}(k)=\sqrt{\frac{\nu|k|^{2}}{1-e^{-\nu|k|^{4}\tau}}}\widehat{u}(k);\quad\widehat{\sqrt{B}u}(k)=\sqrt{\frac{\nu|k|^{6}\tau}{1-e^{-\nu|k|^{4}\tau}}}\widehat{u}(k).$
(3.5)

To proceed from here, it is important to examine the behaviors of these two
operators. Indeed we state the following lemma.

######  Lemma 3.2.

Suppose the operators $\sqrt{A}$ and $\sqrt{B}$ are defined in (3.5) above.
Then the following inequalities hold for any function $u\in
H^{3}(\mathbb{T}^{2})$ and any $v\in H^{1}(\mathbb{T}^{2})$ with zero-mean:

$\displaystyle\left(\frac{1}{\tau}\right)^{\frac{1}{2}}\||\nabla|^{-1}v\|_{2}\leq\|\sqrt{A}v\|_{2},\quad\|\nabla
u\|_{2}\leq\|\sqrt{B}u\|_{2}.$ (3.6)

###### Proof.

The (3.6) can be proved from the Fourier side. In particular, by the
Parseval’s identity we have

$\displaystyle\|\sqrt{A}v\|_{2}^{2}=\frac{1}{(2\pi)^{2}}\sum_{|k|\geq
1}\frac{\nu|k|^{2}}{1-e^{-\nu|k|^{4}\tau}}|\widehat{v}(k)|^{2},$ (3.7)

therefore it suffices to show

$\displaystyle\frac{\nu|k|^{2}}{1-e^{-\nu|k|^{4}\tau}}\geq\frac{1}{\tau|k|^{2}},\quad\mbox{or}\quad\frac{\nu\tau|k|^{4}}{1-e^{-\nu\tau|k|^{4}}}\geq
1.$ (3.8)

To show (3.8), we can define an 1D function for $x\geq 1$, $h(x)=\frac{\nu\tau
x^{4}}{1-e^{-\nu\tau x^{4}}}-1$. By direct computation we can show
$h^{\prime}(x)>0$ for $x\geq 1$ and $h(1)>0$ noting $\nu\tau>0$. Therefore
(3.8) is proved. Then (3.7) implies

$\displaystyle\frac{1}{\tau}\||\nabla|^{-1}v\|_{2}^{2}\leq\|\sqrt{A}v\|_{2}^{2}.$
(3.9)

By the same reason we can show $\|\nabla u\|_{2}\leq\|\sqrt{B}u\|_{2}$. We
hereby conclude (3.6). ∎

We then continue the proof of Theorem 1.5. We now focus on the RHS of (3.1):

$\displaystyle\mbox{RHS}=$
$\displaystyle\left(\frac{e^{-\nu\Delta^{2}\tau}-1}{\tau}u^{n},\frac{-\nu\Delta\tau}{1-e^{-\nu\tau\Delta^{2}}}(u^{n+1}-u^{n})\right)+\left(\frac{1-e^{-\nu\Delta^{2}\tau}}{\nu\tau\Delta}\Pi_{N}f(u^{n}),\frac{-\nu\Delta\tau}{1-e^{-\nu\tau\Delta^{2}}}(u^{n+1}-u^{n})\right)$
(3.10) $\displaystyle=$
$\displaystyle\nu\left(u^{n},\Delta(u^{n+1}-u^{n})\right)-\left(\Pi_{N}f(u^{n}),u^{n+1}-u^{n}\right).$
(3.11)

Note that $u^{n+1},u^{n}\in X_{N}$, we have

$\displaystyle\mbox{RHS}=\nu\left(\nabla
u^{n},\nabla(u^{n}-u^{n+1})\right)-\left(f(u^{n}),u^{n+1}-u^{n}\right).$
(3.12)

Using the identity $2a(a-b)=a^{2}-b^{2}+(a-b)^{2}$, we derive from (3.12) that

$\displaystyle\mbox{RHS}=\frac{\nu}{2}(\|\nabla u^{n}\|_{2}^{2}-\|\nabla
u^{n+1}\|_{2}^{2}+\|\nabla u^{n}-\nabla
u^{n+1}\|_{2}^{2})-(f(u^{n}),u^{n+1}-u^{n}).$ (3.13)

To proceed further, by the fundamental theorem of calculus we arrive at

$\displaystyle F(u^{n+1})-F(u^{n})$
$\displaystyle=f(u^{n})(u^{n+1}-u^{n})+\int_{u^{n}}^{u^{n+1}}f^{\prime}(s)(u^{n+1}-s)\
ds$ (3.14)
$\displaystyle=f(u^{n})(u^{n+1}-u^{n})+\int_{u^{n}}^{u^{n+1}}(3s^{2}-1)(u^{n+1}-s)\
ds$
$\displaystyle=f(u^{n})(u^{n+1}-u^{n})+\frac{1}{4}(u^{n+1}-u^{n})^{2}\left(3(u^{n})^{2}+(u^{n+1})^{2}+2u^{n}u^{n+1}-2\right)\
.$

Combining the estimates (3.4), (3.13) and (3.14) we arrive at

$\displaystyle\|\sqrt{A}(u^{n+1}-u^{n}\|_{2}^{2}+S\|\sqrt{B}(u^{n+1}-u^{n})\|_{2}^{2}+E(u^{n+1})-E(u^{n})$
(3.15) $\displaystyle=$
$\displaystyle\frac{\nu}{2}\|\nabla(u^{n+1}-u^{n})\|_{2}^{2}+\frac{1}{4}\left((u^{n+1}-u^{n})^{2},3(u^{n})^{2}+(u^{n+1})^{2}+2u^{n}u^{n+1}-2\right).$

By further considering (3.6) and (3.15) we indeed have

$\displaystyle\frac{1}{\tau}\||\nabla|^{-1}(u^{n+1}-u^{n})\|_{2}^{2}+(S-\frac{\nu}{2})\|\nabla(u^{n+1}-u^{n})\|_{2}^{2}+E(u^{n+1})-E(u^{n})$
(3.16)
$\displaystyle\leq\frac{1}{4}\left((u^{n+1}-u^{n})^{2},3(u^{n})^{2}+(u^{n+1})^{2}+2u^{n}u^{n+1}-2\right).$

Moreover by Lemma 2.4, we derive that

$\displaystyle(\sqrt{\frac{2S}{\tau}}+\frac{S}{2}-\frac{\nu}{2})\|u^{n+1}-u^{n}\|_{2}^{2}+E(u^{n+1})-E(u^{n})\leq\frac{1}{4}\left((u^{n+1}-u^{n})^{2},3(u^{n})^{2}+(u^{n+1})^{2}+2u^{n}u^{n+1}-2\right).$
(3.17)

Therefore to show the energy dissipation, it suffices to take

$\begin{aligned}
\sqrt{\frac{2S}{\tau}}+\frac{S}{2}-\frac{\nu}{2}\geq\frac{3}{2}\max\left\\{\|u^{n}\|_{\infty}^{2}\
,\ \|u^{n+1}\|_{\infty}^{2}\right\\}\end{aligned}\ .$ (3.18)

We rewrite the numerical scheme (1.4) as follows:

$\displaystyle
u^{n+1}=\frac{S\tau\Delta^{2}+e^{-\nu\Delta^{2}\tau}}{1+S\tau\Delta^{2}}u^{n}+\frac{1-e^{-\nu\Delta^{2}\tau}}{(1+S\tau\Delta^{2})\nu\Delta}\Pi_{N}f(u^{n})$
(3.19)

By the interpolation lemma (Lemma 2.5), to control $\|u^{n+1}\|_{\infty}$ and
$\|u^{n}\|_{\infty}$, we can only focus on the $\dot{H}^{1}$-norm and
$\dot{H}^{\frac{3}{2}}$-norm.

######  Lemma 3.3.

There is an absolute constant $c_{1}>0$ such that for any $n\geq 0$

$\displaystyle\begin{cases}\|u^{n+1}\|_{\dot{H}^{\frac{3}{2}}(\mathbb{T}^{2})}&\leq
c_{1}\cdot\left(\sqrt{\frac{N}{\nu}}+\frac{1}{S}\right)\cdot(E^{n}+1)\\\
\|u^{n+1}\|_{\dot{H}^{1}(\mathbb{T}^{2})}&\leq\left(\frac{1}{S}+\frac{1}{S}\|u^{n}\|^{2}_{\infty}\right)\cdot\|u^{n}\|_{L^{2}(\mathbb{T}^{2})}+\|u^{n}\|_{\dot{H}^{1}(\mathbb{T}^{2})}\
.\end{cases}$ (3.20)

###### Proof.

It suffices to consider Fourier modes $1\leq|k|\leq N$ from the Fourier side.
We then obtain that

$\left\\{\begin{aligned}
&\frac{(S\tau|k|^{4}+e^{-\nu\tau|k|^{4}})|k|^{\frac{3}{2}}}{1+S\tau|k|^{4}}\leq|k|\cdot
N^{\frac{1}{2}}\\\
&\frac{(1-e^{-\nu\tau|k|^{4}})|k|^{\frac{3}{2}}}{(1+S\tau|k|^{4})\nu|k|^{2}}\leq\frac{1}{S|k|^{\frac{1}{2}}},\end{aligned}\right.$
(3.21)

where we have observed that $1-e^{-\nu\tau|k|^{4}}\leq\nu\tau|k|^{4}$. It then
follows that

$\|u^{n+1}\|_{\dot{H}^{\frac{3}{2}}(\mathbb{T}^{2})}\leq
N^{\frac{1}{2}}\|u^{n}\|_{\dot{H}^{1}(\mathbb{T}^{2})}+\frac{1}{S}\|\big{\langle}\nabla\big{\rangle}^{-\frac{1}{2}}f(u^{n})\|_{L^{2}(\mathbb{T}^{2})}.$
(3.22)

Here the notaion
$\big{\langle}\nabla\big{\rangle}^{s}=(1-\Delta)^{\frac{s}{2}}$ , corresponds
to the Fourier side $(1+|k|^{2})^{s/2}$. Note that by the Sobolev inequality
Lemma 2.1 we have

$\displaystyle\|\big{\langle}\nabla\big{\rangle}^{-\frac{1}{2}}f(u^{n})\|_{L^{2}(\mathbb{T}^{2})}$
$\displaystyle\lesssim\|f(u^{n})\|_{L^{\frac{4}{3}}(\mathbb{T}^{2})}=\|(u^{n})^{3}-u^{n}\|_{L^{\frac{4}{3}}(\mathbb{T}^{2})}$
$\displaystyle=\left(\int_{\mathbb{T}^{2}}((u^{n})^{3}-u^{n})^{\frac{4}{3}}\
dx\right)^{\frac{3}{4}}$ $\displaystyle\lesssim E^{n}+1\ .$

Therefore (3.22) can be estimated by the following:

$\|u^{n+1}\|_{\dot{H}^{\frac{3}{2}}(\mathbb{T}^{2})}\lesssim\left(\sqrt{\frac{N}{\nu}}+\frac{1}{S}\right)(E^{n}+1).$

Similarly, we get

$\left\\{\begin{aligned}
&\frac{(S\tau|k|^{4}+e^{-\nu\tau|k|^{4}})|k|}{1+S\tau|k|^{4}}\leq|k|\\\
&\frac{(1-e^{-\nu\tau|k|^{4}})|k|}{(1+S\tau|k|^{4})\nu|k|^{2}}\leq\frac{\nu\tau|k|^{4}}{S\tau\nu|k|^{6}}|k|=\frac{1}{S|k|}\
.\end{aligned}\right.$

This implies for $|k|\geq 1$ that

$\displaystyle\|u^{n+1}\|_{\dot{H}^{1}(\mathbb{T}^{2})}$
$\displaystyle\leq\|u^{n}\|_{\dot{H}^{1}(\mathbb{T}^{2})}+\frac{1}{S}\|f(u^{n})\|_{L^{2}(\mathbb{T}^{2})}$
$\displaystyle\leq\|u^{n}\|_{\dot{H}^{1}(\mathbb{T}^{2})}+\frac{1}{S}\|(1-(u^{n})^{2})\cdot
u^{n}\|_{L^{2}(\mathbb{T}^{2})}$
$\displaystyle\leq\|u^{n}\|_{\dot{H}^{1}(\mathbb{T}^{2})}+(\frac{1}{S}+\frac{\|u^{n}\|^{2}_{\infty}}{S})\|u^{n}\|_{L^{2}(\mathbb{T}^{2})}.$

∎

###### Proof of Theorem 1.5.

Now we will complete the proof for Theorem 1.5 by induction:

Step 1: The induction $n\to n+1$ step. Assuming $E^{n}\leq
E^{n-1}\leq\cdots\leq E^{0}$ and $E^{n}\leq\sup_{N}\ E(\Pi_{N}u_{0})$, we will
show $E^{n+1}\leq E^{n}$. This implies $\|u^{n}\|^{2}_{\dot{H}^{1}}=\|\nabla
u^{n}\|^{2}_{L^{2}}\leq\frac{2E^{n}}{\nu}\leq\frac{2E^{0}}{\nu}$. By Lemma 2.5
and adopting the notation $f\lesssim_{E^{0}}g$ to denote that $f\leq
C(E^{0})\cdot g$ for some constant $C(E^{0})$ depending only on $E^{0}$, we
have

$\displaystyle\|u^{n}\|^{2}_{\infty}$
$\displaystyle\lesssim\|u^{n}\|^{2}_{\dot{H}^{1}}\left(\sqrt{\log(3+c_{1}\left(\sqrt{\frac{N}{\nu}}+\frac{1}{S}\right)(E^{n-1}+1))}\right)^{2}+1$
(3.23)
$\displaystyle\lesssim\frac{2E^{0}}{\nu}\left(1+\log(S)+\log(\frac{1}{\nu})+(\log(N))\right)+1$
$\displaystyle\lesssim_{E^{0}}\nu^{-1}\left(1+\log(S)+\log(\frac{1}{\nu})+\log(N)\right)+1\
.$

Define $m_{0}\coloneqq\nu^{-1}\left(1+\log(S)+|\log(\nu)|+\log(N)\right)$, and
note that $E^{0}\leq\sup_{N}\ E(\Pi_{N}u_{0})\lesssim E_{0}+1$, the inequality
above (3.23) is then rewritten as follows:

$\|u^{n}\|^{2}_{\infty}\lesssim_{E_{0}}m_{0}+1.$

On the other hand by Lemma 3.3,

$\displaystyle\|u^{n+1}\|_{\infty}^{2}$ $\displaystyle\lesssim
1+\|u^{n+1}\|_{\dot{H}^{1}}^{2}\log(3+\|u^{n+1}\|_{\dot{H}^{\frac{3}{2}}})$
(3.24) $\displaystyle\lesssim
1+(1+\frac{1+\|u^{n}\|^{2}_{\infty}}{S})^{2}\|u^{n}\|_{\dot{H}^{1}}^{2}\log(3+\|u^{n+1}\|_{\dot{H}^{\frac{3}{2}}})$
$\displaystyle\lesssim_{E_{0}}1+(1+\frac{m_{0}+1}{S})^{2}\frac{1}{\nu}\log(3+\|u^{n+1}\|_{\dot{H}^{\frac{3}{2}}})$
$\displaystyle\lesssim_{E_{0}}1+(1+\frac{m_{0}}{S})^{2}m_{0}$
$\displaystyle\lesssim_{E_{0}}1+\frac{m_{0}^{3}}{S^{2}}+m_{0}.$

The sufficient condition (3.18) thus becomes

$\left\\{\begin{aligned} &\sqrt{\frac{2S}{\tau}}+\frac{S}{2}-\frac{\nu}{2}\geq
C(E_{0})\left(1+m_{0}+\frac{m_{0}^{3}}{S^{2}}\right),\\\
&m_{0}=\nu^{-1}\left(1+\log(S)+|\log(\nu)|+\log(N)\right)\
.\end{aligned}\right.$

It then suffices for us to choose $S$ such that

$S\gg_{E_{0}}m_{0}+\nu=\nu^{-1}\left(1+\log(S)+|\log(\nu)|+\log(N)\right)+\nu\
,$

where $B\gg_{E_{0}}D$ means there exists a large constant depending only on
$E_{0}$. In fact, for $\nu\gtrsim 1$ , we can take $S\gg_{E_{0}}\nu+\log(N)$;
if $0<\nu\ll 1$, we will choose $S=C_{E_{0}}\cdot\nu^{-1}(|\log\nu|+\log N)$ ,
where $C_{E_{0}}$ is a large constant depending only on $E_{0}$. Therefore it
suffices to choose

$S=C_{E_{0}}\nu^{-1}(\log N+|\log\nu|)+\nu.$ (3.25)

Step 2: We now check the induction base step $n=1$. Clearly we only need to
check

$\sqrt{\frac{2S}{\tau}}+\frac{S}{2}-\frac{\nu}{2}\geq\|\Pi_{N}u_{0}\|^{2}_{\infty}+\frac{1}{2}\|u^{1}\|^{2}_{\infty}.$

By Lemma 3.3,

$\displaystyle\|u^{1}\|_{\dot{H}^{1}}$
$\displaystyle\leq\left(1+\frac{1}{S}+\frac{1}{S}\|\Pi_{N}u_{0}\|^{2}_{\infty}\right)\cdot\|u_{0}\|_{\dot{H}^{1}}$
$\displaystyle\leq\left(1+\frac{1}{S}+\frac{1}{S}\|\Pi_{N}u_{0}\|^{2}_{\infty}\right)\cdot\sqrt{\frac{2E^{0}}{\nu}}\
.$

As a result,

$\displaystyle\|u^{1}\|_{\infty}$ $\displaystyle\lesssim
1+\|u^{1}\|_{\dot{H}^{1}}\sqrt{\log(3+\|u^{1}\|_{\dot{H}^{\frac{3}{2}}})}$
$\displaystyle\lesssim
1+\left(1+\frac{1}{S}+\frac{1}{S}\|\Pi_{N}u_{0}\|^{2}_{\infty}\right)\sqrt{\frac{2E^{0}}{\nu}}\sqrt{\log\left(3+c_{1}\left(\sqrt{\frac{S}{\nu}}+\frac{1}{S}\right)(E_{0}+1)\right)}$
$\displaystyle\lesssim_{E_{0}}1+\left(1+\frac{1}{S}+\frac{1}{S}\|\Pi_{N}u_{0}\|^{2}_{\infty}\right)\cdot\nu^{-\frac{1}{2}}\cdot\sqrt{1+\log(S)+|\log(\nu)|+\log
N}.$

Thus we need to choose $S$ such that

$\displaystyle\sqrt{\frac{2S}{\tau}}+\frac{S}{2}-\frac{\nu}{2}$
$\displaystyle\geq\|\Pi_{N}u_{0}\|^{2}_{\infty}+C_{E_{0}}\cdot\left(1+\frac{1}{S}+\frac{1}{S}\|\Pi_{N}u_{0}\|^{2}_{\infty}\right)^{2}\cdot\nu^{-1}$
$\displaystyle\cdot\left(1+\log(S)+|\log(\nu)|+\log N\right)\ ,$

where $C_{E_{0}}$ is a large constant depending only on $E_{0}$. Note that by
Morrey’s inequality,

$\|\Pi_{N}u_{0}\|_{L^{\infty}(\mathbb{T}^{2})}\lesssim\|\Pi_{N}u_{0}\|_{H^{2}(\mathbb{T}^{2})}\lesssim\|u_{0}\|_{H^{2}(\mathbb{T}^{2})}.$

Then it suffices to take $S$ such that

$S\gg_{E_{0}}\|u_{0}\|^{2}_{H^{2}}+\nu^{-1}(|\log(\nu)|+\log N)+\nu\ .$ (3.26)

This completes the induction and hence proves the theorem.

∎

##  4\. $L^{2}$ error estimate for the scheme (1.4) on the 2D Cahn-Hilliard
equation

In this section, we will like to study the $L^{2}$ error between the semi-
implicit numerical solution and the exact PDE solution in the domain
$\mathbb{T}^{2}$ and eventually prove Theorem 1.7. To start with, we consider
the auxiliary $L^{2}$ error estimate for near solutions.

###  4.1. Auxiliary $L^{2}$ error estimate for near solutions

Consider the following auxiliary system $u^{n}$ and $v^{n}$ in $X_{N}$ for the
first order scheme:

$\left\\{\begin{aligned}
&\frac{u^{n+1}-u^{n}}{\tau}+S\Delta^{2}(u^{n+1}-u^{n})=\frac{e^{-\nu\Delta^{2}\tau}-1}{\tau}u^{n}+\frac{1-e^{-\nu\Delta^{2}\tau}}{\tau\nu\Delta}\Pi_{N}f(u^{n})+\Delta
G^{1}_{n}\\\
&\frac{v^{n+1}-v^{n}}{\tau}+S\Delta^{2}(v^{n+1}-v^{n})=\frac{e^{-\nu\Delta^{2}\tau}-1}{\tau}v^{n}+\frac{1-e^{-\nu\Delta^{2}\tau}}{\tau\nu\Delta}\Pi_{N}f(v^{n})+\Delta
G^{2}_{n}\\\ &u^{0}=u_{0}\ ,\ v^{0}=v_{0}\ .\end{aligned}\right.$ (4.1)

We define that $G_{n}=G^{1}_{n}-G^{2}_{n}$.

######  Proposition 4.1.

For solutions of (4.1), assume for some $N_{1}>0$,

$\sup_{n\geq 0}\|\nabla u^{n}\|_{2}+\sup_{n\geq
0}\|v^{n}\|_{\infty}+\sup_{n\geq 0}\|\nabla v^{n}\|_{2}\leq N_{1}\ .$ (4.2)

Then for any $m\geq 1$,

$\displaystyle\|u^{m}-v^{m}\|^{2}_{2}=\|\varepsilon^{m}\|^{2}_{2}\leq y_{m}$
(4.3) $\displaystyle\leq$
$\displaystyle\exp\left(m\tau\cdot\frac{C(1+N_{1}^{4})}{\nu}\right)\cdot\left(\|Lu_{0}-Lv_{0}\|_{2}^{2}+S\tau\|Mu_{0}-Mv_{0}\|_{2}^{2}+\frac{4}{\nu}\tau\sum_{n=0}^{m-1}\|L^{2}G_{n}\|_{2}^{2}\right).$

where $C>0$ is an absolute constant and $L,M$ are defined from the Fourier
side as below:

$\displaystyle\begin{cases}&\widehat{Lu}(k)=\sqrt{\frac{\nu|k|^{4}\tau}{1-e^{-\nu|k|^{4}\tau}}}\widehat{u}(k)\\\
&\widehat{Mu}(k)=\sqrt{\frac{\nu|k|^{8}\tau}{1-e^{-\nu|k|^{4}\tau}}}\widehat{u}(k).\end{cases}$
(4.4)

The operator $L^{2}=L\circ L$, i.e. $L$ composite with $L$.

######  Lemma 4.2.

We remark here by similar arguments in Lemma 3.2 we have for $g,h\in
C^{\infty}(\mathbb{T}^{2})$ that the following holds:

$\|\Delta g\|_{2}\leq\|Mg\|_{2}\lesssim\|\Delta
g\|_{2}+\sqrt{\tau}\|\Delta^{2}g\|_{2},\quad\|h\|_{2}\leq\|Lh\|_{2}\lesssim\|h\|_{2}+\sqrt{\tau}\|\Delta
h\|_{2}.$ (4.5)

###### Proof of Proposition 4.1.

We assume Lemma 4.2 for the time being. Write $\varepsilon^{n}=u^{n}-v^{n}$.
Then

$\frac{\varepsilon^{n+1}-\varepsilon^{n}}{\tau}+S\Delta^{2}(\varepsilon^{n+1}-\varepsilon^{n})=\frac{e^{-\nu\Delta^{2}\tau}-1}{\tau}\varepsilon^{n}+\frac{1-e^{-\nu\Delta^{2}\tau}}{\tau\nu\Delta}\Pi_{N}(f(u^{n})-f(v^{n}))+\Delta
G_{n}\ .$ (4.6)

Taking $L^{2}$-inner product with
$\frac{\nu\Delta^{2}\tau}{1-e^{-\nu\Delta^{2}\tau}}\varepsilon^{n+1}$ on both
sides of (4.6) and recalling similar computations in previous section, we have
the LHS of (4.6) is

$\displaystyle\mbox{LHS}=$
$\displaystyle\left(\frac{\varepsilon^{n+1}-\varepsilon^{n}}{\tau}+S\Delta^{2}(\varepsilon^{n+1}-\varepsilon^{n}),\
\frac{\nu\Delta^{2}\tau}{1-e^{-\nu\Delta^{2}\tau}}\varepsilon^{n+1}\right)$
(4.7) $\displaystyle=$
$\displaystyle\frac{1}{2\tau}(\|L\varepsilon^{n+1}\|_{2}^{2}-\|L\varepsilon^{n}\|_{2}^{2}+\|L(\varepsilon^{n+1}-\varepsilon^{n})\|_{2}^{2})+\frac{S}{2}(\|M\varepsilon^{n+1}\|_{2}^{2}-\|M\varepsilon^{n}\|_{2}^{2}+\|M(\varepsilon^{n+1}-\varepsilon^{n})\|_{2}^{2}).$

Similarly, the RHS of (4.6) can be given as below

$\displaystyle\mbox{RHS}=$
$\displaystyle\left(\frac{e^{-\nu\Delta^{2}\tau}-1}{\tau}\varepsilon^{n}+\Delta
G_{n}+\frac{1-e^{-\nu\Delta^{2}\tau}}{\tau\nu\Delta}\Pi_{N}(f(u^{n})-f(v^{n})),\
\frac{\nu\Delta^{2}\tau}{1-e^{-\nu\Delta^{2}\tau}}\varepsilon^{n+1}\right)$
(4.8) $\displaystyle\coloneqq$ $\displaystyle I_{1}+I_{2}+I_{3},$

where $(.\ ,\ .)$ denotes the $L^{2}$ inner product. We then estimate
different parts $I_{1},I_{2}$ and $I_{3}$ as follows.

$\displaystyle I_{1}=$
$\displaystyle\left(\frac{e^{-\nu\Delta^{2}\tau}-1}{\tau}\varepsilon^{n},\frac{\nu\Delta^{2}\tau}{1-e^{-\nu\Delta^{2}\tau}}\varepsilon^{n+1}\right)$
(4.9) $\displaystyle=$
$\displaystyle-\nu(\Delta\varepsilon^{n},\Delta\varepsilon^{n+1})$
$\displaystyle=$
$\displaystyle-\nu(\Delta\varepsilon^{n+1},\Delta\varepsilon^{n+1})+\nu(\Delta\varepsilon^{n+1},\Delta\varepsilon^{n+1}-\Delta\varepsilon^{n})$
$\displaystyle=$
$\displaystyle-\frac{\nu}{2}(\|\Delta\varepsilon^{n+1}\|_{2}^{2}+\|\Delta\varepsilon^{n}\|_{2}^{2})+\frac{\nu}{2}\|\Delta\varepsilon^{n+1}-\Delta\varepsilon^{n}\|_{2}^{2}.$

Similarly, by Hölder’s inequality we obtain that

$I_{2}=\left(\Delta
G_{n},\frac{\nu\Delta^{2}\tau}{1-e^{-\nu\Delta^{2}\tau}}\varepsilon^{n+1}\right)\leq\|\Delta\varepsilon^{n+1}\|_{2}\|L^{2}G_{n}\|_{2}\leq\frac{2}{\nu}\|L^{2}G_{n}\|^{2}_{2}+\frac{\nu}{8}\|\Delta\varepsilon^{n+1}\|_{2}^{2},$
(4.10)

where $L^{2}(G_{n})=L(LG_{n})$. On the other hand, note that $\Pi_{N}$ is a
self-adjoint operator $\left(\Pi_{N}f,\ g\right)=\left(f,\ \Pi_{N}g\right)$,
since it is just an $N$-th Fourier mode truncation. Therefore we have

$\displaystyle
I_{3}=\left(\frac{1-e^{-\nu\Delta^{2}\tau}}{\tau\nu\Delta}\Pi_{N}(f(u^{n})-f(v^{n})),\frac{\nu\Delta^{2}\tau}{1-e^{-\nu\Delta^{2}\tau}}\varepsilon^{n+1}\right)=\left(f(u^{n})-f(v^{n}),\Delta\Pi_{N}\varepsilon^{n+1}\right).$
(4.11)

Then by the fundamental theorem of calculus, we have

$\displaystyle f(u^{n})-f(v^{n})$
$\displaystyle=\int_{0}^{1}f^{\prime}(v^{n}+s\varepsilon^{n})ds\cdot\varepsilon^{n}$
(4.12)
$\displaystyle=(a_{1}+a_{2}(v^{n})^{2})\varepsilon^{n}+a_{3}v^{n}(\varepsilon^{n})^{2}+a_{4}(\varepsilon^{n})^{3}\
,$

where $a_{i}$ are exact constants can be computed. Note that we will denote
$C$ to be an absolute constant whose value may vary in different lines. To
proceed further we indeed have by (4.2) that

$\displaystyle|\left((a_{1}+a_{2}(v^{n})^{2})\varepsilon^{n}\ ,\
\Delta\Pi_{N}\varepsilon^{n+1}\right)|$ $\displaystyle\leq
C(1+\|v^{n}\|_{\infty}^{2})\|\Delta\varepsilon^{n+1}\|_{2}\|\varepsilon^{n}\|_{2}$
(4.13)
$\displaystyle\leq\frac{C(1+N_{1}^{4})}{\nu}\|\varepsilon^{n}\|_{2}^{2}+\frac{\nu}{8}\|\Delta\varepsilon^{n+1}\|_{2}^{2}.$

Moreover, the other two terms can be estimated similarly. In particular,
recalling (2.2) that for $f$ defined in $\mathbb{T}^{2}$ with zero mean one
has $\|f\|_{4}^{2}\lesssim\|f\|_{2}\|\nabla f\|_{2}$. Therefore we have

$\displaystyle|\left(a_{3}v^{n}(\varepsilon^{n})^{2}\
,\Delta\Pi_{N}\varepsilon^{n+1}\right)|$ $\displaystyle\leq
C\|v^{n}\|_{\infty}\|\varepsilon^{n}\|_{4}^{2}\|\Delta\varepsilon^{n+1}\|_{2}\leq\frac{CN_{1}^{2}}{\nu}\|\varepsilon^{n}\|_{4}^{4}+\frac{\nu}{8}\|\Delta\varepsilon^{n+1}\|_{2}^{2}$
(4.14) $\displaystyle\leq
C\frac{N_{1}^{2}}{\nu}\|\varepsilon^{n}\|_{2}^{2}\|\nabla\varepsilon^{n}\|_{2}^{2}+\frac{\nu}{8}\|\Delta\varepsilon^{n+1}\|_{2}^{2}\leq\frac{CN_{1}^{4}}{\nu}\|\varepsilon^{n}\|_{2}^{2}+\frac{\nu}{8}\|\Delta\varepsilon^{n+1}\|_{2}^{2}.$

Similarly, recalling (2.3) that we also have

$\displaystyle|\left(a_{4}(\varepsilon^{n})^{3}\
,\Delta\Pi_{N}\varepsilon^{n+1}\right)|$ $\displaystyle\leq
C\|\varepsilon^{n}\|_{6}^{3}\|\Delta\varepsilon^{n+1}\|_{2}\leq\frac{C}{\nu}\|\varepsilon^{n}\|_{6}^{6}+\frac{\nu}{8}\|\Delta\varepsilon^{n+1}\|_{2}^{2}$
(4.15)
$\displaystyle\leq\frac{C}{\nu}\|\varepsilon^{n}\|_{2}^{2}\|\nabla\varepsilon^{n}\|_{2}^{4}+\frac{\nu}{8}\|\Delta\varepsilon^{n+1}\|_{2}^{2}$
$\displaystyle\leq\frac{CN_{1}^{4}}{\nu}\|\varepsilon^{n}\|_{2}^{2}+\frac{\nu}{8}\|\Delta\varepsilon^{n+1}\|_{2}^{2}.$

One can then conclude from (4.12), (4.13), (4.14) and (4.15) that

$|I_{3}|=|\left(f(u^{n})-f(v^{n}),\Delta\Pi_{N}\varepsilon^{n+1}\right)|\leq\frac{C(1+N_{1}^{4})}{\nu}\|\varepsilon^{n}\|_{2}^{2}+\frac{3\nu}{8}\|\Delta\varepsilon^{n+1}\|_{2}^{2}.$
(4.16)

Collecting the estimates (4.8), (4.9), (4.10) and (4.16), we get

$\displaystyle\mbox{RHS}\leq\frac{\nu}{2}\|\Delta\varepsilon^{n+1}-\Delta\varepsilon^{n}\|_{2}^{2}+\frac{2}{\nu}\|L^{2}G_{n}\|_{2}^{2}+\frac{C(1+N_{1}^{4})}{\nu}\|\varepsilon^{n}\|_{2}^{2}.$
(4.17)

Taking the same choice of $S$ as in Theorem 1.5, it then follows by (4.7) and
(4.17) that

$\displaystyle\frac{\|L\varepsilon^{n+1}\|^{2}_{2}-\|L\varepsilon^{n}\|_{2}^{2}}{\tau}+S(\|M\varepsilon^{n+1}\|^{2}_{2}-\|M\varepsilon^{n}\|_{2}^{2})$
(4.18) $\displaystyle\leq$
$\displaystyle\frac{4}{\nu}\|L^{2}G_{n}\|_{2}^{2}+\frac{C(1+N_{1}^{4})}{\nu}\|\varepsilon^{n}\|_{2}^{2},$

where $C$ is an absolute positive constant that can be computed exactly. Then
by defining

$\displaystyle y_{n}=\|L\varepsilon^{n}\|_{2}^{2}+\tau
S\|M\varepsilon^{n}\|_{2}^{2},$ (4.19)
$\displaystyle\alpha=\frac{C(1+N_{1}^{4})}{\nu},$
$\displaystyle\beta_{n}=\frac{4}{\nu}\|L^{2}G_{n}\|_{2}^{2},$

(4.18) thus can be formulated into

$\frac{y_{n+1}-y_{n}}{\tau}\leq\alpha y_{n}+\beta_{n}\ .$

Applying discrete Grönwall’s inequality Lemma 2.6, we arrive at

$\displaystyle\|u^{m}-v^{m}\|^{2}_{2}=\|\varepsilon^{m}\|^{2}_{2}\leq y_{m}$
(4.20) $\displaystyle\leq$
$\displaystyle\exp\left(m\tau\cdot\frac{C(1+N_{1}^{4})}{\nu}\right)\cdot\left(\|Lu_{0}-Lv_{0}\|_{2}^{2}+S\tau\|Mu_{0}-Mv_{0}\|_{2}^{2}+\frac{4}{\nu}\tau\sum_{n=0}^{m-1}\|L^{2}G_{n}\|_{2}^{2}\right).$

∎

######  Remark 4.3.

The proof of Proposition 4.1 only requires the uniform boundedness of
$\|\nabla u^{n}\|_{2}+\|\nabla v^{n}\|_{2}$ rather than $\|\Delta
u^{n}\|_{2}+\|\Delta v^{n}\|_{2}$. It somehow makes sense by considering the
equation (1.4) where the nonlinear part requires lower regularity assumption
and the linear part can be controlled by its sign.

###### Proof of Lemma 4.2.

First of all the same arguments from Lemma 3.2 give

$\displaystyle\|\Delta g\|_{2}\leq\|Mg\|_{2},\quad\|h\|_{2}\leq\|Lh\|_{2}.$

To see the rest, we consider the Fourier series of $Lh$ again:

$\displaystyle\|Lh\|_{2}^{2}=$ $\displaystyle\frac{1}{(2\pi)^{2}}\sum_{|k|\geq
1}\frac{\nu|k|^{4}\tau}{1-e^{-\nu|k|^{4}\tau}}|\widehat{h}(k)|^{2}$ (4.21)
$\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{2}}(\sum_{\nu|k|^{4}\tau\leq\log
2}+\sum_{\nu|k|^{4}\tau\geq\log
2}\frac{\nu|k|^{4}\tau}{1-e^{-\nu|k|^{4}\tau}}|\widehat{h}(k)|^{2})$
$\displaystyle\lesssim$ $\displaystyle\sum_{\nu|k|^{4}\tau\leq\log
2}\frac{\nu|k|^{4}\tau}{1-e^{-\nu|k|^{4}\tau}}|\widehat{h}(k)|^{2}+\tau\sum_{\nu|k|^{4}\tau\geq\log
2}|k|^{4}|\widehat{h}(k)|^{2}.$

The high frequency part can be controlled by $\|\Delta h\|_{2}$ by using the
standard Parseval’s identity. For the lower frequency part, we then consider
the function $\phi(x)=\frac{x}{1-e^{-x}}$. By direct computation, we see that
$\phi^{\prime}(x)>0$ for $x>0$ which shows that $\phi(x)$ is monotone
increasing. Therefore we have

$\displaystyle\sum_{\nu|k|^{4}\tau\leq\log
2}\frac{\nu|k|^{4}\tau}{1-e^{-\nu|k|^{4}\tau}}|\widehat{h}(k)|^{2}\lesssim\sum_{|k|\geq
1}|\widehat{h}(k)|^{2}.$ (4.22)

Collecting the estimates (4.21) and (4.22) we then get the desired result. The
same reason gives

$\displaystyle\|Mg\|_{2}\lesssim\|\Delta
g\|_{2}+\sqrt{\tau}\|\Delta^{2}g\|_{2}.$

∎

###  4.2. $L^{2}$ error estimate of the 2D Cahn-Hilliard equation

In this section, to simplify the notation, we will write $x\lesssim y$ if
$x\leq C(\nu,u_{0})\ y$ for a constant $C$ depending on $\nu$ and $u_{0}$. We
consider the system

$\left\\{\begin{aligned}
&\frac{u^{n+1}-u^{n}}{\tau}+S\Delta^{2}(u^{n+1}-u^{n})=\frac{e^{-\nu\Delta^{2}\tau}-1}{\tau}u^{n}+\frac{1-e^{-\nu\Delta^{2}\tau}}{\tau\nu\Delta}\Pi_{N}f(u^{n})\\\
&\partial_{t}u=-\nu\Delta^{2}u+\Delta(f(u))\\\ &u^{0}=\Pi_{N}u_{0}\ ,\
u(0)=u_{0}\ .\end{aligned}\right.$ (4.23)

In order to prove Theorem 1.7, it is clear that we shall estimate $G_{n}$
introduced in (4.1) from the previous proposition. Note that for a one-
variable function $h(t)$, one has the formula:

$\left\\{\begin{aligned}
&\frac{1}{\tau}\int_{t_{n}}^{t_{n+1}}h(t)=h(t_{n})+\frac{1}{\tau}\int_{t_{n}}^{t_{n+1}}h^{\prime}(t)\cdot(t_{n+1}-t)\
dt\\\
&\frac{1}{\tau}\int_{t_{n}}^{t_{n+1}}h(t)=h(t_{n+1})+\frac{1}{\tau}\int_{t_{n}}^{t_{n+1}}h^{\prime}(t)\cdot(t_{n}-t)\
dt\ .\end{aligned}\right.$ (4.24)

Using the formula (4.24) above and integrating the Cahn-Hilliard equation (CH)
on the time interval $[t_{n},t_{n+1}]$, we get

$\displaystyle\frac{u(t_{n+1})-u(t_{n})}{\tau}=$
$\displaystyle-\nu\Delta^{2}u(t_{n})-S\Delta^{2}\left(u(t_{n+1})-u(t_{n})\right)$
(4.25)
$\displaystyle+\Delta\Pi_{N}f(u)(t_{n})+\Delta\Pi_{>N}f(u)(t_{n})+\Delta
Q_{n}\ ,$

where $\Pi_{>N}=id-\Pi_{N}$, the large mode truncation operator, and

$\displaystyle Q_{n}=-\frac{\nu}{\tau}\int_{t_{n}}^{t_{n+1}}\partial_{t}\Delta
u\cdot(t_{n+1}-t)\
dt+\frac{1}{\tau}\int_{t_{n}}^{t_{n+1}}\partial_{t}(f(u))(t_{n+1}-t)\
dt+S\int_{t_{n}}^{t_{n+1}}\partial_{t}\Delta u\ dt\ .$ (4.26)

We then rewrite (4.25) into the following form:

$\displaystyle\frac{u(t_{n+1})-u(t_{n})}{\tau}+$ $\displaystyle
S\Delta^{2}\left(u(t_{n+1})-u(t_{n})\right)=\frac{e^{-\nu\Delta^{2}\tau}-1}{\tau}u(t_{n})+\frac{1-e^{-\nu\Delta^{2}\tau}}{\tau\nu\Delta}\Pi_{N}f(u(t_{n}))$
(4.27)
$\displaystyle+\left(\frac{1-e^{-\nu\Delta^{2}\tau}}{\tau}-\nu\Delta^{2}\right)u(t_{n})+\left(\Delta+\frac{e^{-\nu\Delta^{2}\tau}-1}{\nu\tau\Delta}\right)\Pi_{N}f(u(t_{n}))$
$\displaystyle+\Delta\Pi_{>N}f(u)(t_{n})+\Delta Q_{n}\ ,$

For the sake of simplicity, we denote $H_{1},H_{2},H_{3}$ to be the following
quantities:

$\displaystyle\begin{cases}&H_{1}\coloneqq\left(\frac{1-e^{-\nu\Delta^{2}\tau}}{\tau\Delta}-\nu\Delta\right)u(t_{n})=\frac{1-\nu\Delta^{2}\tau-e^{-\nu\Delta^{2}\tau}}{\tau\Delta}u(t_{n}),\\\
&H_{2}\coloneqq\left(1+\frac{e^{-\nu\Delta^{2}\tau}-1}{\nu\tau\Delta^{2}}\right)\Pi_{N}f(u)(t_{n}),\\\
&H_{3}\coloneqq\Pi_{>N}f(u)(t_{n}).\end{cases}$ (4.28)

In order to prove Theorem 1.7, we now estimate $\|L^{2}G_{n}\|_{2}$ by
estimating $\|L^{2}Q_{n}\|_{2}$ and $\|L^{2}H_{i}\|_{2}$, $i=1,2,3$. Recall
that $L^{2}u=\frac{\nu\Delta^{2}\tau}{1-e^{-\nu\Delta^{2}\tau}}u$.

###### Proof of Theorem 1.7.

By Theorem 1.5 and Lemma 2.7, $\sup_{n\geq 0}\|u^{n}\|_{H^{1}}\lesssim 1$ and
$\sup_{t>0}\|u(t)\|_{H^{s}}\lesssim 1$ for $s\geq 10$ given in the statement.
Thus the assumptions of Proposition 4.1 (auxiliary $L^{2}$ error estimate
proposition) are satisfied for (4.23) and (4.27). Recall that

$\displaystyle Q_{n}=-\frac{\nu}{\tau}\int_{t_{n}}^{t_{n+1}}\partial_{t}\Delta
u\cdot(t_{n+1}-t)\
dt+\frac{1}{\tau}\int_{t_{n}}^{t_{n+1}}\partial_{t}(f(u))(t_{n+1}-t)\
dt+S\int_{t_{n}}^{t_{n+1}}\partial_{t}\Delta u\ dt\ .$

Recalling (4.5) that we can control
$\|L^{2}Q_{n}\|_{2}\lesssim\|Q_{n}\|_{2}+\sqrt{\tau}\|\Delta
Q_{n}\|_{2}+\tau\|\Delta^{2}Q_{n}\|_{2}$ as follows:

$\displaystyle\|Q_{n}\|_{2}$
$\displaystyle\lesssim\int_{t_{n}}^{t_{n+1}}\|\partial_{t}\Delta u\|_{2}\
dt+\int_{t_{n}}^{t_{n+1}}\|\partial_{t}(f(u))\|_{2}\
dt+S\int_{t_{n}}^{t_{n+1}}\|\partial_{t}\Delta u\|_{2}\ dt$ (4.29)
$\displaystyle\lesssim(1+S)\int_{t_{n}}^{t_{n+1}}\|\partial_{t}\Delta u\|_{2}\
dt+\int_{t_{n}}^{t_{n+1}}\|\partial_{t}u\|_{2}\
dt\cdot\|f^{\prime}(u)\|_{L^{\infty}_{t}L^{\infty}_{x}}.$

Note that by Lemma 2.7,

$\|\partial_{t}u\|_{2}+\|\partial_{t}\Delta
u\|_{2}+\|f^{\prime}(u)\|_{\infty}\lesssim 1\ .$

Therefore we can estimate (4.29) as following

$\displaystyle\int_{t_{n}}^{t_{n+1}}\|\partial_{t}\Delta u\|_{2}\
dt\lesssim\left(\int_{t_{n}}^{t_{n+1}}\|\partial_{t}\Delta u\|_{2}^{2}\
dt\right)^{\frac{1}{2}}\cdot\sqrt{\tau}\ .$

Therefore when $t_{m}\geq 1$, it is not hard to see that

$\displaystyle\sum_{n=0}^{m-1}\|Q_{n}\|_{2}^{2}$
$\displaystyle\lesssim\sum_{n=0}^{m-1}\left(\tau(1+S)^{2}\int_{t_{n}}^{t_{n+1}}\|\partial_{t}\Delta
u\|_{2}^{2}\ dt+\tau\int_{t_{n}}^{t_{n+1}}\|\partial_{t}u\|_{2}^{2}\
dt\right)$ (4.30)
$\displaystyle\lesssim\tau(1+S)^{2}\int_{0}^{t_{m}}\|\partial_{t}\Delta
u\|_{2}^{2}\ dt+\tau\int_{0}^{t_{m}}\|\partial_{t}u\|_{2}^{2}\ dt$
$\displaystyle\lesssim\tau(1+S)^{2}(1+t_{m})+\tau$
$\displaystyle\lesssim(1+S)^{2}\tau\cdot(1+t_{m})\ .$

For the rest terms $\sqrt{\tau}\|\Delta Q_{n}\|_{2}$ and
$\tau\|\Delta^{2}Q_{n}\|_{2}$ we will obtain very similar results as long as
$s\geq 10$ and thereby we can conclude that:

$\displaystyle\sum_{n=0}^{m-1}\|L^{2}Q_{n}\|_{2}^{2}\lesssim(1+S)^{2}\tau\cdot(1+t_{m}).$
(4.31)

We now focus on the $\|L^{2}H_{i}\|_{2}$ terms. Firstly it is easy to see from
the Fourier side that

$\displaystyle\|L^{2}H_{1}\|_{2}^{2}=$
$\displaystyle\frac{1}{(2\pi)^{2}}\sum_{|k|\geq
1}\frac{\nu^{2}|k|^{4}(1-\nu|k|^{4}\tau-e^{-\nu|k|^{4}\tau})^{2}}{(1-e^{-\nu|k|^{4}\tau})^{2}}|\widehat{u}(k)|^{2}(t_{n})$
(4.32) $\displaystyle\lesssim$ $\displaystyle\sum_{\nu|k|^{4}\tau\leq
c}+\sum_{\nu|k|^{4}\tau\geq
c}\frac{|k|^{4}(1-\nu|k|^{4}\tau-e^{-\nu|k|^{4}\tau})^{2}}{(1-e^{-\nu|k|^{4}\tau})^{2}}|\widehat{u}(k)|^{2}(t_{n})$
$\displaystyle\coloneqq$ $\displaystyle J_{1}+J_{2}.$

For the part where $\nu|k|^{4}\tau\geq c$, we indeed have

$\displaystyle J_{2}\coloneqq$ $\displaystyle\sum_{\nu|k|^{4}\tau\geq
c}\frac{|k|^{4}(1-\nu|k|^{4}\tau-e^{-\nu|k|^{4}\tau})^{2}}{(1-e^{-\nu|k|^{4}\tau})^{2}}|\widehat{u}(k)|^{2}(t_{n})$
(4.33) $\displaystyle\lesssim$ $\displaystyle\sum_{\nu|k|^{4}\tau\geq
c}|k|^{12}\tau^{2}|\widehat{u}(k)|^{2}(t_{n})$ $\displaystyle\lesssim$
$\displaystyle\tau^{2}\|\Delta^{3}u\|^{2}_{2}(t_{n}).$

For the lower frequency part, we get

$\displaystyle J_{1}\coloneqq$ $\displaystyle\sum_{\nu|k|^{4}\tau\leq
c}\frac{|k|^{4}(1-\nu|k|^{4}\tau-e^{-\nu|k|^{4}\tau})^{2}}{(1-e^{-\nu|k|^{4}\tau})^{2}}|\widehat{u}(k)|^{2}(t_{n})$
$\displaystyle\lesssim$
$\displaystyle\sum_{1\leq|k|\lesssim(\frac{1}{\tau})^{\frac{1}{4}}}\frac{|k|^{8}(1-\nu|k|^{4}\tau-e^{-\nu|k|^{4}\tau})^{2}}{(1-e^{-\nu|k|^{4}\tau})^{2}}|\widehat{u}(k)|^{2}(t_{n}).$

Recall that $\phi(x)=\frac{x}{1-e^{-x}}$ is monotone increasing for $x>0$ as
in Lemma 4.2. Moreover, we have (by a simple observation) that

$\displaystyle 0<\nu|k|^{4}\tau+e^{-\nu|k|^{4}\tau}-1<\nu^{2}|k|^{8}\tau^{2}.$
(4.34)

Therefore we have

$\displaystyle J_{1}\lesssim$
$\displaystyle\frac{1}{\tau^{2}}\sum_{1\leq|k|\lesssim(\frac{1}{\tau})^{\frac{1}{4}}}(\nu^{2}|k|^{8}\tau^{2})^{2}|\widehat{u}(k)|^{2}(t_{n})\lesssim\tau^{2}\sum_{|k|\geq
1}|k|^{16}|\widehat{u}(k)|^{2}(t_{n})\lesssim\tau^{2}\|\Delta^{4}u\|_{2}^{2}(t_{n}).$
(4.35)

Combing (4.33)-(4.35), we get

$\displaystyle\sum_{n=0}^{m-1}\|L^{2}H_{1}\|_{2}^{2}\lesssim\sum_{n=0}^{m-1}\tau^{2}(\|\Delta^{3}u\|_{2}^{2}(t_{n})+\|\Delta^{4}u\|_{2}^{2}(t_{n}))\lesssim
m\tau^{2}\lesssim(t_{m}+1)\cdot\tau,$ (4.36)

where the uniform Sobolev bound lemma (Lemma 2.7) is applied.

Secondly we estimate $\|L^{2}H_{2}\|_{2}^{2}$. Note that again by considering
the Fourier expansion we arrive at:

$\displaystyle\|L^{2}H_{2}\|_{2}^{2}=$
$\displaystyle\frac{1}{(2\pi)^{2}}\sum_{1\leq|k|\leq
N}\frac{(\nu\tau|k|^{4}+e^{-\nu|k|^{4}\tau}-1)^{2}}{(1-e^{-\nu\tau|k|^{4}})^{2}}|\widehat{f(u)}(k)|^{2}(t_{n})$
(4.37) $\displaystyle\lesssim$ $\displaystyle\sum_{\nu|k|^{4}\tau\leq
c}+\sum_{\nu|k|^{4}\tau\geq
c}\frac{(\nu\tau|k|^{4}+e^{-\nu|k|^{4}\tau}-1)^{2}}{(1-e^{-\nu\tau|k|^{4}})^{2}}|\widehat{f(u)}(k)|^{2}(t_{n})$
$\displaystyle\coloneqq$ $\displaystyle J_{3}+J_{4}.$

Similar to (4.32) above, we have

$\displaystyle J_{4}\coloneqq$ $\displaystyle\sum_{\nu|k|^{4}\tau\geq
c}\frac{(\nu\tau|k|^{4}+e^{-\nu|k|^{4}\tau}-1)^{2}}{(1-e^{-\nu\tau|k|^{4}})^{2}}|\widehat{f(u)}(k)|^{2}(t_{n})$
(4.38) $\displaystyle\lesssim$ $\displaystyle\sum_{|k|\geq
1}|k|^{8}\tau^{2}|\widehat{f(u)}(k)|^{2}(t_{n})$ $\displaystyle\lesssim$
$\displaystyle\tau^{2}\|\Delta^{2}f(u)\|_{2}^{2}(t_{n})\lesssim\tau^{2}.$

Note that here by Lemma 2.7 we have $\sup_{t\geq
0}\|u\|_{H^{s}}(t)\lesssim_{s}1$, which leads to $\sup_{n\geq
0}\|f(u)\|_{H^{s}}(t_{n})\lesssim_{s}1$. For $J_{3}$, we have

$\displaystyle J_{3}\coloneqq$ $\displaystyle\sum_{\nu|k|^{4}\tau\leq
c}\frac{(\nu\tau|k|^{4}+e^{-\nu|k|^{4}\tau}-1)^{2}}{(1-e^{-\nu\tau|k|^{4}})^{2}}|\widehat{f(u)}(k)|^{2}(t_{n})$
(4.39) $\displaystyle\lesssim$
$\displaystyle\sum_{1\leq|k|\lesssim(\frac{1}{\tau})^{\frac{1}{4}}}\frac{|k|^{8}(\nu|k|^{4}\tau+e^{-\nu|k|^{4}\tau}-1)^{2}}{(1-e^{-\nu|k|^{4}\tau})^{2}}|\widehat{f(u)}(k)|^{2}(t_{n})$
$\displaystyle\lesssim$
$\displaystyle\frac{1}{\tau^{2}}\sum_{1\leq|k|\lesssim(\frac{1}{\tau})^{\frac{1}{4}}}\tau^{4}|k|^{16}|\widehat{f(u)}(k)|^{2}(t_{n})\lesssim\tau^{2}\|\Delta^{4}f(u)\|_{2}^{2}(t_{n}),$

where again we use the monotonicity of $\phi(x)=\frac{x}{1-e^{-x}}$ for $x>0$
and (4.34). We can then derive that

$\displaystyle\sum_{n=0}^{m-1}\|L^{2}H_{2}\|_{2}^{2}\lesssim
m\tau^{2}\lesssim(t_{m}+1)\cdot\tau.$ (4.40)

Lastly let us focus on $\|L^{2}H_{3}\|_{2}^{2}$. Note that
$\|L^{2}H_{3}\|_{2}^{2}\lesssim\|H_{3}\|_{2}^{2}+\tau\|\Delta
H_{3}\|_{2}^{2}+\tau^{2}\|\Delta^{2}H_{3}\|_{2}^{2}$, we then estimate
$\|H_{3}\|_{2}^{2}$ below.

$\displaystyle\|H_{3}\|_{2}^{2}=\|\Pi_{>N}f(u)\|_{2}^{2}(t_{n})$
$\displaystyle=\sum_{|k|>N}\left|\widehat{f(u)}(k)\right|^{2}(t_{n})$ (4.41)
$\displaystyle\leq\sum_{|k|>N}|k|^{2s}\left|\widehat{f(u)}(k)\right|^{2}(t_{n})\cdot{|k|^{-2s}}$
$\displaystyle\lesssim N^{-2s}\cdot\|f(u)\|_{H^{s}}^{2}(t_{n})$
$\displaystyle\lesssim N^{-2s}\ .$

Similarly, we get $\tau\|\Delta H_{3}\|_{2}^{2}\lesssim\tau N^{-2s}$ and
$\tau^{2}\|\Delta^{2}H_{3}\|_{2}^{2}\lesssim\tau^{2}N^{-2s}$. As a result,
assuming $\tau\lesssim 1$ or the time step is bounded (which is very natural)
we get

$\sum_{n=0}^{m-1}\|\Pi_{>N}f(u(t_{n}))\|_{2}^{2}\lesssim_{s}m\cdot(1+\tau^{2})N^{-2s}\lesssim\frac{t_{m}N^{-2s}}{\tau}\
.$ (4.42)

Therefore, combining (4.31), (4.36), (4.40) and (4.42) we can conclude that

$\tau\sum_{n=0}^{m-1}\left(\|L^{2}Q_{n}\|_{2}^{2}+\sum_{i=1,2,3}\|L^{2}H_{i}\|_{2}^{2}\right)\lesssim_{s}(1+t_{m})(\tau^{2}+N^{-2s})(1+S)^{2}\
.$ (4.43)

In order to apply Proposition 4.1, it remains to control
$\|Lu^{0}-Lu(0))\|_{2}$ and $\|Mu^{0}-Mu(0)\|_{2}$. Note that by (4.5) in
Lemma 4.2 and similar arguments in (4.41) we have

$\displaystyle\|Lu^{0}-Lu(0)\|_{2}^{2}=\|L\Pi_{N}u_{0}-Lu_{0}\|_{2}^{2}\lesssim\|\Pi_{N}u_{0}-u_{0}\|_{2}^{2}+\tau\|\Delta(\Pi_{N}u_{0}-u_{0})\|_{2}^{2}\lesssim
N^{-2s}+\tau N^{-2(s-2)}.$ (4.44)

Similarly, we have

$\|Mu^{0}-Mu(0)\|_{2}^{2}=\|M\Pi_{N}u_{0}-Mu_{0}\|_{2}^{2}\lesssim
N^{-2(s-2)}+\tau N^{-2(s-4)}\ .$ (4.45)

We now apply the auxiliary solutions estimate in Proposition 4.1. Noting the
estimates (4.43), (4.44) and (4.45) we can get

$\|u^{m}-u(t_{m})\|_{2}^{2}\lesssim_{s}(1+S)^{2}e^{Ct_{m}}\left(N^{-2s}+\tau^{2}\cdot
N^{-2(s-4)}+\tau^{2}\right)\ .$

Thus

$\|u^{m}-u(t_{m})\|_{2}\leq(1+S)\cdot C_{2}\cdot
e^{C_{1}t_{m}}\left(N^{-s}+\tau+\tau\cdot N^{-s+4}\right)\ ,$ (4.46)

where $C_{1}>0$ is a constant depending on $\nu$ and $u_{0}$; $C_{2}>0$ is a
constant depending on $s\ ,\nu$ and $u_{0}$. This completes the proof of
$L^{2}$ error estimate. ∎

##  5\. Numerical experiments

In this section, we present several numerical results. To begin with we
present numerical evidence that show the necessity of the stabilizers.

###  5.1. A benchmark computation with different $S$-values

In the first several experiments we vary the choice of the stabilizers, namely
put $S=0,0.01,0.1$. Moreover for the rest of the parameters, we choose
$\nu=0.01$, $\tau=0.1,~{}N_{x}=N_{y}=256$ and the initial data $u_{0}$ is
given basically “supported” in seven circles as below:

$\displaystyle
u_{0}(x,y)=-1+\sum_{i=1}^{7}f_{0}\left(\sqrt{(x-x_{i})^{2}+(y-y_{i})^{2}}-r_{i}\right),$
(5.1)

where

$\displaystyle f_{0}(s)=\begin{cases}&2e^{-\frac{\nu}{s^{2}}},\quad\mbox{if}\
s<0;\\\ &0,\quad\mbox{otherwise}.\end{cases}$

The centers and radii of the chosen circles are given in the Table 1 below.

$x_{i}$ | $-\frac{\pi}{2}$ | $-\frac{3\pi}{4}$ | $-\frac{\pi}{2}$ | 0 | $\frac{\pi}{2}$ | 0 | $\frac{\pi}{2}$
---|---|---|---|---|---|---|---
$y_{i}$ | $-\frac{\pi}{2}$ | $-\frac{\pi}{4}$ | $\frac{\pi}{4}$ | $-\frac{3\pi}{4}$ | $-\frac{3\pi}{4}$ | 0 | $\frac{\pi}{2}$
$r_{i}$ | $\frac{\pi}{5}$ | $\frac{2\pi}{15}$ | $\frac{2\pi}{15}$ | $\frac{\pi}{10}$ | $\frac{\pi}{10}$ | $\frac{\pi}{4}$ | $\frac{\pi}{4}$
Table 1. Choice of centers and radii

The cases $S=0$ (left) and $S=0.01$ (right) can be found in Fig 1 and the case
$S=0.1$ can be found in the left of Fig 3. As you can tell, the choice of
$S=0.1$ can guarantee the energy dissipation already. Recall the discussion in
Theorem 1.5 and Remark 1.6, we see that the optimal size of $S$ could be much
smaller than $\nu^{-1}\log\nu$ that we will not dig further in this direction.

Figure 1. Dynamics of 2D Cahn-Hilliard equation using EI scheme (1.4) where
$\nu=0.01$, $\tau=0.1,~{}N_{x}=N_{y}=256$ and the initial data $u_{0}$ is
given in (5.1). $S=0$ (Left), $S=0.01$ (right).

Figure 2. Dynamics of 2D Cahn-Hilliard equation using EI scheme (1.4) where
$\nu=0.01$, $S=0.1$, $\tau=0.1,~{}N_{x}=N_{y}=256$ and the initial data
$u_{0}$ is given in (5.1).

Figure 3. Dynamics of 2D Cahn-Hilliard equation using EI scheme (1.4) where
$\nu=0.01$, $S=0.1$, $\tau=0.1,~{}N_{x}=N_{y}=256$ and the initial data
$u_{0}=0.5\sin(x)\sin(y)$ .

###  5.2. Convergence test

In this subsection we consider a benchmark computation test with initial data
being $u_{0}=0.5*\sin(x)\sin(y)$. We take $\nu=1$, $S=0.1$ and
$N_{x}=N_{y}=256$. Then we consider the exact solution
$u_{e}=0.5*e^{-t}\sin(x)\sin(y)$ corresponding to certain forcing term that
can be computed explicitly. With these settings we perform our numerical
experiments with various time steps $\tau=\frac{0.01}{2^{k}}$ with
$k=0,1,...,6$. The relative $L^{2}$-errors and $L^{\infty}$-errors at time
$T=0.5$ are presented below in Table 2. As usual, the experimental order of
convergence is computed by comparing the errors of two consecutive
refinements. Indeed the rate of convergence indicates the order of the error
is $O(\tau)$.

$\tau=0.01$ |  $L^{2}$-error | Rate |  $L^{\infty}$-error | Rate
---|---|---|---|---
$\tau$ | 9.099e-05 | - | 1.156e-03 | -
$\tau/2$ | 4.523e-05 | 2.012 | 5.745e-04 | 2.012
$\tau/4$ | 2.255e-05 | 2.006 | 2.864e-04 | 2.006
$\tau/8$ | 1.126e-05 | 2.003 | 1.430e-04 | 2.003
$\tau/16$ | 5.624e-06 | 2.002 | 7.143e-05 | 2.002
$\tau/32$ | 2.811e-06 | 2.001 | 3.570e-05 | 2.001
Table 2. Errors and orders of convergence

###  5.3. More dynamics and energy evolution

In this section, we present more dynamics of our scheme 1.4 with patterns
indicating the energy dissipation. In Fig 3 we present the dynamics of 2D
Cahn-Hilliard equation using EI scheme (1.4) where $\nu=0.01$, S=0.1,
$\tau=0.1,~{}N_{x}=N_{y}=256$ and the initial data $u_{0}=0.5\sin(x)\sin(y)$.
In Fig 5 we present the dynamics choosing $\nu=0.01$, $S=0.1$,
$\tau=0.001,~{}N_{x}=N_{y}=256$ and the initial data
$u_{0}=0.05\sin(x)\sin(y)$. In Fig 5 we present the dynamics where $\nu=0.01$,
$S=0.1$, $\tau=0.01,~{}N_{x}=N_{y}=256$. The values of the initial data
$u_{0}$ are random between $-1$ and $1$ .

Figure 4. Dynamics of 2D Cahn-Hilliard equation using EI scheme (1.4) where
$\nu=0.001$, $S=0.1$, $\tau=0.1,~{}N_{x}=N_{y}=256$ and the initial data
$u_{0}=0.05\sin(x)\sin(y)$.

Figure 5. Dynamics of 2D Cahn-Hilliard equation using EI scheme (1.4) where
$\nu=0.01$, $S=0.1$, $\tau=0.01,~{}N_{x}=N_{y}=256$. The values of the initial
data $u_{0}$ are random between $-1$ and $1$.

##  6\. Concluding remark

To conclude, we give a systematic approach on applying EI-type schemes to
models where maximum principle no longer exists by solving the Cahn-Hilliard
equation with a first order EI scheme and showing the energy dissipation. In
fact second order EI schemes will be handled similarly in a subsequent paper.
In particular, we prove the energy dissipation of EI schemes for the Cahn-
Hilliard equation without assuming any strong Lipschitz condition or
$L^{\infty}$ boundedness. Furthermore, we also analyze the $L^{2}$ error and
present some numerical simulations to demonstrate the dynamics. Indeed, the
analysis framework can be applied to more general models and higher order
schemes. We leave the discussion in a subsequent work.

##  Appendix A Proof of Lemma 2.7

######  Lemma A.1 ($H^{k}$ boundedness of the exact CH solution).

Assume $u(x,t)$ is a smooth solution to the Cahn-Hilliard equation in
$\mathbb{T}^{2}$ and the initial data $u_{0}\in H^{k}(\mathbb{T}^{2})$ for
$k\geq 2$. Then,

$\sup_{t\geq 0}\|u(t)\|_{H^{k}(\mathbb{T}^{2})}\lesssim_{k}1$ (A.1)

where we omit the dependence on $\nu$ and $u_{0}$.

###### Proof.

The proof is standard and we only sketch the details. To start with we can
write the solution $u$ in the mild form

$u(t)=e^{-\nu
t\Delta^{2}}u_{0}+\int_{0}^{t}e^{-\nu(t-s)\Delta^{2}}\Delta(u^{3}-u)\ ds\ .$

Note that by the energy dissipation we indeed have
$\|u\|_{H^{1}(\mathbb{T}^{2})}\lesssim 1$ for any $t>0$. Therefore it suffices
to this argument inductively, namely we show $\|u\|_{H^{2}}\lesssim 1$ for any
$t\geq 0$. Then by taking any second order spatial derivative and $L^{2}$ norm
in the formula above, we derive

$\displaystyle\|D^{2}u\|_{2}\leq\|D^{2}e^{-\nu
t\Delta^{2}}u_{0}\|_{2}+\int_{0}^{t}\|D^{2}e^{-\nu(t-s)\Delta^{2}}\Delta(u^{3}-u)\|_{2}\
ds,$

where $D^{m}u$ denotes any differential operator $D^{\alpha}u$ for any
$|\alpha|=m$.

Firstly, we consider the nonlinear part. We indeed observe that

$\displaystyle\|D^{2}e^{-\nu(t-s)\Delta^{2}}\Delta(u^{3}-u)\|_{2}\lesssim\|K_{1}*(u^{3}-u)\|_{2},$

where $K_{1}$ is the kernel corresponding to $\Delta
D^{2}e^{-\nu(t-s)\Delta^{2}}$. (It is easy to see that here $D$ can
interchange with $\Delta$.) Let $\gamma=t-s$, it then suffices to estimate the
following quantity:

$\displaystyle\int_{0}^{t}\|K_{1}*(u^{3}-u)\|_{2}\
d\gamma=\int_{0}^{1}\|K_{1}*(u^{3}-u)\|_{2}\
d\gamma+\int_{1}^{t}\|K_{1}*(u^{3}-u)\|_{2}\ d\gamma.$ (A.2)

Here we assume $t>1$ with no loss since the case $t<1$ follows easily as a
special case. For the region when $\gamma\geq 1$ we get

$\displaystyle\|K_{1}*(u^{3}-u)\|_{2}$
$\displaystyle\lesssim\|K_{1}*(u^{3}-u)\|_{\infty}$
$\displaystyle\lesssim\|K_{1}\|_{2}\cdot\|u^{3}-u\|_{2}$
$\displaystyle\lesssim\|K_{1}\|_{2}\cdot\|u\|_{6}^{3}$
$\displaystyle\lesssim\|K_{1}\|_{2},$

by the standard Sobolev embedding $\|u\|_{6}\lesssim\|u\|_{H^{1}}\lesssim 1$.
Note that for $\gamma>1$ we have

$\displaystyle\|K_{1}\|_{2}$ $\displaystyle\lesssim\left(\sum_{|k|\geq
1}|k|^{8}e^{-2\nu\gamma|k|^{4}}\right)^{\frac{1}{2}}$ (A.3)
$\displaystyle\lesssim\left(\int_{1}^{\infty}e^{-2\nu\gamma r^{4}}r^{9}\
dr\right)^{\frac{1}{2}}$
$\displaystyle\lesssim\gamma^{-\frac{1}{2}}e^{-\nu\gamma}.$

As a result, we have

$\displaystyle\int_{1}^{t}\|K_{1}*(u^{3}-u)\|_{2}\
d\gamma\lesssim\int_{1}^{t}\gamma^{-\frac{1}{2}}e^{-\nu\gamma}\
d\gamma\lesssim 1.$ (A.4)

The arguments above indeed indicate a smoothing effect. The region
$0<\gamma<1$ follows from the standard local theory. This shows
$\int_{0}^{t}\|De^{\nu(t-s)\Delta}u\|_{2}\ ds\lesssim 1$; the linear part
follows from similar arguments. As a remark our proof works at least for
dimension $d=1,2,3$.

∎

## Acknowledgement

X. Cheng was partially supported by the Shanghai “Super Postdoc” Incentive
Plan (No. 2021014), the International Postdoctoral Exchange Fellowship Program
(No. YJ20220071) and China Postdoctoral Science Foundation (Grant No.
2022M710796, 2022T150139).

## References

  * [1] N. Alikakos, P. Bates, and X. Chen.  Convergence of the Cahn-Hilliard equation to the Hele-Shaw model.  Arch. Ration. Mech. An., 128(2): 165–205, 1994.
  * [2] S. M. Allen and J. W. Cahn.  Ground state structures in ordered binary alloys with second neighbor interactions.  Acta Metallurgica., 20(3): 423–433, 1972.
  * [3] G. Bai, B. Li, and Y. Wu.  A constructive low-regularity integrator for the 1d cubic nonlinear Schrödinger equation under the Neumann boundary condition.  IMA J. Numer. Anal., to appear.
  * [4] A. Bertozzi, N. Ju, and H.-W. Lu.  A biharmonic-modified forward time stepping method for fourth order nonlinear diffusion equations.  Discrete Cont. Dyn. Syst., 29(4): 1367–1391, 2011\.
  * [5] J. Bourgain and D. Li.  Strong ill-posedness of the incompressible Euler equation in borderline Sobolev spaces.  Invent. Math., 201(1): 97–157, 2015.
  * [6] J. Bourgain and D. Li.  Strong illposedness of the incompressible Euler equation in integer $C^{m}$ spaces.  Geom. Funct. Anal., 25(1): 1–86, 2015.
  * [7] J. W. Cahn and J. E. Hilliard.  Free Energy of a Nonuniform System. I. Interfacial Free Energy.  J. Chem. Phys., 28(2): 258–267, 1958.
  * [8] L.-Q. Chen and J. Shen.  Applications of semi-implicit Fourier-spectral method to phase field equations.  Comput. Phys. Comm., 108(2–3): 147–158, 1998.
  * [9] X. Cheng. Energy stable semi-implicit schemes for the 2D Allen-Cahn and fractional Cahn-Hilliard equations. Preprint.
  * [10] A. Christlieb, J. Jones, K. Promislow, B. Wetton, and M. Willoughby.  High accuracy solutions to energy gradient flows from material science models.  J. Comput. Phys., 257: 193–215, 2014.
  * [11] J M. Church, Z. Guo, P K. Jimack, A. Madzvamuse, K. Promislow, B. Wetton, S M. Wise, F. Yang.  High accuracy benchmark problems for Allen-Cahn and Cahn-Hilliard dynamics.  Communications in computational physics., 26(4): 947-972, 2019.
  * [12] D. J. Eyre.  Unconditionally gradient stable time marching the Cahn-Hilliard equation.  MRS Proc. 529(39), 1998.
  * [13] M. Fei, F. Lin, W. Wang, and Z. Zhang.  Matrix-valued Allen–Cahn equation and the Keller–Rubinstein–Sternberg problem.  Inventiones mathematicae, 1-80, 2023.
  * [14] X. Feng and A. Prohl.  Error analysis of a mixed finite element method for the Cahn-Hilliard equation.  Numerische Mathematik., 99(1): 47–84, 2004.
  * [15] G. Gavish, J. Jones, Z. X. Xu, A. Christlieb, and K. Promislow.  Variational Models of Network Formation and Ion Transport: Applications to Perfluorosulfonate Ionomer Membranes.  Polymers., 4(1): 630–655, 2012.
  * [16] T. Ilmanen.  Convergence of the Allen-Cahn equation to Brakke’s motion by mean curvature.  J. Differ. Geom., 38: 417–461, 1993.
  * [17] B. Li.  Maximal regularity of multistep fully discrete finite element methods for parabolic equations.  IMA J. Numer. Anal., 42: 1700–1734, 2022.
  * [18] B. Li, S. Ma and K. Schratz. A Semi-implicit Exponential Low-Regularity Integrator for the Navier–Stokes Equations. SIAM J. Numer. Anal., 60(4): 2273-2292, 2022.
  * [19] D. Li.  Effective maximum principles for spectral methods.  Ann. Appl. Math., 37:131–290, 2021.
  * [20] D. Li.  Why large time-stepping methods for the Cahn-Hilliard equation is stable.  Math. Comp., 91: 2501–2515, 2022.
  * [21] Q. Liu, J. Jing, M. Yuan and W. Chen. A Positivity-preserving, energy stable BDF2 scheme with variable steps for the Cahn–Hilliard equation with logarithmic potential. J. Sci. Comput., 95(37), 2023.
  * [22] R. Pego.  Front migration in the nonlinear Cahn-Hilliard equation.  Proc. Roy. Soc. London., A(422): 261–278, 1989.
  * [23] F. Rousset and K. Schratz. A general framework of low regularity integrators, SIAM J. Numer. Anal., 59(3): 1735–1768, 2021.
  * [24] H. N. Wu and X. Yuan. Breaking quadrature exactness: A spectral method for the Allen–Cahn equation on spheres. Preprint, arXiv:2305.04820, 2023.