<figure><img src="image_1.png" /></figure>

HAL
open science

# Highly efficient schemes for time-fractional Allen-Cahn equation using extended SAV approach

Dianming Hou, Hongyi Zhu, Chuanju Xu

► To cite this version:

Dianming Hou, Hongyi Zhu, Chuanju Xu. Highly efficient schemes for time-fractional Allen-Cahn equation using extended SAV approach. Numerical Algorithms, 2021, 88 (3), pp.1077-1108. 10.1007/s11075-021-01068-y. hal-03480635

HAL Id: hal-03480635

https://hal.science/hal-03480635

Submitted on 14 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
HIGHLY EFFICIENT SCHEMES FOR TIME FRACTIONAL ALLEN-CAHN EQUATION USING EXTENDED SAV APPROACH*

DIANMING HOU<fnref n="1" /> HONGYI ZHU<fnref n="2" /> CHUANJU XU<fnref n="2" />,<fnref n="3" />,<fnref n="4" />

ABSTRACT. In this paper, we propose and analyze high order efficient schemes for the time fractional Allen-Cahn equation. The proposed schemes are based on the L1 discretization for the time fractional derivative and the extended scalar auxiliary variable (SAV) approach developed very recently to deal with the nonlinear terms in the equation. The main contributions of the paper consist in: 1) constructing first and higher order unconditionally stable schemes for different mesh types, and proving the unconditional stability of the constructed schemes for the uniform mesh; 2) carrying out numerical experiments to verify the efficiency of the schemes and to investigate the coarsening dynamics governed by the time fractional Allen-Cahn equation. Particularly, the influence of the fractional order on the coarsening behavior is carefully examined. Our numerical evidence shows that the proposed schemes are more robust than the existing methods, and their efficiency is less restricted to particular forms of the nonlinear potentials.

1. INTRODUCTION

As a useful modelling tool, gradient flows has been used to model many physical problems, particularly dissipative systems, which are systems driven by dissipation of free energy. This variational modelling framework generally leads to partial differential equations having the mathematical structure of gradient flows. In its abstract form, variational modelling of such systems consists of choosing a state space, a driving functional, and a dissipation mechanism. Precisely, models of gradient flows take the general form:

$$ \frac{\partial \phi}{\partial t} = -\operatorname{grad}_H E(\phi), \qquad (1.1) $$

where $\phi$ is the state function (also called phase function in many cases), $E[\phi(x, t)]$ is the free energy driving functional associated to the physical problem, and $\operatorname{grad}_H E(\cdot)$ is the functional derivative of $E$ in the Sobolev space $H$. Multiplying both sides of (1.1) by $\delta E/\delta \phi$ and integrating

<footnote>Date: November 12, 2020.</footnote>

<footnote n="2020">2020 Mathematics Subject Classification. 65N35, 65M70, 45K05, 41A05, 41A10, 41A25.</footnote>

<footnote>Key words and phrases. time fractional Allen-Cahn, time-stepping scheme, unconditional stability, spectral method.</footnote>

<footnote>*This research is partially supported by NSF of China (Grant numbers 11971408, 51661135011, 11421110001, and 91630204).</footnote>

<footnote n="1">¹School of Mathematics and Statistics, Jiangsu Normal University, 221116 Xuzhou, China.</footnote>

<footnote n="2">²School of Mathematical Sciences and Fujian Provincial Key Laboratory of Mathematical Modeling and High Performance Scientific Computing, Xiamen University, 361005 Xiamen, China.</footnote>

<footnote n="3">³Bordeaux INP, Laboratoire I2M UMR 5295, 33607 Pessac, France.</footnote>

<footnote n="4">⁴Corresponding author. Email: cjxu@xmu.edu.cn (C. Xu).</footnote>
the resulting equation gives the energy dissipation law:

$$ \frac{d}{dt} E(\phi) = (\mathrm{grad}_H E(\phi), \frac{\partial \phi}{\partial t}) = -\|\mathrm{grad}_H E(\phi)\|_0^2, \quad (1.2) $$

where $(\cdot, \cdot)$ is the $L^2$-inner product, and $\| \cdot \|_0$ is the associated norm. This means that the state function $\phi$ evolves in such a way as to have the energy functional $E$ decrease in time (least action principle from physics), i.e., in the opposite direction to the gradient of $E$ at $\phi$.

The energy-based variational framework makes the equations a thermodynamically-consistent and physically attractive approach to model multi-phase flows. It has long been used in many fields of science and engineering, particularly in materials science and fluid dynamics, see, e.g., [10, 6, 24, 13, 19, 7, 16, 39] and the references therein. Typical examples include the Cahn-Hilliard and Allen-Cahn equations for multi-phase flows, for which the evolution PDE system is resulted from the energetic variation of the action functional of the total free energy in different Sobolev spaces.

In this paper, we are interested in models deriving from gradient flows having modified dissipation mechanism. That is, we consider the gradient flows of the form

$$ C_0^\alpha D_t^\alpha \phi = -\mathrm{grad}_H E(\phi), \quad (1.3) $$

where $0 < \alpha < 1$, and $C_0^\alpha D_t^\alpha$ is the Caputo fractional derivative defined as follows:

$$ C_0^\alpha D_t^\alpha \phi(t) := \frac{1}{\Gamma(1-\alpha)} \int_0^t (t-s)^{-\alpha} \phi'(s) ds. $$

It is seen from the definition that the fractional derivative is some kind of weighted average of the traditional derivative in the history. This means that the change rate, i.e., the derivative, at the current time is affected by the historical rates. This property is often used to describe the memory effect which can be present in some materials such as viscoelastic materials or polymers. Intuitively, the gradient flows model (1.1) is useful for describing the systems in which dissipation of the associated free energy has memory effect in some circumstances.

Mathematically, such fractional-type gradient flows have been studied in the last few years and have quickly attracted more and more attentions [12, 44, 14, 1, 2, 3, 42, 28, 34, 25, 21]. Wang et al. [42, 28] numerically studied the time fractional phase field models, including the Cahn-Hilliard equations with different variable mobilities and molecular beam epitaxy models. Their numerical tests indicated that the effective free energy and roughness of the time fractional phase field models obey a universal power law scaling dynamics during coarsening. Tang et al. [36] proved that the time-fractional phase field models indeed admit a modified energy dissipation law of an integral type. Moreover, they applied the L1 approach to discretize the time fractional derivative on the uniform mesh, leading to a scheme of first order accuracy. Recently, Du et al. [15] developed several time stepping schemes, i.e., convex splitting scheme, weighted convex splitting scheme, and linear weighted stabilized scheme, for the time fractional Allen-Cahn equation. They proved that the convergence rates of the proposed time-stepping schemes are $O(\tau^\alpha)$ for typical solutions of the underlying equation. Very recently, Liao et al. [26] developed a second-order nonuniform time-stepping scheme for the time-fractional Allen-Cahn equation. They showed that the proposed scheme preserves the discrete maximum principle and presented a sharp maximum-norm error estimate which reflects the temporal regularity.
The goal of this paper is to construct and analyze some more efficient schemes for the time fractional Allen-Cahn equation as follows:

$$C_0 D_t^\alpha \phi - \varepsilon^2 \Delta \phi + F'(\phi) = 0. \qquad (1.4)$$

Precisely, our aim is to propose unconditionally stable schemes which have provable higher order convergence rate than the existing schemes in the literature. The main idea is to combine existing efficient discretizations for the time fractional derivative and SAV approaches for the traditional gradient flows. In fact, there has been significant progress in the numerical development for both the time fractional differential equation and gradient flows in recent years. Among the existing approaches for the time fractional derivative, the so-called L1 scheme is probably the best known; see, e.g., Sun and Wu [35] and Lin and Xu [27]. The L1 scheme makes use of a piecewise linear approximation and attains $(2 - \alpha)$-order convergence rate for the $\alpha$-order fractional derivative [27, 18, 29]. Based on piecewise linear approximation at the closest interval to the time instant and piecewise quadratic interpolation at the previous time intervals, Alikhanov [5] constructed a L2-1$_\sigma$ scheme and proved it is second order accurate when applied to the time fractional diffusion equation. Gao et al. [17] and Lv and Xu [30] developed $(3 - \alpha)$-order schemes by using piecewise quadratic interpolation. In order to reduce the high storage requirement, several authors have proposed to approximate the weakly singular kernel function by a sum-of-exponentials (SOE), resulting in some new stepping schemes with reduced storage; see, e.g., [8, 9, 22, 41, 37, 40]. Concerning numerical approaches for gradient flows, we would like to mention the most recent results in the past few years. There appeared two novel strategies: the invariant energy quadratization (IEQ) method [38, 43] and the so-called scalar auxiliary variable (SAV) approach [32, 33]. Based on these two approaches, second-order unconditionally stable schemes have been successfully constructed for a large class of gradient flow models. Very recently, Hou et al. [20] proposed an extension of the SAV approach, which leads to more robust unconditionally stable schemes and achieves more ability to guarantee the dissipation law of the total free energy.

The purpose of the current study is twofold. The first is to make use of the above mentioned results to construct efficient numerical schemes for the time fractional Allen-Cahn equation, and analyze the stability properties of the proposed methods. The convergence behavior will be investigated by means of numerical examples. The second is to carry out numerical experiments by using the proposed schemes to study coarsening dynamics governed by the fractional Allen-Cahn equation. In particular, we are interested to the effect of the fractional order on the coarsening behavior. The paper is organized as follows: In the next section, we propose the extended SAV reformulation for the fractional Allen-Cahn equation. In Section 3, we describe the numerical schemes for both uniform and nonuniform meshes, and carry out the stability analysis for the proposed schemes on the uniform mesh. In Section 4, we construct and analyze a higher order scheme. The numerical experiments are carried out in Section 5, not only to validate the proposed methods, but also to numerically investigate the coarsening dynamics. Finally, the paper ends with some concluding remarks.

## 2. EXTENDED SAV REFORMULATION

We first introduce some notations which will be used throughout the paper. Let $H^m(\Omega)$ and $\| \cdot \|_m, m = 0, \pm 1, \dots$, denote the standard Sobolev spaces and their norms, respectively.
In particular, the norm and inner product of $L^2(\Omega) := H^0(\Omega)$ are denoted by $\| \cdot \|_0$ and $(\cdot, \cdot)$ respectively.

For the sake of simplicity we consider the time fractional Allen-Cahn equation on a bounded domain $\Omega \in \mathbb{R}^n$ ($n = 1, 2, 3$) subject to the periodic boundary condition or Neumann boundary condition, i.e.,

$$ 
\begin{aligned}
& C_0^\alpha D_t^\alpha \phi - \varepsilon^2 \Delta \phi + F'(\phi) = 0, \text{ in } (0, T] \times \Omega; \\
& \frac{\partial \phi}{\partial n}|_{\partial \Omega} = 0 \text{ or periodic boundary condition;} \\
& \phi(\boldsymbol{x}, 0) = \phi_0(\boldsymbol{x}), \text{ in } \Omega.
\end{aligned} 
\quad (2.1) $$

It has been proved in [36] that this problem admits the following dissipation law:

$$ E(\phi(\boldsymbol{x}, t)) - E(\phi(\boldsymbol{x}, 0)) = - \int_{\Omega} \mathcal{A}_{\alpha}(\phi_t, \phi_t) d\boldsymbol{x} \le 0, \text{ for all } t \ge 0, \quad (2.2) $$

where

$$ E(\phi) := \int_{\Omega} \left[ \frac{\varepsilon^2}{2} |\nabla \phi|^2 + F(\phi) \right] d\boldsymbol{x}, \quad \mathcal{A}_{\alpha}(g, h) := \frac{1}{\Gamma(1-\alpha)} \int_0^t \int_0^s (s-\sigma)^{-\alpha} g(\sigma) h(s) d\sigma ds. $$

We also define the partial free energy $E_{\theta}(\phi)$ by

$$ E_{\theta}(\phi) := \int_{\Omega} \left[ \frac{\theta}{2} |\nabla \phi|^2 + F(\phi) \right] d\boldsymbol{x}, \quad 0 \le \theta \le \varepsilon^2. $$

Throughout this paper, we will assume that $E_{\theta}(\phi)$ is bounded from below, i.e., there always exists a positive constant $C_0$ such that $E_{\theta}(\phi) + C_0 > 0$.

We first derive a $H^2$-regularity result for the time fractional Allen-Cahn equation (2.1) with a generalized nonlinear free energy density $F(\phi)$. For nonlinear equations the $H^2$-regularity plays a crucial role in both theoretical and numerical analysis for the reason that the boundedness of the solution can be guaranteed by the $H^2$-regularity (for equations of dimension less than four). We would like to point out that the $H^2$-regularity for the classical Allen-Cahn equation, i.e., $\alpha = 1$, has recently been established [31] by employing a type of Galerkin approach. We will make use of a similar approach to study the regularity of the fractional version.

**Theorem 2.1.** Assume $\phi_0 \in H^2(\Omega)$ and $F \in C^3(\mathbb{R})$ satisfying

$$ |F''(x)| < c(|x|^p + 1), \quad p > 0 \text{ if } n = 1, 2; \quad 0 < p < 4 \text{ if } n = 3. \quad (2.3) $$

Then for any $T > 0$, the problem (2.1) has a unique solution $\phi$ satisfying

$$ \phi \in C([0, T]; H^2(\Omega)) \cap L^2([0, T]; H^3(\Omega)). $$

*Proof.* Denote by $\{\omega_j\}$ the orthonormal basis in $L^2(\Omega)$ consisting of the eigenfunctions of $-\Delta$, i.e.

$$ -\Delta\omega_j = \lambda_j\omega_j. $$

We first consider the approximate solution under the form

$$ \phi_m(\cdot, t) = \sum_{j=1}^{m} g_{jm}(t)\omega_j $$
to the problem

$$ ( {}_0^C D_t^\alpha \phi_m, \omega_j ) + \varepsilon^2 ( \nabla \phi_m, \nabla \omega_j ) + ( F'(\phi), \omega_j ) = 0, \quad (2.4) $$

where $\phi_m(x, 0)$ is the $L^2$-projection of $\phi_0(x)$ on the space spanned by $\{\omega_j\}$.

Multiplying both sides of (2.4) by $g'_{jm}(t)$ then summing up for $j = 1, \dots, m$, we get the following energy inequality

$$ E(\phi_m(x, t)) - E(\phi_m(x, 0)) = - \int_{\Omega} \mathcal{A}_{\alpha}(\partial_t \phi_m, \partial_t \phi_m) d\mathbf{x} \le 0. $$

Under the assumption $\phi_0 \in H^2(\Omega) \subseteq L^\infty(\Omega)$, $E_\theta(\phi) + C_0 > 0$, $F \in C^3(\mathbb{R})$, and using the definition of $E(\phi)$, we deduce from the above energy inequality the following $H^1$-bound

$$ \|\phi_m\|_{H^1} \le cE(\phi_m(0)) + c \le M_0, $$

where $M_0$ is positive constant depending on $\|\phi_0\|_{H^2}$.

To derive the $H^2$-bound, we multiply (2.4) by $\lambda_j^2 g_{jm}(t)$ and sum up for $j = 1, \dots, m$. This gives

$$ \int_{\Omega} {}_0^C D_t^\alpha \Delta \phi_m \Delta \phi_m d\mathbf{x} + \|\nabla \Delta \phi_m\|_0^2 = -(\nabla F'(\phi_m), \nabla \Delta \phi_m) \le \frac{1}{2} \|\nabla F'(\phi_m)\|_0^2 + \frac{1}{2} \|\nabla \Delta \phi_m\|_0^2. \quad (2.5) $$

Using Lemma 1 in [4] and Lemma 2.3 in [31], we get

$$ \int_{\Omega} {}_0^C D_t^\alpha \Delta \phi_m \Delta \phi_m d\mathbf{x} \ge \frac{1}{2} {}_0^C D_t^\alpha \|\Delta \phi_m\|_0^2, \quad \|\nabla F'(\phi_m)\|_0^2 \le \frac{1}{2} \|\nabla \Delta \phi_m\|_0^2 + c(M_0). \quad (2.6) $$

Combining (2.5) and (2.6) gives

$$ {}_0^C D_t^\alpha \|\Delta \phi_m\|_0^2 + \frac{1}{2} \|\nabla \Delta \phi_m\|_0^2 \le c(M_0). $$

Then, applying the Riemann-Liouville fractional integral operator ${}_0I_t^\alpha$ to the above inequality, we obtain

$$ \|\Delta \phi_m\|_0^2 - \|\Delta \phi_m(0)\|_0^2 + \frac{1}{2} {}_0I_t^\alpha \|\nabla \Delta \phi_m\|_0^2 \le c(M_0) T^\alpha. $$

Furthermore, we have

$$ \|\Delta \phi_m\|_0^2 \le c(M_0) T^\alpha, \quad \int_0^t \|\nabla \Delta \phi_m\|_0^2 dt \le c(M_0)(T + T^{1-\alpha}). $$

This means $\phi_m \in L^\infty([0, T]; H^2(\Omega)) \cap L^2([0, T]; H^3(\Omega))$.

Following a similar discussion as in Theorem 2.6 in [31], we can prove $\phi_m \to \phi$ strongly in $L^2([0, T]; H^1(\Omega))$ as $m \to \infty$, and $\phi \in C([0, T]; H^2(\Omega)) \cap L^2([0, T]; H^3(\Omega))$.

The uniqueness of the solution can also be established by following the same lines as in [31], Lemma 1 in [4], and Gronwall lemma. This completes the proof. □

The key to construct our schemes for (2.1) is to rewrite the original equation into the following equivalent form, inspired by the idea of the extended SAV approach introduced in [20]:

$$ {}_0^C D_t^\alpha \phi - \varepsilon^2 \Delta \phi + \left(1 - \frac{R(t)}{\sqrt{E_\theta(\phi) + C_0}}\right) \theta \Delta \phi + \frac{R(t)}{\sqrt{E_\theta(\phi) + C_0}} F'(\phi) = 0, \quad (2.7) $$

where

$$ R(t) = \sqrt{E_{\theta}(\phi) + C_0}. $$
At the discrete level, we need to find a suitable way to evaluate the auxiliary variable $R(t)$. To this end, we formally take the derivative of $R(t)$ with respect to $t$ to obtain the following auxiliary equation:

$$ \frac{dR}{dt} = \frac{1}{2\sqrt{E_{\theta}(\phi) + C_0}} \int_{\Omega} (-\theta\Delta\phi + F'(\phi)) \frac{\partial\phi}{\partial t} dx. \quad (2.8) $$

Then we have two coupled equations (2.7) and (2.8) governing two unknown functions, i.e., the original two variable function $\phi$ and the auxiliary one variable function $R(t)$. These two equations will be used to construct our efficient time-stepping schemes to numerically compute the approximate solutions. Before that, we first realize, by taking $L^2$-inner products and integrating from 0 to $t$ of the two equations (2.7) and (2.8) with $\frac{\partial\phi}{\partial t}$ and $2R(t)$ respectively, that

$$ [R^2(t) + \frac{\varepsilon^2 - \theta}{2} \|\nabla\phi(\mathbf{x}, t)\|_0^2] - [R^2(0) + \frac{\varepsilon^2 - \theta}{2} \|\nabla\phi(\mathbf{x}, 0)\|_0^2] = - \int_{\Omega} \mathcal{A}_{\alpha}(\phi_t, \phi_t) d\mathbf{x} \le 0. \quad (2.9) $$

Noticing that

$$ R^2 + \frac{\varepsilon^2 - \theta}{2} \|\nabla\phi\|_0^2 = E(\phi) + C_0, $$

we obtain exactly the same dissipation law as (2.2). This is a natural result since the two coupled equations (2.7) and (2.8), together with corresponding boundary and initial conditions, are strictly equivalent to the original problem (2.1). However, as we are going to see, at the discrete level the energy dissipation law will take the form (2.9) more than (2.2).

### 3. A FIRST ORDER SCHEME

Starting with the equivalent equations (2.7) and (2.8), we are able to construct various time stepping schemes to calculate the solution $\phi$. Let $n = 0, 1, \dots, M$ be the time step, $\Delta t_n := t_n - t_{n-1}$ be the time step size, and $\tau = \max\{\Delta t_n, n = 1, 2, \dots\}$ be the maximum time step size of the temporal mesh.

3.1. **First order scheme.** We first recall the very popular L1 approximation [27] to discretize the Caputo fractional derivative $C_0^\alpha D_t^\alpha \phi$:

$$ \begin{aligned} C_0^\alpha D_{t_{n+1}}^\alpha \phi &= \frac{1}{\Gamma(1-\alpha)} \int_0^{t_{n+1}} (t_{n+1}-s)^{-\alpha} \phi_s(\mathbf{x}, s) ds \\ &= \sum_{k=0}^n \frac{\phi(\mathbf{x}, t_{k+1}) - \phi(\mathbf{x}, t_k)}{\Delta t_{k+1}} \frac{1}{\Gamma(1-\alpha)} \int_{t_k}^{t_{k+1}} (t_{n+1}-s)^{-\alpha} ds + e_{\tau}^{n+1} \\ &:= L_{n+1}^{\alpha} \phi + e_{\tau}^{n+1}, \end{aligned} $$

where

$$ L_{n+1}^{\alpha} \phi = \sum_{k=0}^{n} b_{n-k} \frac{\phi(\mathbf{x}, t_{k+1}) - \phi(\mathbf{x}, t_k)}{\Delta t_{k+1}}, \quad b_{n-k} = \frac{1}{\Gamma(1-\alpha)} \int_{t_k}^{t_{k+1}} (t_{n+1}-s)^{-\alpha} ds, \quad (3.1) $$

and the truncation error $e_{\tau}^{n+1}$ is defined by

$$ e_{\tau}^{n+1} = \frac{1}{\Gamma(1-\alpha)} \left[ \int_{0}^{t_{n+1}} (t_{n+1}-s)^{-\alpha} \phi_s(\mathbf{x}, s) ds - \sum_{k=0}^{n} \frac{\phi(\mathbf{x}, t_{k+1}) - \phi(\mathbf{x}, t_k)}{\Delta t_{k+1}} \int_{t_k}^{t_{k+1}} (t_{n+1}-s)^{-\alpha} ds \right]. $$
- In the uniform mesh case, i.e., $t_n = \frac{n}{M}T, n = 0, 1, \dots, M$, we have

$$b_j = \frac{\tau^{1-\alpha}}{\Gamma(2-\alpha)} [(j+1)^{1-\alpha} - j^{1-\alpha}], \quad j = 0, 1, \dots, n. \qquad (3.2)$$

In this case the truncation error $e_{\tau}^{n+1}$ can be bounded by $c_{\phi}\tau^{2-\alpha}$ if the exact solution is sufficiently smooth [27]. The following property has also been known; see, e.g., [36]: for any $(u_0, u_1, \dots, u_n)^T \in \mathbb{R}^{n+1}$, it holds

$$\sum_{k=0}^{n} \sum_{j=0}^{k} b_{k-j} u_j u_k > 0. \qquad (3.3)$$

- In the graded mesh case, i.e., $t_n = (\frac{n}{M})^r T, r \ge 1, n = 0, 1, \dots, M$, we have

$$b_j = \frac{T^{1-\alpha}}{\Gamma(2-\alpha)M^{(1-\alpha)r}} \left[ [(n+1)^r - (n-j)^r]^{1-\alpha} - [(n+1)^r - (n-j+1)^r]^{1-\alpha} \right].$$

Applying the above L1 approximation and some other first order finite difference approximations to (2.7) and (2.8), we propose the first scheme as follows:

$$L_{n+1}^{\alpha} \phi - \varepsilon^2 \Delta \phi^{n+1} + \left(1 - \frac{R^{n+1}}{\sqrt{E_{\theta}^{n} + C_{0}}\right) \theta \Delta \phi^{n} + \frac{R^{n+1}}{\sqrt{E_{\theta}^{n} + C_{0}}} F'(\phi^{n}) = 0, \qquad (3.4a)$$

$$\frac{R^{n+1} - R^n}{\Delta t_{n+1}} = \frac{1}{2\sqrt{E_{\theta}^n + C_0}} \int_{\Omega} (-\theta \Delta \phi^n + F'(\phi^n)) \frac{\phi^{n+1} - \phi^n}{\Delta t_{n+1}} dx, \qquad (3.4b)$$

where $\phi^n$ is an approximation to $\phi(t^n)$ and $E_{\theta}^n := E_{\theta}(\phi^n)$. Intuitively, this scheme is of first-order convergence with respect to the time step size. Although a rigorous proof remains an open question, this intuitive convergence rate will be investigated through the numerical experiments to be presented later on.

3.2. **Stability.** Next we will analyze the stability property of the scheme in the following theorem.

**Theorem 3.1.** If the mesh is uniform, then the scheme (3.4) is unconditionally stable in the sense that the following discrete energy law holds

$$\frac{E_{\varepsilon, \theta}^{n+1} - E_{\varepsilon, \theta}^{0}}{\tau} \le - \sum_{k=0}^{n} (L_{k+1}^{\alpha} \phi, \frac{\phi^{k+1} - \phi^k}{\tau}) < 0, \qquad (3.5)$$

where

$$E_{\varepsilon, \theta}^{n} := \frac{1}{2}(\varepsilon^2 - \theta) \|\nabla \phi^{n}\|_{0}^{2} + |R^{n}|^{2}.$$

*Proof.* For the uniform mesh, we deduce from taking the $L^2$-inner products of (3.4a) and (3.4b) with $\frac{\phi^{n+1}-\phi^n}{\tau}$ and $2R^{n+1}$ respectively:

$$\left(L_{n+1}^{\alpha} \phi, \frac{\phi^{n+1} - \phi^n}{\tau}\right) + \frac{\varepsilon^2 - \theta}{\tau} \left(\nabla \phi^{n+1}, \nabla (\phi^{n+1} - \phi^n)\right) + \frac{\theta}{\tau} \|\nabla (\phi^{n+1} - \phi^n)\|_0^2 \\ + \frac{R^{n+1}}{\sqrt{E_{\theta}^n + C_0}} \left(-\theta \Delta \phi^n + F'(\phi^n), \frac{\phi^{n+1} - \phi^n}{\tau}\right) = 0, \\ \frac{2R^{n+1}(R^{n+1} - R^n)}{\tau} = \frac{R^{n+1}}{\sqrt{E_{\theta}^n + C_0}} \left(-\theta \Delta \phi^n + F'(\phi^n), \frac{\phi^{n+1} - \phi^n}{\tau}\right).$$
Then using the identity

$$2a^{n+1}(a^{n+1} - a^n) = |a^{n+1}|^2 - |a^n|^2 + |a^{n+1} - a^n|^2. \qquad (3.6)$$

and dropping some non-essential positive terms, we obtain

$$\frac{E_{\varepsilon, \theta}^{n+1} - E_{\varepsilon, \theta}^{n}}{\tau} \leq - \left( L_{n+1}^{\alpha} \phi, \frac{\phi^{n+1} - \phi^n}{\tau} \right), \quad n = 0, 1, 2, \dots,$$

Summing up the above equality from 0 to n and using (3.1) and (3.3), we obtain (3.5). □

**Remark 3.1.** Theorem 3.1 indicates that the scheme (3.4) is unconditionally stable in the sense that the modified energy $E_{\varepsilon, \theta}^n$ is a decreasing function of the time step. Unfortunately, as all other SAV-based schemes, there is no proof for the decay of the original energy $E(\phi^n)$, although many (but not all) numerical investigation have shown that the discrete original energy is also decreasing. Our numerical experiments will also confirm this point. It is notable that $E_{\varepsilon, \theta}^n$ is a first order approximation to $E(\phi^n)$. Furthermore, the energy stability (3.5) implies the $H^1$-boundedness of $\phi^n$, i.e., there exists a positive constant $M_0$, which may depend on $\Omega$ and $\phi_0$, such that

$$\|\phi^n\|_{H^1} + |R^n| \leq M_0, \text{ for all } n = 1, 2, \dots, \text{ if } \phi_0 \in H^1(\Omega). \qquad (3.7)$$

The next theorem gives the $H^2$ bound, consequently $L^\infty$ bound, of the numerical solutions.

**Theorem 3.2.** Assume $\phi_0 \in H^3(\Omega)$ and $F \in C^3(\mathbb{R})$ satisfying (2.3). Let $M_0$ be the bound given in (3.7). Then for the uniform mesh and any $(n+1) < T/\tau$, the scheme (3.4) satisfies

$$\sum_{k=0}^{n} b_{n-k} \|\Delta\phi^{k+1}\|_{0}^{2} + \tau\varepsilon^{2} \|\nabla\Delta\phi^{n+1}\|_{0}^{2} \leq c(M_{0}, \varepsilon)T + \frac{T^{1-\alpha}}{\Gamma(2-\alpha)} \|\Delta\phi_{0}\|_{0}^{2} + \tau\varepsilon^{2} \|\nabla\Delta\phi_{0}\|_{0}^{2}.$$

*Proof.* We deduce from taking the $L^2$-inner products of (3.4a) with $2\tau\Delta^2\phi^{n+1}$,

$$2\tau \left( L_{n+1}^{\alpha} \Delta\phi, \Delta\phi^{n+1} \right) + 2\tau\varepsilon^2 \|\nabla\Delta\phi^{n+1}\|_0^2 = 2\tau \left( 1 - \frac{R^{n+1}}{\sqrt{E_{\theta}^n + C_0}} \right) \theta(\nabla\Delta\phi^n, \nabla\Delta\phi^{n+1}) \\ + \frac{2\tau R^{n+1}}{\sqrt{E_{\theta}^n + C_0}} (\nabla F'(\phi^n), \nabla\Delta\phi^{n+1}).$$

For the terms in the right hand side, with a similar discussion as in [31] Lemma 2.4, we can obtain that

$$2\tau \left( L_{n+1}^{\alpha} \Delta\phi, \Delta\phi^{n+1} \right) + 2\tau\varepsilon^2 \|\nabla\Delta\phi^{n+1}\|_0^2 \leq \tau c(M_0, \varepsilon) + \tau\varepsilon^2 \|\nabla\Delta\phi^n\|_0^2 + \tau\varepsilon^2 \|\nabla\Delta\phi^{n+1}\|_0^2. \quad (3.8)$$

Moreover, using Lemma 1 in [5] gives

$$(L_{n+1}^{\alpha} \Delta \phi, \Delta \phi^{n+1}) \geq \frac{1}{2} L_{n+1}^{\alpha} \| \Delta \phi \|_{0}^{2}. \qquad (3.9)$$

Combining (3.8) and (3.9), we have

$$\tau L_{n+1}^{\alpha} \|\Delta\phi\|_{0}^{2} + \tau\varepsilon^{2} \|\nabla\Delta\phi^{n+1}\|_{0}^{2} \leq \tau c(M_{0}, \varepsilon) + \tau\varepsilon^{2} \|\nabla\Delta\phi^{n}\|_{0}^{2}.$$

Taking the sum of the above inequality from 0 to n and noticing that

$$\sum_{k=0}^{n} \tau L_{k+1}^{\alpha} \|\Delta\phi\|_{0}^{2} = \sum_{k=0}^{n} b_{n-k} \|\Delta\phi^{k+1}\|_{0}^{2} - \sum_{k=0}^{n} b_{k} \|\Delta\phi_{0}\|_{0}^{2} = \sum_{k=0}^{n} b_{n-k} \|\Delta\phi^{k+1}\|_{0}^{2} - \frac{[(n+1)\tau]^{1-\alpha}}{\Gamma(2-\alpha)} \|\Delta\phi_{0}\|_{0}^{2},$$
we get

$$ \sum_{k=0}^{n} b_{n-k} \|\Delta\phi^{k+1}\|_{0}^{2} + \tau\varepsilon^{2} \|\nabla\Delta\phi^{n+1}\|_{0}^{2} \le (n+1)\tau c(M_{0}, \varepsilon) + \frac{[(n+1)\tau]^{1-\alpha}}{\Gamma(2-\alpha)} \|\Delta\phi_{0}\|_{0}^{2} + \tau\varepsilon^{2} \|\nabla\Delta\phi_{0}\|_{0}^{2} \\ \le c(M_{0}, \varepsilon)T + \frac{T^{1-\alpha}}{\Gamma(2-\alpha)} \|\Delta\phi_{0}\|_{0}^{2} + \tau\varepsilon^{2} \|\nabla\Delta\phi_{0}\|_{0}^{2}. $$

The proof is completed.

<figure><img src="image_3.png" /></figure>

3.3. **Implementation.** Beside of its unconditional stability, another notable property of the new scheme is that it can be efficiently implemented. To see that, let's consider the case of uniform mesh. We first eliminate $R^{n+1}$ from (3.4a) by using (3.4b) to obtain

$$ \sum_{k=0}^{n} b_{n-k} \frac{\phi^{k+1} - \phi^{k}}{\tau} - \varepsilon^2 \Delta\phi^{n+1} + \theta \Delta\phi^{n} + \left[ \frac{R^n}{\sqrt{E_{\theta}^n + C_0}} + \frac{1}{2(E_{\theta}^n + C_0)} (\gamma^n, \phi^{n+1} - \phi^n) \right] \gamma^n = 0, \quad (3.10) $$

where

$$ \gamma^n := -\theta \Delta \phi^n + F'(\phi^n). $$

Then reformulating (3.10) gives

$$ \left( \frac{b_0}{\tau} - \varepsilon^2 \Delta \right) \phi^{n+1} + (\gamma^n, \phi^{n+1}) \frac{\gamma^n}{2(E_{\theta}^n + C_0)} \quad (3.11) $$

$$ = -\frac{1}{\tau} \left[ \sum_{k=0}^{n-1} (b_{n-k} - b_{n-k-1})\phi^{k+1} + b_n\phi^0 \right] - \theta\Delta\phi^n - \left[ \frac{R^n}{\sqrt{E_{\theta}^n + C_0}} - \frac{1}{2(E_{\theta}^n + C_0)}(\gamma^n, \phi^n) \right] \gamma^n. $$

Denoting the right hand side in (3.11) by $g(\phi^n)$, we see that problem (3.11) can be solved in two steps as follows:

$$ \left\{ \begin{array}{l} \left( \frac{b_0}{\tau} - \varepsilon^2 \Delta \right) \phi_1^{n+1} = - \frac{\gamma^n}{2(E_\theta^n + C_0)}, \\ \text{Neumann boundary condition or periodic boundary condition on } \phi_1^{n+1}; \end{array} \right. \qquad (3.12a) $$

$$ \left\{ \begin{array}{l} \left( \frac{b_0}{\tau} - \varepsilon^2 \Delta \right) \phi_2^{n+1} = g(\phi^n), \\ \text{Neumann boundary condition or periodic boundary condition on } \phi_2^{n+1}; \end{array} \right. \qquad (3.12b) $$

$$ \phi^{n+1} = (\gamma^n, \phi^{n+1})\phi_1^{n+1} + \phi_2^{n+1}. \qquad (3.12c) $$

In a first look it seems that (3.12c) satisfied by $\phi^{n+1}$ is an implicit equation. However a careful examination shows that $(\gamma^n, \phi^{n+1})$ in (3.12c) can be determined explicitly. In fact, taking the inner product of (3.12c) with $\gamma^n$ yields

$$ (\gamma^n, \phi^{n+1}) + \sigma^n(\gamma^n, \phi^{n+1}) = (\gamma^n, \phi_2^{n+1}), \qquad (3.13) $$

where

$$ \sigma^n = -(\gamma^n, \phi_1^{n+1}) = (\gamma^n, A^{-1} \frac{\gamma^n}{2(E_\theta^n + C_0)}) \quad \text{with} \quad A = \frac{b_0}{\tau} - \varepsilon^2 \Delta. \qquad (3.14) $$

Note that $A^{-1}$ is a positive definite operator and $E_\theta^n + C_0 > 0$. Thus $\sigma^n \ge 0$. Then it follows from (3.13) that

$$ (\gamma^n, \phi^{n+1}) = \frac{(\gamma^n, \phi_2^{n+1})}{1 + \sigma^n}. \qquad (3.15) $$
Using this expression, $\phi^{n+1}$ can be explicitly computed from (3.12c).

In detail, the scheme (3.4) results in the following algorithm at each time step:

(i) Calculation of $\phi_1^{n+1}$ and $\phi_2^{n+1}$: solving the elliptic problems (3.12a) and (3.12b) respectively, which can be realized in parallel.

(ii) Evaluation of $(\gamma^n, \phi^{n+1})$ using (3.15), then $\phi^{n+1}$ using (3.12c).

Thus the overall computational complexity at each time step is essentially solving two second-order elliptic problems with constant coefficients, for which there exist different fast solvers depending on the spatial discretization method.

## 4. A HIGHER ORDER SCHEME

### 4.1. 2-α order scheme.
Here we propose a semi-implicit scheme, called L1-CN scheme hereafter, by applying the L1 discretization to the fractional derivative and the Crank-Nicolson discretization to the remaining terms for the system (2.7) and (2.8):

$$L_{n+\frac{1}{2}}^{\alpha} \phi - \varepsilon^2 \frac{\Delta(\phi^{n+1} + \phi^n)}{2} + \left(1 - \frac{R^{n+1} + R^n}{2\sqrt{E_{\theta}^{n+\frac{1}{2}} + C_0}}\right) \theta \Delta \phi^{n+\frac{1}{2}} + \frac{R^{n+1} + R^n}{2\sqrt{E_{\theta}^{n+\frac{1}{2}} + C_0}} F'(\phi^{n+\frac{1}{2}}), \quad (4.1a)$$

$$\frac{R^{n+1} - R^n}{\Delta t_{n+1}} = \frac{1}{2\sqrt{E_{\theta}^{n+\frac{1}{2}} + C_0}} \int_{\Omega} (-\theta \Delta \phi^{n+\frac{1}{2}} + F'(\phi^{n+\frac{1}{2}})) \frac{\phi^{n+1} - \phi^n}{\Delta t_{n+1}} dx, \quad (4.1b)$$

where $\phi^{n+\frac{1}{2}} := \phi^n + \frac{\Delta t_{n+1}}{2\Delta t_n} [\phi^n - \phi^{n-1}]$ is an explicit approximation to $\phi(t^{n+\frac{1}{2}})$, $E_{\theta}^{n+\frac{1}{2}} := E_{\theta}(\phi^{n+\frac{1}{2}})$, for $n = 0, 1, \dots$ with $\phi^{-1} := \phi^0$, and

$$\begin{aligned} L_{n+\frac{1}{2}}^{\alpha} \phi &:= \sum_{k=0}^{n-1} \frac{\phi^{k+1} - \phi^k}{\Delta t_{k+1}} \frac{1}{\Gamma(1-\alpha)} \int_{t_k}^{t_{k+1}} (t_{n+\frac{1}{2}} - s)^{-\alpha} ds \\ &\quad + \frac{\phi^{n+1} - \phi^n}{\Delta t_{n+1}} \frac{1}{\Gamma(1-\alpha)} \int_{t_n}^{t_{n+\frac{1}{2}}} (t_{n+\frac{1}{2}} - s)^{-\alpha} ds \\ &:= \sum_{k=0}^{n} \tilde{b}_{n-k} \frac{\phi^{k+1} - \phi^k}{\Delta t_{k+1}}, \quad \text{with } t_{n+\frac{1}{2}} = (t_n + t_{n+1})/2. \end{aligned} \quad (4.2)$$

- In the uniform mesh case, the coefficient $\tilde{b}_k$ can be expressed as follows:

$$\tilde{b}_0 = \frac{\tau^{1-\alpha}}{\Gamma(2-\alpha)2^{1-\alpha}}, \quad \tilde{b}_k = \frac{\tau^{1-\alpha}}{\Gamma(2-\alpha)} \left[ \left(k + \frac{1}{2}\right)^{1-\alpha} - \left(k - \frac{1}{2}\right)^{1-\alpha} \right], \quad k = 1, 2, \dots, n. \quad (4.3)$$

- In the graded mesh case, a direct calculation using $t_n = (\frac{n}{M})^r T$ gives $\tau = \Delta t_M$

$$\tilde{b}_0 = \frac{T^{1-\alpha}}{\Gamma(2-\alpha)(2Mr)^{1-\alpha}[(n+1)^r - n^r]^\alpha},$$

$$\tilde{b}_k = \frac{T^{1-\alpha}}{\Gamma(2-\alpha)(2Mr)^{1-\alpha}} \frac{[(n+1)^r + n^r - 2(n-k)^r]^{1-\alpha} - [(n+1)^r + n^r - 2(n-k+1)^r]^{1-\alpha}}{(n-k+1)^r - (n-k)^r}, \quad k = 1, 2, \dots, n.$$
The scheme (4.1) is expected to have $2-\alpha$ order convergence, since intuitively the approximation $L_{n+\frac{1}{2}}^{\alpha} \phi$ to $C_0^1 D_{t_{n+\frac{1}{2}}}^{\alpha} \phi$ is of $2-\alpha$ order, and the remaining approximation is of second order. Indeed the approximation order of the operator $L_{n+\frac{1}{2}}^{\alpha}$ can be derived rigorously, as shown in the following lemma.

**Lemma 4.1.** For the graded mesh, $\alpha \in (0,1)$ and $\phi \in C^2[0, t_{n+1}]$, we have

$$|C_0^1 D_{t_{n+\frac{1}{2}}}^{\alpha} \phi - L_{n+\frac{1}{2}}^{\alpha} \phi| = O(\tau^{2-\alpha}).$$

*Proof.* On each small interval $[t_k, t_{k+1}]$ ($0 \le k \le n$), denoting the linear interpolation function of $\phi(t)$ as $\Pi_{1,k}\phi$:

$$\Pi_{1,k}\phi = \frac{t_{k+1} - t}{\Delta t_k}\phi(\mathbf{x}, t_k) + \frac{t - t_k}{\Delta t_k}\phi(\mathbf{x}, t_{k+1}),$$

it follows from the linear interpolation theorem that

$$\phi - \Pi_{1,k}\phi = \frac{\partial_t^2 \phi(\mathbf{x}, \xi_k)}{2}(t - t_k)(t - t_{k+1}), \quad t \in [t_k, t_{k+1}], \xi_k \in (t_k, t_{k+1}), 0 \le k \le n.$$

Let $C_0^1 D_{t_{n+\frac{1}{2}}}^{\alpha} \phi - L_{n+\frac{1}{2}}^{\alpha} \phi = \hat{R}^n + \hat{R}^{n+\frac{1}{2}}$, where

$$\begin{aligned}
\hat{R}^n &= \frac{1}{\Gamma(1-\alpha)} \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} (t_{n+\frac{1}{2}} - s)^{-\alpha} \phi_s(\mathbf{x}, s) ds - \frac{1}{\Gamma(1-\alpha)} \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} (t_{n+\frac{1}{2}} - s)^{-\alpha} \partial_s (\Pi_{1,k}\phi) ds \\
&= \frac{1}{\Gamma(1-\alpha)} \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} (t_{n+\frac{1}{2}} - s)^{-\alpha} \partial_s (\phi(\mathbf{x}, s) - \Pi_{1,k}\phi) ds \\
&= \frac{-\alpha}{\Gamma(1-\alpha)} \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} (t_{n+\frac{1}{2}} - s)^{-\alpha-1} (\phi(\mathbf{x}, s) - \Pi_{1,k}\phi) ds \\
&= \frac{-\alpha}{2\Gamma(1-\alpha)} \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \partial_t^2 \phi(\mathbf{x}, \xi_k) (s-t_k)(s-t_{k+1})(t_{n+\frac{1}{2}}-s)^{-\alpha-1} ds,
\end{aligned}$$

$$\begin{aligned}
\hat{R}^{n+\frac{1}{2}} &= \frac{1}{\Gamma(1-\alpha)} \int_{t_n}^{t_{n+\frac{1}{2}}} (t_{n+\frac{1}{2}} - s)^{-\alpha} \phi_s(\mathbf{x}, s) ds \\
&\quad - \frac{1}{\Gamma(1-\alpha)} \int_{t_n}^{t_{n+\frac{1}{2}}} (t_{n+\frac{1}{2}} - s)^{-\alpha} \frac{\phi(\mathbf{x}, t_{n+1}) - \phi(\mathbf{x}, t_n)}{\Delta t_n} ds \\
&= \frac{1}{\Gamma(1-\alpha)} \int_{t_n}^{t_{n+\frac{1}{2}}} (t_{n+\frac{1}{2}} - s)^{-\alpha} (\phi_t(\mathbf{x}, t_{n+\frac{1}{2}}) - \frac{u^{n+1} - u^n}{\Delta t_n}) ds \\
&\quad + \frac{1}{\Gamma(1-\alpha)} \int_{t_n}^{t_{n+\frac{1}{2}}} (t_{n+\frac{1}{2}} - s)^{-\alpha} \partial_t^2 \phi(\mathbf{x}, \theta_s) (s - t_{n+\frac{1}{2}}) ds, \quad \theta_s \in (s, t_{n+\frac{1}{2}}).
\end{aligned}$$
Noting $\Delta t_{k-1} \le \Delta t_k$ for graded mesh, we have following estimates

$$
\begin{align*}
|\hat{R}^n| &\le \frac{\alpha \max_{0<t<t_{n+1}} |\partial_t^2 \phi(\mathbf{x}, t)|}{2\Gamma(1-\alpha)} \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} (s-t_k)(t_{k+1}-s)(t_{n+\frac{1}{2}}-s)^{-\alpha-1} ds \\
&\le \frac{\alpha \max_{0<t<t_{n+1}} |\partial_t^2 \phi(\mathbf{x}, t)| \Delta t_{n-1}^2}{8\Gamma(1-\alpha)} \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} (t_{n+\frac{1}{2}}-s)^{-\alpha-1} ds \\
&= \frac{\alpha \max_{0<t<t_{n+1}} |\partial_t^2 \phi(\mathbf{x}, t)| \Delta t_{n-1}^2}{8\Gamma(1-\alpha)} \int_0^{t_n} (t_{n+\frac{1}{2}}-s)^{-\alpha-1} ds \\
&= \frac{\max_{0<t<t_{n+1}} |\partial_t^2 \phi(\mathbf{x}, t)| \Delta t_{n-1}^2}{8\Gamma(1-\alpha)} \left( \frac{2^\alpha}{\Delta t_n^\alpha} - \frac{1}{t_{n+\frac{1}{2}}^\alpha} \right) \le \frac{\max_{0<t<t_{n+1}} |\partial_t^2 \phi(\mathbf{x}, t)| \Delta t_n^{2-\alpha}}{2^{3-\alpha}\Gamma(1-\alpha)}
\end{align*}
$$

and

$$
\begin{align*}
|\hat{R}^{n+\frac{1}{2}}| &\le \frac{\max_{0<t<t_{n+1}} |\partial_t^2 \phi(\mathbf{x}, t)|}{\Gamma(1-\alpha)} \int_{t_n}^{t_{n+\frac{1}{2}}} (t_{n+\frac{1}{2}}-s)^{1-\alpha} ds + \mathcal{O}(\Delta t_n^{3-\alpha}) \\
&= \frac{\max_{0<t<t_{n+1}} |\partial_t^2 \phi(\mathbf{x}, t)| \Delta t_n^{2-\alpha}}{2^{2-\alpha}(2-\alpha)\Gamma(1-\alpha)} + \mathcal{O}(\Delta t_n^{3-\alpha}).
\end{align*}
$$

Then we complete the proof. $\square$

It is readily seen that the approximation $L_{n+\frac{1}{2}}^{\alpha} \phi$ to $\mathcal{C}^1 D_{n+\frac{1}{2}}^{\alpha} \phi$ is of $2-\alpha$ order, and the remaining approximation is of second order. Thus the scheme (4.1) is expected to have $2-\alpha$ order convergence.

**4.2. Stability.** We will provide a stability proof for the scheme (4.1) in the uniform mesh case. The case of the graded mesh seems to be much technique, and requires further investigation. The following lemma plays a key role in proving the stability for the uniform mesh.

**Lemma 4.2.** If the mesh is uniform, it holds for any $(u_0, u_1, \dots, u_n)^T \in \mathbb{R}^{n+1}$:

$$
\sum_{k=0}^{n} \sum_{j=0}^{k} \tilde{b}_{k-j} u_j u_k > 0, \qquad (4.4)
$$

where the coefficient set $\{\tilde{b}_k\}_0^n$ is defined in (4.3).

*Proof.* We first define a piecewise constant function $u(t)$ as follows:

$$
u(t) := \begin{cases} u_{\lfloor t/\tau \rfloor}, & 0 \le t < (n+1)\tau, \\ 0, & \text{otherwise,} \end{cases}
$$

where $\lfloor t/\tau \rfloor$ denote the integer part of the real number $t/\tau$. Obviously $u(t)$ is a function in $L^2(0, (n+1)\tau)$. Let's define

$$
\mathcal{A}_{\alpha}^{n}(u, u) := \frac{1}{\Gamma(1-\alpha)} \int_{0}^{(n+1)\tau} \int_{0}^{t} (t-s)^{-\alpha} u(s)u(t) dsdt.
$$
It has been proved in [36] Lemma 2.1 that

$$ \mathcal{A}_{\alpha}^{n}(u, u) \geq 0. \qquad (4.5) $$

Furthermore, a direct calculation shows

$$ \begin{aligned} \mathcal{A}_{\alpha}^{n}(u, u) &= \frac{1}{2\Gamma(1-\alpha)} \int_{0}^{(n+1)\tau} \int_{0}^{(n+1)\tau} |t-s|^{-\alpha} u(s)u(t) dsdt \\ &= \frac{1}{2\Gamma(1-\alpha)} \sum_{k=0}^{n} u_k \sum_{j=0}^{n} u_j \int_{k\tau}^{(k+1)\tau} \int_{j\tau}^{(j+1)\tau} |t-s|^{-\alpha} dsdt \quad (4.6) \\ &:= \frac{\tau}{2} \left[ \hat{b}_0 \sum_{k=0}^{n} u_k^2 + \sum_{k=0}^{n} \sum_{j=0}^{n} \hat{b}_{|k-j|} u_k u_j \right], \end{aligned} $$

where $\hat{b}_0 = \frac{1}{\Gamma(3-\alpha)} \tau^{1-\alpha}$ and

$$ \hat{b}_k = \frac{\tau^{1-\alpha}}{\Gamma(3-\alpha)} \left[ (k+1)^{2-\alpha} - 2k^{2-\alpha} + (k-1)^{2-\alpha} \right], \quad k=1,2,\dots. \qquad (4.7) $$

Let

$$ B := (\tilde{b}_0 - \hat{b}_0) \sum_{k=0}^{n} u_k^2 + \sum_{k=0}^{n} \sum_{j=0}^{n} (\tilde{b}_{|k-j|} - \hat{b}_{|k-j|}) u_j u_k. $$

Then it follows from (4.6) and (4.5)

$$ 2 \sum_{k=0}^{n} \sum_{j=0}^{k} \tilde{b}_{k-j} u_j u_k = \tilde{b}_0 \sum_{k=0}^{n} u_k^2 + \sum_{k=0}^{n} \sum_{j=0}^{n} \tilde{b}_{k-j} u_j u_k = B + \frac{2}{\tau} \mathcal{A}_{\alpha}^{n}(u, u) > B. $$

Therefore (4.4) is true if $B > 0$, which we want to prove below.

Let $c_k = \tilde{b}_{|k|} - \hat{b}_{|k|}$, $k = \pm 1, \pm 2, \dots, n$. Define the function

$$ g(x) := \frac{(x+1)^{2-\alpha} - x^{2-\alpha}}{2-\alpha}, \quad x \ge 0. $$

We deduce from the Taylor expansion and the definitions of $\tilde{b}_k$ and $\hat{b}_k$ in (4.3) and (4.7) that $c_k = c_{-k}$, and, for $k \ge 1$,

$$ \begin{aligned} c_k &= \tilde{b}_{|k|} - \hat{b}_{|k|} = \frac{\tau^{1-\alpha}}{\Gamma(2-\alpha)} \left[ \left(k + \frac{1}{2}\right)^{1-\alpha} - \left(k - \frac{1}{2}\right)^{1-\alpha} - \frac{(k+1)^{2-\alpha} - 2k^{2-\alpha} + (k-1)^{2-\alpha}}{2-\alpha} \right] \\ &= \frac{\tau^{1-\alpha}}{\Gamma(2-\alpha)} \left[ g'(k - \frac{1}{2}) - g(k) + g(k-1) \right] \\ &= \frac{\tau^{1-\alpha}}{\Gamma(2-\alpha)} \left[ g'(k - \frac{1}{2}) - \left( g(k - \frac{1}{2}) + \frac{g'(k - \frac{1}{2})}{2} + \frac{g''(\zeta_1)}{8} \right) + \left( g(k - \frac{1}{2}) - \frac{g'(k - \frac{1}{2})}{2} + \frac{g''(\zeta_2)}{8} \right) \right] \\ &= \frac{\tau^{1-\alpha}}{\Gamma(2-\alpha)} \frac{g''(\zeta_2) - g''(\zeta_1)}{8}, \end{aligned} $$
where $\zeta_1 \in (k - \frac{1}{2}, k)$ and $\zeta_2 \in (k - 1, k - \frac{1}{2})$. Since $g'''(x) > 0$ for all $x > 0$, then it holds $g''(\zeta_2) < g''(\zeta_1)$. Therefore, we have $c_k = c_{-k} < 0$ for $k = 1, 2, \dots, n$, and

$$ 
\begin{aligned}
\sum_{k=1}^{n} |c_k| &= \sum_{k=1}^{n} \frac{\tau^{1-\alpha}}{\Gamma(2-\alpha)} \left[ \frac{(k+1)^{2-\alpha} - 2k^{2-\alpha} + (k-1)^{2-\alpha}}{2-\alpha} - \left( \left(k+\frac{1}{2}\right)^{1-\alpha} - \left(k-\frac{1}{2}\right)^{1-\alpha} \right) \right] \\
&= \frac{\tau^{1-\alpha}}{\Gamma(2-\alpha)} \left[ \frac{(n+1)^{2-\alpha} - 2n^{2-\alpha}}{2-\alpha} - \frac{1}{2-\alpha} - \left(n+\frac{1}{2}\right)^{1-\alpha} + \frac{1}{2^{1-\alpha}} \right].
\end{aligned} 
\quad (4.8) $$

Making use of Taylor expansion with integral remainder and differential mean value theorem, we obtain

$$ 
\begin{aligned}
& \frac{(n+1)^{2-\alpha} - 2n^{2-\alpha}}{2-\alpha} - (n+\frac{1}{2})^{1-\alpha} \\
&= \frac{1}{2-\alpha} \left\{ \left[ (n+\frac{1}{2})^{2-\alpha} + \frac{2-\alpha}{2}(n+\frac{1}{2})^{1-\alpha} + (2-\alpha)(1-\alpha) \int_{n+1/2}^{n+1} (n+1-s)s^{-\alpha}ds \right] \right. \\
& \qquad \left. - \left[ (n+\frac{1}{2})^{2-\alpha} - \frac{2-\alpha}{2}(n+\frac{1}{2})^{1-\alpha} + (2-\alpha)(1-\alpha) \int_{n+1/2}^{n} (n-s)s^{-\alpha}ds \right] \right\} - (n+\frac{1}{2})^{1-\alpha} \\
&= (1-\alpha) \left[ \int_{n+1/2}^{n+1} (n+1-s)s^{-\alpha}ds - \int_{n+1/2}^{n} (n-s)s^{-\alpha}ds \right] \\
&= (1-\alpha) \int_{0}^{1/2} s[(n+1-s)^{-\alpha} - (n+s)^{-\alpha}]ds \leq -\alpha(1-\alpha)(n+1)^{-\alpha-1} \int_{0}^{1/2} s(1-2s)ds \\
&= -\frac{\alpha(1-\alpha)}{24}(n+1)^{-\alpha-1}.
\end{aligned} 
\quad (4.9) $$

Combining (4.8) and (4.9) gives

$$ \sum_{k=1}^{n} |c_k| \leq \frac{\tau^{1-\alpha}}{\Gamma(2-\alpha)} \left[ \frac{1}{2^{1-\alpha}} - \frac{1}{2-\alpha} - \frac{\alpha(1-\alpha)}{24} (n+1)^{-\alpha-1} \right]. $$

Let $c_0 = \frac{\tau^{1-\alpha}}{\Gamma(2-\alpha)}[2^\alpha - \frac{2}{2-\alpha} - \frac{\alpha(1-\alpha)}{12}(n+1)^{-\alpha-1}]$. We construct the matrix $C := \{c_{k-j}\}_{k,j=0}^n$. Then

$$ \sum_{j=0, j \neq k}^{n} |C_{k,j}| = \sum_{j=0, j \neq k}^{n} |c_{k-j}| < 2 \sum_{l=1}^{n} |c_l| \leq c_0, \text{ for any } k = 0, 1, 2, \dots, n. $$

This means $C$ is a symmetric positive definite matrix. Thus $\sum_{k=0}^{n} \sum_{j=0}^{n} C_{k,j}u_j u_k > 0$, for any $(u_0, u_1, \dots, u_n)^T \in \mathbb{R}^{n+1}$.

It then follows from the definitions of $B$ and $C_{k,j}$ that

$$ 
\begin{aligned}
B &= \sum_{k=0}^{n} (2\tilde{b}_0 - \hat{b}_0 - c_0)u_k^2 + \sum_{k=0}^{n} \sum_{j=0}^{n} C_{k,j}u_j u_k \\
&= \frac{\tau^{1-\alpha}}{\Gamma(2-\alpha)} \frac{\alpha(1-\alpha)}{12} (n+1)^{-\alpha-1} \sum_{k=0}^{n} u_k^2 + \sum_{k=0}^{n} \sum_{j=0}^{n} C_{k,j}u_j u_k > 0.
\end{aligned} $$

This completes the proof.

<figure><img src="image_4.png" /></figure>

Now we are in a position to establish the unconditional stability of the L1-CN scheme (4.1).
**Theorem 4.1.** In the case of uniform mesh, the L1-CN scheme (4.1) is unconditionally stable in the sense that a discrete energy is always bounded during the time stepping. More precisely, it holds

$$ \tilde{E}_{\varepsilon, \theta}^{n} \leq \tilde{E}_{\varepsilon, \theta}^{0}, \quad n = 1, 2, \dots, \qquad (4.10) $$

where the discrete energy $E_{\varepsilon, \theta}^n$ is defined as:

$$ \tilde{E}_{\varepsilon, \theta}^{n} := \frac{\varepsilon^2 - \theta}{2} \|\nabla \phi^{n}\|_{0}^{2} + \frac{\theta}{4} \|\nabla(\phi^{n} - \phi^{n-1})\|_{0}^{2} + |R^{n}|^{2}. \qquad (4.11) $$

*Proof.* By taking the inner products of (4.1a) and (4.1b) with $\frac{\phi^{n+1}-\phi^n}{\tau}$ and $R^{n+1}+R^n$ respectively, we have

$$ \left(L_{n+\frac{1}{2}}^{\alpha} \phi, \frac{\phi^{n+1}-\phi^n}{\tau}\right) + \frac{\varepsilon^2 - \theta}{2\tau} \left(\|\nabla \phi^{n+1}\|^2 - \|\nabla \phi^n\|^2\right) + \frac{\theta}{2\tau} \left(\nabla(\phi^{n+1} - 2\phi^n + \phi^{n-1}), \nabla(\phi^{n+1} - \phi^n)\right) \\ + \frac{R^{n+1} + R^n}{2\sqrt{E_{\theta}^{n+\frac{1}{2}} + C_0}} \left(-\theta \Delta \phi^{n+\frac{1}{2}} + F'(\phi^{n+\frac{1}{2}}), \frac{\phi^{n+1}-\phi^n}{\tau}\right) = 0, $$

$$ \frac{|R^{n+1}|^2 - |R^n|^2}{\tau} = \frac{R^{n+1} + R^n}{2\sqrt{E_{\theta}^{n+\frac{1}{2}} + C_0}} \left( -\theta \Delta \phi^{n+\frac{1}{2}} + F'(\phi^{n+\frac{1}{2}}), \frac{\phi^{n+1} - \phi^n}{\tau} \right). $$

Then we apply the identities

$$ 2(a^{n+1} - a^n)(a^{n+1} - 2a^n + a^{n-1}) = (a^{n+1} - a^n)^2 - (a^n - a^{n-1})^2 + (a^{n+1} - 2a^n + a^{n-1})^2. $$

to the above equations and drop some non-essential positive terms to obtain

$$ \frac{\tilde{E}_{\varepsilon, \theta}^{n+1} - \tilde{E}_{\varepsilon, \theta}^{n}}{\tau} \leq - \left( L_{n+\frac{1}{2}}^{\alpha} \phi, \frac{\phi^{n+1} - \phi^n}{\tau} \right), \quad n = 0, 1, 2, \dots. $$

Summing up the above inequality from 0 to n yields

$$ \frac{\tilde{E}_{\varepsilon, \theta}^{n+1} - \tilde{E}_{\varepsilon, \theta}^{0}}{\tau} \leq - \sum_{k=0}^{n} \left( L_{k+\frac{1}{2}}^{\alpha} \phi, \frac{\phi^{k+1} - \phi^{k}}{\tau} \right), \quad n = 0, 1, 2, \dots. $$

We know from (4.2) that the right hand side is equal to $-\sum_{k=0}^n \sum_{j=0}^k \tilde{b}_{k-j} \frac{\phi^{j+1}-\phi^j}{\tau} \frac{\phi^{k+1}-\phi^k}{\tau}$, which is non-positive according to Lemma 4.2. Thus $\tilde{E}_{\varepsilon, \theta}^n \leq \tilde{E}_{\varepsilon, \theta}^0, n = 1, 2, \dots$. This proves the theorem. $\square$

Since the L1-CN scheme (4.1) also enjoys the unconditional energy stability similar to (3.5) for the uniform mesh, we can derive $H^2$ bound for $\phi^n$ similar to Theorem 3.2.

**Theorem 4.2.** Assume $\phi_0 \in H^3(\Omega)$ and $F \in C^3(\mathbb{R})$ satisfying (2.3). Then for the uniform mesh and any $(n+1) < T/\tau$, the scheme (4.1) satisfies

$$ \sum_{k=0}^{n} \tilde{b}_{n-k} \|\Delta \phi^{k+1}\|_{0}^{2} + \frac{\tau \varepsilon^2}{2} \|\nabla \Delta \phi^{k+1}\|_{0}^{2} \leq c(M_0, \varepsilon)T + \frac{T^{1-\alpha}}{\Gamma(2-\alpha)} \|\Delta \phi_0\|_{0}^{2} + \frac{\tau \varepsilon^2}{2} \|\nabla \Delta \phi_0\|_{0}^{2}, $$

where $M_0$ is a positive constant depending only on $\Omega$ and $\phi_0$.

The L1-CN scheme (4.1) can be efficiently implemented by following the lines similar to the algorithm (i)-(ii) in the previous section.
## 5. NUMERICAL RESULTS

In this section, we present some numerical examples to demonstrate the efficiency of the proposed schemes in terms of stability and accuracy. In all tests that follow, we always set $\theta = \varepsilon^2$ in the schemes. The spatial discretization is the Fourier method or Legendre spectral method using numerical quadratures. In order to test the accuracy, the error is measured by the maximum norm, i.e., $\max_{1 \le n \le M} \|\phi^n - \phi(t_n)\|_{\infty}$ or $\|\phi^M - \phi(T)\|_{\infty}$ depending on whether or not the exact solution is exactly known. It is worth to mention that a fast evaluation technique based on the so-called sum-of-exponentials approach to the time fractional derivative is used to accelerate the calculation (also to reduce the storage); see, e.g., [23, 37].

### 5.1. Convergence order test.

**Example 5.1.** Consider the following time fractional Allen-Cahn equation with the periodic boundary condition:

$$C_0^\alpha D_t^\alpha \phi - \varepsilon^2 \Delta \phi - \phi(1 - \phi^2) = s(\mathbf{x}, t), \quad (\mathbf{x}, t) \in (0, 2\pi)^2 \times (0, T], \quad (5.1)$$

where $s(\mathbf{x}, t)$ is a fabricated source term chosen such that the exact solution is

$$\phi(\mathbf{x}, t) = 0.2t^5 \sin(x) \cos(y).$$

The spatial discretization used in the calculation is the Fourier method with 128 × 128 modes. It has been checked that this Fourier mode number is large enough so that the spatial discretization error is negligible compared to the temporal discretization. We present in Figure 1 the error at $T = 1$ as functions of the time step sizes in log-log scale. It is clearly observed that the first order scheme (3.12) and L1-CN scheme (4.1) achieve the expected convergence rate for all tested $\alpha$, i.e., first order and $2 - \alpha$ order, respectively. Note that there is no numerical instability observed during the calculation for all the time step sizes we tested.

**Example 5.2.** Consider the same equation as in Example 5.1, but the Neumann boundary condition, with the exact solution

$$\phi(\mathbf{x}, t) = 0.2(t^\mu + 1) \cos(\pi x) \cos(\pi y), \quad (\mathbf{x}, t) \in (-1, 1)^2 \times (0, T],$$

which has limited regularity at the initial time $t = 0$.

For this Neumann problem, the space variable is discretized by the Legendre Galerkin spectral method using polynomials of degree 32 in each spatial direction. The purpose of this test is to not only verify the accuracy of the schemes, but also investigate how the regularity of the solution affects the accuracy. In particular, we are interested to study the impact of the graded mesh parameter $r$ on the convergence rate. The calculation is performed by using the L1-CN scheme with $M = 2^k, k = 4, 5, \dots, 13$. In Figure 2, we plot the $L^\infty$ errors in log-log scale with respect to the maximum time step sizes $\tau$ for each fixed $M$. It is shown in the figure that the numerical solution achieves the convergence rate $\min\{\mu r, 2 - \alpha\}$. It is worth to point out that the error curves can be misleading if one only looks at the error behaviour for relatively larger time step sizes. For example, the numerical results given in Figure 2 (d) and (f) for the case $\mu r > 2 - \alpha$ seem to make us expect a higher order accuracy. More precisely, the error curve in Figure 2 (d) indicates a convergence rate like $\mu r$ in the range of larger $\tau$. However when $\tau$ decreases we clearly
<figure><img src="image_5.png" /><figcaption>(a) first order scheme for $\alpha = 0.1$ and $0.5$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) L1-CN scheme for $\alpha = 0.6$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) L1-CN scheme for $\alpha = 0.9$</figcaption></figure>

FIGURE 1. Example 5.1: Error decay at $T = 1$ versus the time step sizes for the first order scheme and L1-CN scheme with $C_0 = 0$.

observe the convergence order $2-\alpha$, which is what we can expect from a truncation error analysis for the L1 scheme.

**Example 5.3.** Consider the time fractional Allen-Cahn Neumann problem in the domain $(-1, 1) \times (-1, 1)$ with the double well potential $F(\phi) := \frac{(\phi^2 - 1)^2}{4}$ and the initial condition $\phi(x, 0) = \cos(4\pi x) \cos(4\pi y)$.

Again, we use the Legendre Galerkin method for the spatial discretization with high enough mode number to avoid possible spatial error contamination. The error behavior is investigated by comparing the computed solutions to the one obtained by using the L1-CN scheme (4.1) with $r = 3$, $M = 5 \times 10^4$ and Legendre polynomial space of degree $32 \times 32$, considered as the “exact” solution. In Figure 3, we plot the $L^\infty$ errors at $T = 1$ in log-log scale with respect to the time discretization parameter $\tau$. Once again the observed error behavior in Figure 3 demonstrates
<figure><img src="image_8.png" /><figcaption>FIGURE 2. (Example 5.2) Error history for the L1-CN scheme with $C_0 = 0$ for different values of the mesh parameter $r$: the left side figures correspond to $\alpha = 0.5$ and $\mu = 0.4$; the right side are for $\alpha = 0.9$ and $\mu = 0.9$.</figcaption></figure>
that the proposed L1-CN scheme (4.1) attains the convergence rate $\min\{\alpha r, 2 - \alpha\}$ for all tested $\alpha$ and $r$.

The stability of the schemes is investigated using this example through a long time calculation up to $T = 100$. In view of the possible singularity feature of the solution, we split the interval $[0, T]$ into two parts: $[0, 1]$ and $(1, T]$. We first compute the solution by using the graded mesh with $r = \frac{2-\alpha}{\alpha}$ in the first subinterval $[0, 1]$, then use the uniform mesh with the time step size $\Delta t$ in the second subinterval $(1, T]$. The computed discrete energies, both $E(\phi^n)$ and $\tilde{E}_{\varepsilon, \theta}^n$, by using the L1-CN scheme are presented in Figure 4 (b), (d), and (f), showing dissipative feature during the running time for all tested values of $\alpha, M$ and $\Delta t$. This confirms the unconditional stability of the proposed scheme. We would like to emphasize that while the modified energy $\tilde{E}_{\varepsilon, \theta}^n$ keeps dissipative all the time, the original energy $E(\phi^n)$ may exhibit oscillation at some time instants, especially when large time step size is used. This has been a well known fact in the SAV approach; see, e.g., [20] and the reference therein. Thus it is suggested that, in order to avoid undesirable oscillation of the original energy in real applications, the time step size should not be taken too large, even though such oscillations do not necessarily make the calculation unstable.

### 5.2. Order sensibility of a benchmark problem.
Now we consider an interface moving problem governed by the following fractional Allen-Cahn equation in the rectangular domain $\Omega := (-32, 32)^2$

$$ egin{cases} 
\frac{\partial}{\partial t} D_t^\alpha \phi - \Delta \phi - \phi(1 - \phi^2) = 0, \\ 
\frac{\partial \phi}{\partial n} \bigg|_{\partial \Omega} = 0, \phi(x, 0) = \phi_0(x). 
\end{cases} 
\quad (5.2) $$

The initial state of the interface is a circle of the radius $R_0$. In the classical case, i.e., $\alpha = 1$, this interface moving problem has been frequently served as benchmark problem to validate numerical methods. In this case it has been known [11] that such a circle interface is unstable, and will shrink and disappear due to the driving force. In fact, it was shown in [38] that the radius $R(t)$ of the circle at the given time $t$ evolves as

$$ R(t) = \sqrt{R_0^2 - 2t}. \quad (5.3) $$

The purpose of this example is twofold: demonstration of the efficiency of the proposed schemes, and sensibility investigation of how the fractional order $\alpha$ effects the interface evolution. To this end we consider the circle interface of the radius $R_0 = 8$. We first transform the computational domain into the reference domain $(-1, 1)^2$, then the equation to be solved reads:

$$ \frac{\partial}{\partial t} D_t^\alpha \phi - \varepsilon^2 \Delta \phi - \phi(1 - \phi^2) = 0, \quad (5.4) $$

where $\varepsilon = 0.0313$. The initial condition becomes

$$ \phi_0(x) = \begin{cases} 1, & |x|^2 < \left(\frac{8}{32}\right)^2, \\ -1, & |x|^2 \ge \left(\frac{8}{32}\right)^2. \end{cases} $$

We numerically solve this problem by using the L1-CN time scheme (4.1) based on the graded mesh with $r = \frac{2-\alpha}{\alpha}$ in the subinterval $[0, 1]$ and the uniform mesh with the time step size $\Delta t$ in the subinterval $(1, T]$. The Legendre spectral method for the spatial discretization make use of polynomials of degree 128 in each direction. In Figure 5 we investigate the impact of the
<figure><img src="image_9.png" /><figcaption>FIGURE 3. (Example 5.3) Effect of the mesh parameter on the convergence rate of the L1-CN scheme: $\alpha = 0.7$ (left) and $\alpha = 0.9$ (right).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIGURE 4. (Example 5.3) Dissipation feature of the computed energy $E(\phi^n)$ and modified energy $\tilde{E}_{\varepsilon,\theta}^n$ with different time step sizes using the L1-CN scheme.</figcaption></figure>
fractional order $\alpha$ on the interface evolution. First, the result for $\alpha = 1$ indicates that the circle interface almost disappears at $T = 32$, which is in a perfect agreement with the prediction (5.3). This partially demonstrates the accuracy of the scheme. The second observation is that the interface shirking dynamics exists for all $\alpha$, but the shirking speed slows down significantly when the fractional order decreases. This is not really surprising since the fractional derivative represents some kind of memory effect. More details about the evolution of the circle radius are presented in Figure 6(a), which confirms the result reported in [15]. That is, the circle shrinking rate seems to exhibit a power-law scaling, although the exact rate remains unknown. Also shown in Figure 6(b) is the evolution of the total free energy, from which we see that the free energy is always dissipative for all tested $\alpha$. This further increases the credibility of the simulation results.

<figure><img src="image_11.png" /><figcaption>FIGURE 5. Snapshots of the interface evolution simulated by using the L1-CN scheme for $\alpha = 1, 0.9$, and 0.4.</figcaption></figure>

5.3. **Coarsening dynamics.** Finally we consider the application of the proposed schemes in investigating coarsening dynamics of two-phase problems, governed by the time fractional Allen-Cahn equation (1.4), and starting with a random initial data ranging from -0.05 to 0.05.
<figure><img src="image_12.png" /><figcaption>FIGURE 6. The evolution of computed $R^2$ and total free energy for time fractional Allen-Cahn equation (5.4) with $\alpha = 1, 0.9$ and $0.4$.</figcaption></figure>

We set the computational domain $\Omega$ to be $(-1, 1)^2$, subject to the Neumann boundary condition, $\varepsilon^2 = 0.001$. The simulation is performed by using L1-CN scheme (4.1) for the time discretization in the graded mesh with $r = \frac{2-\alpha}{\alpha}$ in the first subinterval $[0, 1]$ and the uniform mesh with $\Delta t = 0.01$ in the second subinterval $(1, T]$. The spatial discretization uses the Legendre spectral method with $128 \times 128$ basis functions. In Figure 7 we compare the simulation results for different time fractional order $\alpha$. As we can observe from this figure, there is no distinguishable difference at early time during the phase separation period (say, before $t = 5$) for all tested $\alpha$. After that the solutions apparently start to deviate, the dynamics described by the fractional Allen-Cahn equation enters into long time scale phase coarsening. This is further verified by the energy evolution in Figure 8. Numerically, we can observe that the equations with different fractional orders have similar dynamics for the phase separation but the coarsening dynamics becomes slower for smaller fractional orders.
<figure><img src="image_13.png" /><figcaption>FIGURE 7. Snapshots of the simulated phase field evolution by the fractional Allen-Cahn equation with a random initial data for $\alpha = 1, 0.9, 0.7, 0.5$ with $M = 100, \Delta t = 0.01$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIGURE 8. Decay history of the computed free energy by the fractional Allen-Cahn equation for $\alpha = 1, 0.9, 0.7, 0.5$.</figcaption></figure>

## 6. CONCLUDING REMARKS

We have proposed two schemes for the $\alpha$ order time fractional Allen-Cahn equation: a first order scheme and a $2 - \alpha$ order scheme, both constructed on the uniform mesh and graded mesh. Mathematically, such fractional gradient flows can be regarded as a variation of the phase field models, used to describe the memory effect of some materials, which have attracted attentions in recent few years. Our proposed schemes are based on the L1 discretization for the time fractional derivative and extended SAV approach to deal with the nonlinear energy potential. The stability property of the schemes was rigorously established, while the accuracy was carefully examined through a series of numerical tests. Precisely, our stability analysis showed that the schemes constructed on the uniform mesh are unconditionally stable in the sense that the associated discrete energy is always decreasing. As far as we know this is the first proof for the unconditional stability of a $2 - \alpha$ order scheme. The numerical experiments carried out in the paper have demonstrated that the proposed schemes have desirable accuracy. In particular, the initial singularity of the solution was effectively resolved by using the graded mesh. From the coarsening simulation, we can conclude that the fractional Allen-Cahn equation have similar dynamics for the phase separation but the coarsening dynamics slows down with decreasing fractional order. Numerically, it was observed that the time to the steady state becomes longer for smaller fractional order. Notice that this phenomena has already been observed in some other existing works, but theoretical justification remains open. It is also worth to mention that the proposed schemes can be directly
applied to some other gradient flows, such as the Cahn-Hilliard equation and molecular beam epitaxial growth models, and the stability results can be derived similarly.

To summarize, the proposed SAV-based schemes enjoy a few good properties: unconditional stability; only require solving linear decoupled systems with constant coefficients at each time step; only require that the free energy is bounded from below. However there exist several theoretical issues such as rigorous proof of the maximum principle of the numerical solution, error estimation, etc. The difficulty is not only due to the presence of the fractional derivative, but also comes from the SAV approach. In SAV-type schemes, the nonlinear term is treated explicitly, and an auxiliary variable is introduced in front of the nonlinear term. This makes the numerical analysis harder than for other approaches. Possible answer to these issues certainly needs more effort.

## REFERENCES

[1] M. Ainsworth and Z. Mao. Analysis and approximation of a fractional Cahn-Hilliard equation. *SIAM J. Numer. Anal.*, 55(4):1689-1718, 2017.

[2] M. Ainsworth and Z. Mao. Well-posedness of the Cahn-Hilliard equation with fractional free energy and its Fourier Galerkin approximation. *Chaos Soliton. Fract.*, 102:264-273, 2017.

[3] G. Akagi, G. Schimperna, and A. Segatti. Fractional Cahn-Hilliard, Allen-Cahn and porous medium equations. *J. Differ. Equations*, 261(6):2935-2985, 2016.

[4] A. A. Alikhanov. A priori estimates for solutions of boundary value problems for fractional-order equations. *Differ. Equ.*., 46(5):660-666, 2010.

[5] A.A. Alikhanov. A new difference scheme for the time fractional diffusion equation. *J. Comput. Phys.*, 280(C):424-438, 2015.

[6] S. M. Allen and J. W. Cahn. A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. *Acta Metall. Mater.*, 27(6):1085-1095, 1979.

[7] D. M. Anderson, G. B. Mcfadden, and A. A. Wheeler. Diffuse-interface methods in fluid mechanics. *Annu. Rev. Fluid Mech.*, 30:139-165, 1997.

[8] D. Baffet and J.S. Hesthaven. High-order accurate adaptive kernel compression time-stepping schemes for fractional differential equations. *J. Sci. Comput.*, 72(3):1169-1195, 2017.

[9] D. Baffet and J.S. Hesthaven. A kernel compression scheme for fractional differential equations. *SIAM J. Numer. Anal.*, 55(2):496-520, 2017.

[10] J. W. Cahn and J. E. Hilliard. Free Energy of a Nonuniform System. I. Interfacial Free Energy. *J. Chem. Phys.*, 28:258-267, 1958.

[11] L. Chen and J. Shen. Applications of semi-implicit Fourier-spectral method to phase field equations. *Comput. Phys. Commun.*, 108(2-3):147-158, 1998.

[12] L. Chen, J. Zhao, W. Cao, H. Wang, and J. Zhang. An accurate and efficient algorithm for the time-fractional molecular beam epitaxy model with slope selection. *Comput. Phys. Commun.*, 2018.

[13] M. Doi and S. F. Edwards. *The theory of polymer dynamics*, volume 73. Oxford University Press, 1988.

[14] Q. Du, L. Ju, X. Li, and Z. Qiao. Stabilized linear semi-implicit schemes for the nonlocal Cahn-Hilliard equation. *J. Comput. Phys.*, 363:39-54, 2018.

[15] Q. Du, J. Yang, and Z. Zhou. Time-fractional Allen-Cahn equations: analysis and numerical methods. *arXiv:1906.06584v1*, pages 1-24, 2019.

[16] K. R. Elder, M. Katakowski, M. Haataja, and M. Grant. Modeling elasticity in crystal growth. *Phys. Rev. Lett.*, 88(24):245701, 2002.

[17] G. Gao, Z. Sun, and H. Zhang. A new fractional numerical differentiation formula to approximate the caputo fractional derivative and its applications. *J. Comput. Phys.*, 259(2):33-50, 2014.

[18] G. Gao, Z. Sun, and Y. Zhang. A finite difference scheme for fractional sub-diffusion equations on an unbounded domain using artificial boundary conditions. *J. Comput. Phys.*, 231(7):2865-2879, 2012.

[19] M. E. Gurtin and D. Polignone. Two-phase binary fluids and immiscible fluids described by an order parameter. *Math. Models Methods Appl. Sci.*, 6:815-831, 1996.
[20] D. Hou, M. Azaiez, and C. Xu. A variant of scalar auxiliary variable approaches for gradient flows. *J. Comput. Phys.*, 395:307–332, 2019. https://doi.org/10.1016/j.jcp.2019.05.037.

[21] B. Ji, H. Liao, Y. Gong, and L. Zhang. Adaptive linear second-order energy stable schemes for time-fractional Allen-Cahn equation with volume constraint. *arXiv:1909.10936*, pages 1–21, 2019.

[22] S. Jiang, J. Zhang, Q. Zhang, and Z. Zhang. Fast evaluation of the caputo fractional derivative and its applications to fractional diffusion equations. *Commun. Comput. Phys.*, 21(3):650–678, 2017.

[23] S. Jiang, J. Zhang, Q. Zhang, and Z. Zhang. Fast evaluation of the Caputo fractional derivative and its applications to fractional diffusion equations. *Commun. Comput. Phys.*, 21(3):650–678, 2017.

[24] F. M. Leslie. Theory of flow phenomena in liquid crystals. *Adv. Liq. Cryst.*, 4:1–81, 1979.

[25] Z. Li, H. Wang, and D. Yang. A space-time fractional phase-field model with tunable sharpness and decay behavior and its efficient numerical simulation. *J. Comput. Phys.*, 347:20–38, 2017.

[26] H. Liao, T. Tang, and T. Zhou. A second-order and nonuniform time-stepping maximum-principle preserving scheme for time-fractional Allen-Cahn equations. *arXiv:1909.10216*, pages 1–22, 2019.

[27] Y. Lin and C. Xu. Finite difference/spectral approximations for the time-fractional diffusion equation. *J. Comput. Phys.*, 225(2):1533–1552, 2007.

[28] H. Liu, A. Cheng, H. Wang, and J. Zhao. Time-fractional Allen-Cahn and Cahn-Hilliard phase-field models and their numerical investigation. *Comput. Math. Appl.*, 76:1876–1896, 2018.

[29] C. Lv and C. Xu. Improved error estimates of a finite difference/spectral method for time-fractional diffusion equations. *Int. J. Numer. Anal. Mod.*, 12(2):384–400, 2015.

[30] C. Lv and C. Xu. Error analysis of a high order method for time-fractional diffusion equation. *SIAM J. Sci. Comput.*, 38(5):A2699–A2724, 2016.

[31] J. Shen and J. Xu. Convergence and error analysis for the scalar auxiliary variable(SAV) schemes to gradient flows. *SIAM J. Numer. Anal.*, 56(5):2895–2912, 2018.

[32] J. Shen, J. Xu, and J. Yang. The scalar auxiliary variable (SAV) approach for gradient flows. *J. Comput. Phys.*, 353:407–416, 2018.

[33] J. Shen, J. Xu, and J. Yang. A new class of efficient and robust energy stable schemes for gradient flows. *SIAM Rev.*, 61(3):474–506, 2019.

[34] F. Song, C. Xu, and G. Em Karniadakis. A fractional phase-field model for two-phase flows with tunable sharpness: Algorithms and simulations. *Comput. Methods Appl. Mech. Engrg.*, 305:376–404, 2016.

[35] Z. Sun and X. Wu. A fully discrete difference scheme for a diffusion-wave system. *Appl. Numer. Math.*, 56(2):193–209, 2006.

[36] T. Tang, H. Yu, and T. Zhou. On energy dissipation theory and numerical stability for time-fractional phase field equations. *SIAM J. Sci. Comput.*, 41(6):A3757–A3778, 2019.

[37] Y. Yan, Z. Sun, and J. Zhang. Fast evaluation of the caputo fractional derivative and its applications to fractional diffusion equations: A second-order scheme. *Commun. Comput. Phys.*, 22(4):1028–1048, 2017.

[38] X. Yang. Linear, first and second-order, unconditionally energy stable numerical schemes for the phase field model of homopolymer blends. *J. Comput. Phys.*, 327:294–316, 2016.

[39] P. Yue, J. Feng, C. Liu, and J. Shen. A diffuse-interface method for simulating two-phase flows of complex fluids. *J. Fluid Mech.*, 515:293–317, 2004.

[40] F. Zeng, I. Turner, and K. Burrage. A stable fast time-stepping method for fractional integral and derivative operators. *J. Sci. Comput.*, 77(1):283–307, 2018.

[41] Q. Zhang, J. Zhang, S. Jiang, and Z. Zhang. Numerical solution to a linearized time fractional KdV equation on unbounded domains. *Math. Comp.*, 87(310), 2018.

[42] J. Zhao, L. Chen, and H. Wang. On power law scaling dynamics for time-fractional phase field models during coarsening. *Commun. Nonlinear Sci. Numer. Simulat.*, 70:257–270, 2019.

[43] J. Zhao, Q. Wang, and X. Yang. Numerical approximations for a phase field dendritic crystal growth model based on the invariant energy quadratization approach. *Int. J. Numer. Meth. Eng.*, 110(3):279–300, 2017.

[44] G. Zhen, J. Lowengrub, C. Wang, and S. Wise. Second order convex splitting schemes for periodic nonlocal Cahn-Hilliard and Allen-Cahn equations. *J. Comput. Phys.*, 277:48–71, 2014.