Adaptive second-order Crank-Nicolson time-stepping schemes for time fractional molecular beam epitaxial growth models

TL;DR

This paper introduces adaptive second-order Crank-Nicolson schemes using the SAV approach for efficient and stable simulation of time-fractional molecular beam epitaxial growth models, especially on nonuniform meshes.

Contribution

It develops a novel second-order approximation for Caputo's fractional derivative compatible with nonuniform meshes and demonstrates unconditional energy stability of the schemes.

Findings

Schemes are unconditionally energy stable.

Effective in multiscale time simulations.

Accurately captures initial singularity and fast dynamics.

Abstract

Adaptive second-order Crank-Nicolson time-stepping methods using the recent scalar auxiliary variable (SAV) approach are developed for the time-fractional Molecular Beam Epitaxial models with Caputo's derivative. Based on the piecewise linear interpolation, the Caputo's fractional derivative is approximated by a novel second-order formula, which is naturally suitable for a general class of nonuniform meshes and essentially preserves the positive semi-definite property of integral kernel. The resulting Crank-Nicolson SAV time-stepping schemes are unconditional energy stable on nonuniform time meshes, and are computationally efficient in multiscale time simulations when combined with adaptive time steps, such as are appropriate for accurately resolving the intrinsically initial singularity of solution and for efficiently capturing fast dynamics away from the initial time. Numerical examples are presented to show the effectiveness of our methods.

Equations

$${\partial }_t\varphi =-M\mu ,$$

$$\frac{\mathrm{d}E}{\mathrm{d}t}=\left(\frac{\delta E}{\delta \varphi },\frac{\partial \varphi }{\partial t}\right)=-\frac{1}{M}{\Vert {\partial }_t\varphi \Vert }^2⩽0,$$

$${\partial }_t^{\alpha }\varphi =-M\mu ,$$

$$({\partial }_t^{\alpha }v)(t):=({\mathcal{I}}_t^{1-\alpha }{v}^{\prime })(t)={\int }_0^t{\omega }_{1-\alpha }(t-s)v^{\prime }(s)\mathrm{d}s,0<\alpha <1,$$

$$({\mathcal{I}}_t^{\beta }v)(t):={\int }_0^t{\omega }_{\beta }(t-s)v(s)\mathrm{d}s,\text{where}{\omega }_{\beta }(t):=t^{\beta -1}/\Gamma (\beta ).$$

$$\displaystyle E\left[\phi(T)\right]-E\left[\phi(0)\right]=-\frac{1}{M}\int_{\Omega}\int_{0}^{T}\big{(}\partial_{t}^{\alpha}\phi\big{)}\partial_{t}\phi\,\mathrm{d}{t}\,\mathrm{d}{\mathbf{x}}\leqslant{0}.$$

$${\Pi }_{1}v:={\Pi }_{1,k}v\text{so that}({\Pi }_{1}v{)}^{\prime }(t)={\partial }_{\tau }{v}^{k-\frac{1}{2}}\text{for tk−1<t⩽tk and k≥1}.$$

$$\log |u_t(\mathbf{x},t)|\approx (\alpha -1)\log (t)+C(\mathbf{x})\text{as t→0.}$$

$$E[\varphi ]={\int }_{\Omega }\left(\frac{{\epsilon }^2}{2}{\left|\Delta \varphi \right|}^2+F\left(\nabla \varphi \right)\right)\mathrm{d}\mathbf{x}.$$

$${\partial }_t^{\alpha }\varphi =-M\frac{\delta E\left(\varphi \right)}{\delta \varphi }=-M\left({\epsilon }^2{\Delta }^2\varphi +f\left(\nabla \varphi \right)\right),$$

$$\displaystyle f\left(\nabla\phi\right)=-\nabla\cdot\big{(}\big{(}|\nabla\phi|^{2}-1\big{)}\nabla\phi\big{)},$$

$$\displaystyle f\left(\nabla\phi\right)=\nabla\cdot\big{(}\frac{1}{1+|\nabla\phi|^{2}}\nabla\phi\big{)}.$$

$$u(t)=\sqrt{{\int }_{\Omega }\frac{1}{4}{\left({\left|\nabla \varphi \right|}^2-1-\beta \right)}^2\mathrm{d}\mathbf{x}+C_0},$$

$$\displaystyle E\left[\phi,u\right]=\int_{\Omega}\left(\frac{\varepsilon^{2}}{2}\left|\Delta\phi\right|^{2}+\frac{\beta}{2}\left|\nabla\phi\right|^{2}\right)\,\mathrm{d}{\mathbf{x}}+u^{2}-\big{(}\frac{\beta}{2}+\frac{\beta^{2}}{4}\big{)}\left|\Omega\right|-C_{0}.$$

$${\partial }_t^{\alpha }\varphi$$

$$u_t$$

$$\displaystyle U\left(\phi\right)=\frac{\nabla\cdot\big{(}\big{(}\left|\nabla\phi\right|^{2}-1-\beta\big{)}\nabla\phi\big{)}}{\sqrt{\int_{\Omega}\frac{1}{4}\left(\left|\nabla\phi\right|^{2}-1-\beta\right)^{2}\,\mathrm{d}{\mathbf{x}}+C_{0}}}.$$

$$\varphi \left(\mathbf{x},0\right)={\varphi }_0\left(\mathbf{x}\right)\text{and}u\left(0\right)=u\left({\varphi }_0\left(\mathbf{x}\right)\right),$$

$$\displaystyle E\left[\phi(T),u(T)\right]-E\left[\phi(0),u(0)\right]=-\frac{1}{M}\int_{\Omega}\int_{0}^{T}\big{(}\partial_{t}^{\alpha}\phi\big{)}\partial_{t}\phi\,\mathrm{d}{t}\,\mathrm{d}{\mathbf{x}}\leqslant{0}.$$

$$\displaystyle v(t)=\sqrt{\int_{\Omega}\left(\frac{1}{2}\ln\big{(}1+\left|\nabla\phi\right|^{2}\big{)}+\frac{\beta}{2}\left|\nabla\phi\right|^{2}\right)\,\mathrm{d}{\mathbf{x}}+C_{0}},$$

$$E\left[\varphi ,v\right]={\int }_{\Omega }\left(\frac{{\epsilon }^2}{2}{\left|\Delta \varphi \right|}^2+\frac{\beta }{2}{\left|\nabla \varphi \right|}^2\right)\mathrm{d}\mathbf{x}-v^2+C_0.$$

$${\partial }_t^{\alpha }\varphi$$

$$v_t$$

$$\displaystyle V(\phi)=\frac{\nabla\cdot\left(\big{(}\frac{1}{1+\left|\nabla\phi\right|^{2}}+\beta\big{)}\nabla\phi\right)}{\sqrt{\int_{\Omega}\left(\frac{1}{2}\ln\big{(}1+\left|\nabla\phi\right|^{2}\big{)}+\frac{\beta}{2}\left|\nabla\phi\right|^{2}\right)\,\mathrm{d}{\mathbf{x}}+C_{0}}},$$

$$\varphi \left(\mathbf{x},0\right)={\varphi }_0\left(\mathbf{x}\right)\text{and}v(0)=v\left({\varphi }_0\left(\mathbf{x}\right)\right).$$

$$\displaystyle E\left[\phi(T),v(T)\right]-E\left[\phi(0),v(0)\right]=-\frac{1}{M}\int_{\Omega}\int_{0}^{T}\big{(}\partial_{t}^{\alpha}\phi\big{)}\partial_{t}\phi\,\mathrm{d}{t}\,\mathrm{d}{\mathbf{x}}\leqslant{0}.$$

$$({\partial }_{\tau }^{\alpha }v{)}^n:={\int }_{t_0}^{t_n}{\omega }_{1-\alpha }(t_n-s)({\Pi }_{1}v{)}^{\prime }(s)\mathrm{d}s=\sum _{k=1}^n{a}_{n-k}^{(n)}{▽}_{\tau }{v}^k,$$

$$a_{n-k}^{(n)}:=\frac{1}{{\tau }_k}{\int }_{t_{k-1}}^{t_k}{\omega }_{1-\alpha }(t_n-s)\mathrm{d}s\text{for 1⩽k⩽n.}$$

$$a_{n-k}^{(n)}>0\text{and}{a}_{n-k-1}^{(n)}>a_{n-k}^{(n)}\text{for 1⩽k⩽n−1.}$$

$${\left({\partial }_{\tau }^{\alpha }\varphi \right)}^n=-M\left({\epsilon }^2{\Delta }^2{\varphi }^n+f\left(\nabla {\varphi }^{n-1}\right)\right)+S(\Delta {\varphi }^n-\Delta {\varphi }^{n-1}).$$

$$a_{n-k}^{(n)}=a_{n-k}=\frac{1}{{\tau }^{\alpha }}\left[{\omega }_{2-\alpha }(n-k+1)-{\omega }_{2-\alpha }(n-k)\right]\text{for 1⩽k⩽n,}$$

Full text

Keywords:   Time-fractional molecular beam epitaxial; novel L1-type formula; scalar auxiliary variable approach; unconditional energy stability; adaptive time-stepping strategy

AMS subject classiffications.   35Q99, 65M06, 65M12, 74A50

1 Introduction

The Molecular Beam Epitaxial (MBE) growth models, in recent years, have become a powerful new technique in material science, such as making compound semi-conductor manufacture with great precision and high purity. Also, this technique is broadly applied to investigate thin-film deposition of single crystal. Roughly speaking, the mathematical models used in previous works to study dynamics of the MBE growth process can be classified into three broad categories known as: atomistic models that are performed using the form of molecular dynamics [1]; continuum models in the form of partial differential systems [2]; and hybrid models [3] which can be regarded as a compromise in the light of the models mentioned above.

In current work, we are interested in the continuum model for the evolution of the MBE growth that is derived by energy variational strategy. More precisely, given the effective free energy of the model $E[\phi]$ and associated with $L^{2}$ gradient flow, the height evolution model could be written as

[TABLE]

in which $\phi$ is a scaled height function of a thin film, positive constant $M$ is the mobility coefficient and $\mu=\frac{\delta{E}}{\delta\phi}$ is the variational derivative of the free energy $E$. One intrinsic property of above system is the energy dissipation law, that is,

[TABLE]

where the $L^{2}$ inner product is defined by $\left(f,g\right)=\int_{\Omega}fg\,\mathrm{d}{\mathbf{x}}$, and the $L^{2}$ norm $\left\|f\right\|=\sqrt{\int_{\Omega}f^{2}\,\mathrm{d}{\mathbf{x}}}$, for all $f,g\in{L}^{2}(\Omega)$. Notice that the periodic boundary condition or any other proper boundary that can satisfy the flux free condition at the boundary $\partial_{\mathbf{n}}\phi|_{\partial\Omega}=0$ and $\partial_{\mathbf{n}}\Delta\phi|_{\partial\Omega}=0$ are chose to ensure that the boundary integrals resulted during the integration by parts vanish, where $\mathbf{n}$ is the outward normal on the boundary.

There are a great amount of works contributed to the investigation of numerical approximations for the solution of the phase field models, for instance the so-called convex splitting technique [4] and the stabilized semi-implicit method [5]. More recently, there are two novel strategies that are used to design second-order, unconditionally energy stable numerical schemes to solve the phase field models: the invariant energy quadratization (IEQ) method [6] and the scalar auxiliary variable (SAV) approach [7]. The latter is inspired by the former while leads to numerical schemes that only decoupled equations with constant coefficients need to be solved at each time step. Notwithstanding, the common goal of IEQ and SAV methods is to transform the primitive system into a new equivalent system with a quadratic energy functional and the corresponding modified energy dissipation property. Rigorous analysis and comparisons between IEQ and SAV might be out of this article’s scope, we refer to [6, 8, 7] for more details.

In comparison with the bright achievement of classical phase field models, there are many researches on building fractional phase field models to better model the anomalously diffusive effects. For instance, the time, space and time-space fractional Allen-Cahn equations were suggested [9, 10, 11, 12] to accurately describe anomalous diffusion problems. Li et al. [9] investigated a space-time fractional Allen-Cahn phase field model that describes the transport of the fluid mixture of two immiscible fluid phases. They concluded that the alternative model could provide more accurate description of anomalous diffusion processes and sharper interfaces than the classical model. Hou et al. [10] showed that the space-fractional Allen-Cahn equation could be viewed a $L^{2}$ gradient flow for the fractional analogue version of Ginzburg-Landau free energy function. Also, the authors proved the energy decay property and the maximum principle of continuous problem. Tang et al. [13] proved the time-fractional phase field models indeed admit an energy dissipation law of an integral type. Meanwhile, they applied the uniform L1 formula to construct a class of finite difference schemes, which can inherit the theoretical energy dissipation property. Along the numerical front, Zhao et al. [11, 12] studied a series of the time-fractional phase field models numerically, covering the time-fractional Cahn-Hilliard equation with different types of variable mobilities and time-fractional molecular beam epitaxy model. The considerable numerical evidences indicate that the effective free energy or roughness of the time-fractional phase field models during coarsening obeys a similar power scaling law as the integer ones, where the power is linearly proportional to the fractional index $\alpha$. In other words, the main difference between the time-fractional phase field models and integer ones lie in the timescales of coarsening.

The multi-scale nature of time-fractional phase field models prompts us to construct reliable time-stepping methods on general nonuniform meshes. In this paper, nonuniform time-stepping schemes are investigated for the time-fractional MBE model

[TABLE]

where the notation $\partial_{t}^{\alpha}:={}_{0}^{C}\!D_{t}^{\alpha}$ in (1.3) denotes the fractional Caputo derivative of order $\alpha$ with respect to $t$,

[TABLE]

involving the fractional Riemann-Liouville integral $\mathcal{I}_{t}^{\beta}$ of order $\beta>0$, that is,

[TABLE]

Specifically, in comparison to the decay property (1.2) of the classical model, the energy dissipation law of the time-fractional MBE model (1.3), see also [13], is given by

[TABLE]

To our knowledge, there are few results in the literature on the discrete energy decay laws of numerical approaches for the time-fractional phase field models, especially on nonuniform time meshes. One of our interests in this paper is to build nonuniform time-stepping methods preserving the energy dissipation law of the problem (1.3) in discrete sense.

We consider the nonuniform time levels $0=t_{0}<t_{1}<\cdots<t_{k-1}<t_{k}<\cdots<t_{N}=T$ with the time-step sizes $\tau_{k}:=t_{k}-t_{k-1}$ for $1\leqslant{k}\leqslant{N}$ and the maximum time-step size $\tau:=\max_{1\leqslant{k}\leqslant{N}}\tau_{k}$. Also, let the local time-step ratio $\rho_{k}:=\tau_{k}/\tau_{k+1}$ and the maximum step ratio $\rho:=\max_{k\geq 1}\rho_{k}$. Given a grid function $\{v^{k}\}$, put $\triangledown_{\tau}v^{k}:=v^{k}-v^{k-1}$, $\partial_{\tau}v^{k-\frac{1}{2}}:=\triangledown_{\tau}v^{k}/\tau_{k}$ and $v^{k-\frac{1}{2}}:=(v^{k}+v^{k-1})/2$ for $k\geq{1}$. Always, let $(\Pi_{1,k}v)(t)$ denote the linear interpolant of a function $v(t)$ at two nodes $t_{k-1}$ and $t_{k}$, and define a piecewise linear approximation

[TABLE]

As an essential mathematical feature of linear and nonlinear subdiffusion problems including the time-fractional MBE model (1.3), the solution always lacks the smoothness near the initial time although it would be smooth away from $t=0$, see [14, 15]. Actually, in any numerical methods for solving time-fractional diffusion equations, a key consideration is the initial singularity of solution, see the recent works [16, 17, 18]. More directly, we apply a new L1-type formula to the time-fractional problem (1.3), see more details in subsection 3.1 and Example 4.3. Figure 1 plots the discrete time derivative $\partial_{\tau}u^{k-\frac{1}{2}}$ near $t=0$ on the graded mesh $t_{k}=(k/N)^{\gamma}$ when $\alpha=0.4$. They suggest that

[TABLE]

It says that the solution possesses weakly singularity like $u_{t}=O(t^{\alpha-1})$ near the initial time, which can be alleviated by using the graded meshes. Thus, the second interest of this paper is to resolve the essentially weak singularity in the equation (1.3) by refining time mesh. Actually, we will show that the graded mesh can recover the optimal time accuracy when the solution $u$ is non-smooth near $t=0$.

In the next section, we reformulate the time-fractional MBE model (1.3) by using SAV technique, which provides an elegant platform for constructing energy stable schemes. Although the L1 formula was applied in [13] to handle the time-fractional phase filed models and the discrete energy decay laws were also established on uniform time mesh, we are not able to exploit any discrete energy dissipation laws under the adaptive time-steps and leave it as an open problem (see Remark 1). In section 3, this open problem guides us to develop a novel numerical Caputo derivative (3.7), called L1*+* formula, which preserves inherently the positive semi-definite property of integral kernel, see (3.9). This new formula is as simple as the classical L1 formula because it is also constructed by the piecewise linear interpolation; however it can achieve second-order accurate for all fractional orders $\alpha\in(0,1)$ and substantially improves the accuracy of L1 formula, especially when the fractional order $\alpha\rightarrow 1$. Applying the new formula to the time-fractional MBE model (1.3), we obtain two Crank-Nicolson SAV schemes (3.12)-(3.13) and (3.17)-(3.18), respectively. Thus it is straightforward to confirm that both of them preserve the discrete energy dissipation law well, see Theorems 3.1-3.2, so that they are unconditionally stable in the energy norm. At last, we also apply the sum-of-exponentials (SOE) technique to develop a fast L1*+* algorithm in section 4. By using an accuracy criterion based adaptive time-stepping strategy, extensive numerical experiments are curried out to show the effectiveness of our numerical approaches and to support our analysis.

In summary, the main contributions of this paper for the time-fractional MBE model (1.3) are the following: suggest a simple L1*+* formula of Caputo derivative having second-order accuracy for any fractional order $\alpha\in(0,1)$; apply it to build two Crank-Nicolson SAV schemes preserving the discrete energy dissipation law; and develop a fast L1*+* algorithm incorporating adaptive time steps to speed up the long-time, multi-scale simulations.

2 Equivalent PDE systems

Let $\varepsilon>0$ is a positive constant and $F\left(\mathbf{v}\right)$ is a nonlinear, smooth function of its argument $\mathbf{v}\in\mathbb{R}^{d}\,(d=1,2,3)$. The Ehrlich-Schwoebel energy $E[\phi]$ is given by

[TABLE]

There are two popular choices of the nonlinear bulk potential: the double well potential, $F\left(\nabla\phi\right)=\frac{1}{4}(|\nabla\phi|^{2}-1)^{2}$, for the case of slope selection model; and the logarithmic potential, $F\left(\nabla\phi\right)=-\frac{1}{2}\ln(1+|\nabla\phi|^{2})$, for the case without slope selection model. Correspondingly, the governing equation of (1.3) for the height function $\phi$ reads:

[TABLE]

in which the case of slope selection model

[TABLE]

while the case without slope selection model

[TABLE]

In what follows, we denote the time fractional MBE model (2.2) with and without slope selection by “the Slope-Model” and “the No-Slope-Model”, respectively, for simplicity. Additionally, boundary conditions are set to be periodic so as not to complicate the analysis with unwanted details. We then reformulate the MBE models into the equivalent PDE systems by borrowing the ideas from the IEQ and SAV methods [6, 8, 7].

2.1 The Slope-Model and its equivalent system

We introduce a scalar auxiliary function $u(t)$ in term of original variable $\phi$ given by

[TABLE]

where $C_{0}>0$ is constant that ensures the radicand positive. Compared with those in [6, 7], the scalar auxiliary variable $u(t)$ adds an artificial parameter $\beta$ to regularize the numerical approach. As a consequence, the free energy of the Slope-Model is transformed into a quadratic form

[TABLE]

Correspondingly, the Slope-Model can be reformulated to the following equivalent form

[TABLE]

where the expression of notation $U\left(\phi\right)$ given by

[TABLE]

We should note in passing that the new system is subjected to the initial conditions

[TABLE]

and the boundary conditions are same as the primitive model. By taking the $L^{2}$ inner product of (2.7) with $\phi$, of (2.8) with $u$, and then making time integrations on both sides, we see the equivalent system preserves the energy dissipation law

[TABLE]

The non-positive of the right part of above equality is determined by [13, Lemma 2.1].

2.2 The No-Slope-Model and its equivalent system

For the No-Slope-Model, we introduce a scalar auxiliary function $v(t)$ in term of original variable $\phi$ as follows

[TABLE]

where $\beta$ and $C_{0}$ are similar to the previous ones. Therefore, the free energy of the No-Slope-Model could be rewritten into

[TABLE]

We then could rewrite the No-Slope-Model as an equivalent form

[TABLE]

in which

[TABLE]

with the following initial conditions

[TABLE]

Similarly, the new system admits the following energy dissipation law

[TABLE]

3 Novel L1*+* formula and Crank-Nicolson SAV schemes

The well-known L1 formula of Caputo derivative is given by

[TABLE]

where the corresponding discrete convolution kernels $a_{n-k}^{(n)}$ are defined by

[TABLE]

Obviously, the discrete L1 kernels $a_{n-k}^{(n)}$ are positive and decreasing, see also [16, 17],

[TABLE]

Based on the above L1 formula on uniform time mesh, a linearized scheme by using the stabilized technique via a stabilized term $S(\Delta\phi^{n}-\Delta\phi^{n-1})$ for a properly large scalar parameter $S>0$ is presented. We refer to [13] for more details. The resulting stabilized semi-implicit scheme for the problem (2.2) reads

[TABLE]

It is to mention that, when the time mesh is uniform such that the discrete L1 kernels

[TABLE]

the semi-implicit scheme (3.4) inherits a discrete energy dissipation law, see more details in [13, Theorem 3.3]. As seen, the proof of discrete version of energy dissipation law (1.6) relies on the positive semi-definite property of a discrete quadratic form $\sum_{k=1}^{n}w_{k}\sum_{j=1}^{k}a_{k-j}w_{j}\geq 0$.

Remark 1

It seems rather difficult to extend the positive semi-definite property to a general class of nonuniform meshes. More precisely, we are not able to verify the positive semi-definite property of the following quadratic form (by taking $w_{k}=\triangledown_{\tau}v^{k}$)

[TABLE]

The open problem in Remark 1 motivates us to design a novel discrete Caputo formula such that it naturally possesses the energy decay law (1.6) in the discrete sense for the time-fractional MBE model. Alternatively, the novel discrete Caputo formula should inherit the positive semi-definite property of a quadratic form like (3.5).

Fortunately, as pointed out early in [19, 20], we know that the weakly singular kernel $\omega_{1-\alpha}$ is positive semi-definite, that is,

[TABLE]

for $t>0$ and any $w\in{C}[0,T]$. Actually, one has $(-1)^{k}\omega_{1-\alpha}^{(k)}(t)\geq{0}$ for any integer $k\geq 1$ and $\int_{0}^{\infty}\omega_{1-\alpha}(t)\cos(\theta{t})\,\mathrm{d}{t}\geq{0}$ for any $\theta>0$, then the Plancherel’s theorem implies the positive semi-definite property (3). We will see that a discrete counterpart of (3) yields a discrete Caputo approximation preserving the desired energy dissipation property (1.6) when it is applied to time-fractional MBE model (2.2).

3.1 The L1*+* formula

The L1*+* formula for the Caputo derivative (1.4) is defined at time $t=t_{n-\frac{1}{2}}$ as follows

[TABLE]

where the discrete convolution kernels $\bar{a}_{n-k}^{(n)}$ are defined by

[TABLE]

Obviously, the naturally nonuniform L1*+* approximation $(\partial_{\tau}^{\alpha}v)^{n-\frac{1}{2}}\approx\frac{1}{\tau_{n}}\int_{t_{n-1}}^{t_{n}}\left(\partial_{t}^{\alpha}v\right)(t)\,\mathrm{d}{t}$ ensures the positive semi-definite property (3) by taking $w=\Pi_{1}v$, that is,

[TABLE]

The definition (3.8) and the arbitrariness of function $v$ yield the following result.

Lemma 3.1

The discrete convolution kernels $\bar{a}_{n-k}^{(n)}$ in (3.8) are positive, and for any real sequence $\{w_{k}\}_{k=1}^{n}$ with $n$ entries, it holds that

[TABLE]

Before applying the L1*+* formula (3.7) to the time-fractional MBE model (2.2), we show that it is a second-order approximation for the Caputo derivative (1.4) numerically. Consider a simple fractional ODE problem $\partial_{t}^{\alpha}u(t)=f(t)$ for $0<t<1$, we run a Crank-Nicolson-type scheme $(\partial_{\tau}^{\alpha}u)^{n-\frac{1}{2}}=f(t_{n-\frac{1}{2}})$ with uniform time-steps $\tau_{k}=\tau$, by choosing a smooth solution $u=\omega_{1+\sigma}(t)$ with the regularity parameter $\sigma=2.5$. The discrete $l^{\infty}$ norm errors $e(N)=\max_{1\leqslant{n}\leqslant{N}}|u(t_{n})-u^{n}|$ are listed in Table 1. It seems that the numerical accuracy of (3.7) is second-order accurate for any fractional order $\alpha\in(0,1)$. If the solution has an initial singularity, the L1*+* formula (3.7) can also achieve the second-order accuracy by properly refining the mesh near $t=0$, see more tests in Example 4.1.

3.2 Crank-Nicolson SAV schemes preserving energy dissipation

In what follows, we concern only with the time discretization of the equivalent systems, while the spatial approximation can be diverse, examples as finite difference, finite element or spectral methods. Integrating the equations (2.7)-(2.8) from $t=t_{n-1}$ to $t_{n}$, respectively, leads to

[TABLE]

Applying the L1*+* formula (3.7), the trapezoidal formula, we have the following Crank-Nicolson SAV (CN-SAV) time-stepping scheme for the Slope-Model

[TABLE]

where $\hat{\phi}^{n-\frac{1}{2}}:=\phi^{n-1}+\triangledown_{\tau}\phi^{n-1}/(2\rho_{n-1})$ is the local extrapolation.

Note that, the construction of L1*+* formula (3.7) implies that the above CN-SAV scheme (3.12)-(3.13) is naturally suitable for a general class of nonuniform meshes. Moreover, next result shows that it is unconditionally energy stable.

Theorem 3.1

The CN-SAV scheme (3.12)-(3.13) preserves the energy dissipation law,

[TABLE]

such that it is unconditionally stable, where

[TABLE]

**Proof ** Taking the inner product of (3.12)-(3.13) with $\triangledown_{\tau}\phi^{n}$ and $2\tau_{n}u^{n-\frac{1}{2}}$, respectively, and adding the resulting two identities, we obtain

[TABLE]

which implies that

[TABLE]

By summing the superscript $k$ from $1$ to $n$, we apply the property (3.9) to derive that

[TABLE]

It completes the proof.

Analogously, for the No-Slope-Model (2.13)-(2.14), we present the following CN-SAV time-stepping scheme

[TABLE]

For the above numerical scheme, we shall only state the energy dissipation law below, as their proofs are essentially the same as Theorem 3.1.

Theorem 3.2

The CN-SAV scheme (3.17)-(3.18) satisfies the following discrete energy dissipation law,

[TABLE]

where

[TABLE]

3.3 Further notes on L1*+* formula

Remark 2** (Multi-term and distributed time-fractional problems)**

As seen in Table 1 and more tests for Example presented in next section, the accuracy of the L1+ formula (3.7) is dependent on the regularity of solution, but would be independent of the fractional order $\alpha\in(0,1)$. It is quite different from some exiting numerical approaches, such as L1, Alikhanov [21, 22] and BDF2-like formulas [23, 24, 25], of which the consistency errors are dependent on the fractional order $\alpha$. This feature is attractive for further applications in developing second-order approximations for multi-term and distributed-order fractional diffusion equations. For an example, consider a simple multi-term fractional diffusion problem,

[TABLE]

where $0<\alpha_{i}<1$ for $1\leqslant i\leqslant m$, and $w_{i}$ is the corresponding weights. One can construct the following second-order Crank-Nicolson-type time-stepping scheme

[TABLE]

Obviously, the L1+ formula (3.7) would be also useful in approximating the distributed-order Caputo derivative since it can be approximated by certain multi-term derivative via some proper quadrature rule, see [25].

The suggested L1*+* formula (3.7) seems very promising in further applications for other time-fractional field phase models and other time-fractional differential equations. They make the rigorous theoretical analysis on consistency, stability and convergence very important, especially on a general class of nonuniform meshes.

However, the established theory [16, 26, 17, 18] for nonuniform L1 and L1-2σ (Alikhanov) formulas can not be applied to the L1*+* formula directly because the corresponding discrete convolution kernels $\bar{a}_{n-k}^{(n)}$ in (3.8) do not have the uniform monotonicity like (3.2). Actually, the definition (3.8) and the integral mean-value theorem yield the following result.

Lemma 3.2

The positive discrete kernels $\bar{a}_{n-k}^{(n)}$ in (3.8) fulfill

[TABLE]

Notwithstanding, it ought to be emphasized that

[TABLE]

It is easily seen that $\bar{a}_{0}^{(n)}<\bar{a}_{1}^{(n)}$ as $\alpha\rightarrow{0}$ and $\bar{a}_{0}^{(n)}>\bar{a}_{1}^{(n)}$ as $\alpha\rightarrow{1}$. So the value of $\bar{a}_{0}^{(n)}-\bar{a}_{1}^{(n)}$ may change the sign when the fractional order $\alpha$ varies over $(0,1)$. At the same time, this situation is no worse than the case in the BDF2-like formulas [24, 25] in which the second kernel would be negative when $\alpha\rightarrow{1}$, see more details and a potential remedy technique in [26, Remark 6]. The theoretical investigations, including the consistency, stability and convergence, of nonuniform L1*+* formula (3.7) will be addressed in a separate technical report.

4 Numerical algorithms and examples

4.1 A fast version of L1*+* formula

It is evident that the approximations (3.1) or (3.7) are prohibitively expensive for long time simulations due to the long-time memory. Therefore, to reduce the computational cost and storage requirements, we apply the sum-of-exponentials (SOE) technique to speed up the evaluation of the L1*+* formula. A core result is to approximate the kernel function $\omega_{1-\alpha}(t)$ efficiently on the interval $[\Delta{t},\,T]$, see [27, Theorem 2.5].

Lemma 4.1

For the given $\alpha\in(0,\,1)$, an absolute tolerance error $\epsilon\ll{1}$, a cut-off time $\Delta{t}>0$ and a finial time $T$, there exists a positive integer $N_{q}$, positive quadrature nodes $\theta^{\ell}$ and corresponding positive weights $\varpi^{\ell}\,(1\leqslant{\ell}\leqslant{N_{q}})$ such that

[TABLE]

The Caputo derivative (1.4) is split into the sum of a history part (an integral over $[0,\,t_{n-1}]$) and a local part (an integral over $[t_{n-1},\,t_{n}]$) at the time $t_{n}$, see also [17]. Then, the local part will be approximated by linear interpolation directly, the history part can be evaluated via the SOE technique, that is,

[TABLE]

where

[TABLE]

By utilizing the linear interpolation and a recursive formula, we can approximate $\mathcal{H}^{\ell}(t_{k})$ by

[TABLE]

where the positive coefficients

[TABLE]

Having taken this excursion through (4.1)-(4.1), we arrive at the fast L1*+* formula

[TABLE]

in which $\mathcal{H}^{\ell}(t_{k})$ is computed by using the recursive relationship

[TABLE]

4.2 Adaptive time-stepping strategy

In the previous sections, we have proved that the numerical schemes are unconditionally energy stable which implies large time steps are allowed. Indeed, in simulating the phase field problems such as the coarsening dynamics problems discussed in Example 4.4, adaptive time-stepping strategy is necessary to efficiently resolve widely varying time scales and to significantly reduce the computational cost. Roughly speaking, the adaptive time steps can be selected by using an accuracy criterion example as [28], or the time evolution of the total energy such as [29]. We focus on the former and update the time step size by using the formula

[TABLE]

where $\rho$ is a default safety coefficient, $tol$ is a reference tolerance, and $e$ is the relative error at each time level. The details of the adaptive time steps strategy are presented in Algorithm 1. Here, the first-order SAV and second-order SAV schemes refer to the backward Euler-L1 method and Crank-Nicolson SAV method proposed in this article, respectively.

4.3 Numerical examples

The CN-SAV methods (3.12)-(3.13) and (3.17)-(3.18) are examined in this section for the time-fractional MBE model (2.2). Specially, the model fractional ODE problem $\partial_{t}^{\alpha}u(t)=f(t)$ is used in Example 4.1 to test the accuracy of L1*+* formula (3.7). The time interval $[0,T]$ is always divided into two parts $[0,T_{0}]$ and $[T_{0},T]$ with total $N$ subintervals. We will take $T_{0}=\min\{1/\gamma,T\}$, and apply the graded grid $t_{k}=T_{0}(k/N_{0})^{\gamma}$ in $[0,T_{0}]$ to resolve the initial singularity. In the remainder interval $[T_{0},T]$, we put $N_{1}:=N-N_{0}$ cells with random time-steps

[TABLE]

where $\epsilon_{k}\in(0,1)$ are the random numbers.

The time accuracy of the proposed methods is mainly focused on, so for simplicity, and the Fourier pseudo-spectral method is always applied to approximate the space variables using the same spacing in each spatial direction. To examine the CN-SAV schemes, the maximum norm error $e(N):=\max_{1\leqslant{n}\leqslant{N}}\|U^{n}-u^{n}\|_{\infty}$ is recorded in each run, and the experimental convergence order in time is computed by

[TABLE]

where $\tau(N)$ denotes the maximum time-step size for total $N$ subintervals

Example 4.1

Consider the model fractional ODE problem $\partial_{t}^{\alpha}u(t)=f(t)$ with an exact solution $u=\omega_{1+\sigma}(t)$, where the parameter $\sigma$ determines the initial regularity of $u$. The accuracy of L1+ formula (3.7) is examined carefully using the following three scenarios:

1. (a)

$\sigma>{2}$* and $\alpha=0.1,\,0.5$ and $0.9$ on uniform mesh, see Table 1 in subsection 3.1;* 2. (b)

$\sigma<{2}$* and $\alpha=0.1,\,0.5$ and $0.9$ using the graded parameter $\gamma=1$, see Table 2;* 3. (c)

$\sigma<{2}$* and $\alpha=0.3$ and $0.7$ using different graded parameters $\gamma\geq 1$, see Tables 3-4.*

Table 1 lists the case (a) having smooth solution, while Tables 2-4 record the two cases (b)-(c) having non-smooth solutions. From these numerical results in Tables 1-4, one sees that it is accurate of $O(\tau^{\max\{\gamma\sigma,2\}})$ via the following observations: (i) The numerical accuracy is independent of the fractional order $\alpha\in(0,1)$, and it is second-order accurate for smooth solutions with $\sigma\geq{2}$; (ii) On the uniform mesh, the numerical accuracy degenerates to $O(\tau^{\sigma})$ when the regularity parameter $\sigma\in(0,2)$; (iii) When the solution is non-smooth, the numerical accuracy reaches $O(\tau^{\gamma\sigma})$ by the graded mesh, and the second-order accuracy would be recovered by choosing $\gamma\geq{2/\sigma}$.

Example 4.2

To examine the temporal accuracy of our CN-SAV schemes, consider the time-fractional MBE model $\partial_{t}^{\alpha}\phi=-M\left(\varepsilon^{2}\Delta^{2}\phi+f\left(\nabla\phi\right)\right)+g(\mathbf{x},t)$ for $\mathbf{x}\in(0,2\pi)^{2}$ and $0<t<1$ such that it has an exact solution $\phi=\omega_{1+\sigma}(t)\sin(x)\sin(y)$.

The parameters are taken as $M=0.1,\beta=1,C_{0}=1,\varepsilon^{2}=0.5$, and the space is discretized by $128\times{128}$ meshes. We run the CN-SAV schemes (3.12)-(3.13) and (3.17)-(3.18) by setting a variety of regularity parameters. Numerical results are tabulated in Tables 5-6, respectively. Tables 5 and 6 show the numerical results in the worse case of $\sigma<\alpha$. It seen that it is accurate of order $O(\tau^{\gamma\sigma})$ on the graded mesh, and the second-order accuracy can be recovered by taking $\gamma\geq\gamma_{\mathrm{opt}}=2/\sigma$. The computational results suggest that it is convergent of $O(\tau^{\min\{\gamma\sigma,2\}})$ in time although no theoretical proof is available up to now.

Example 4.3

We examine here the performance of the adaptive time-stepping method for the CN-SAV schemes (3.12)-(3.13) and (3.17)-(3.18) with parameter value $\epsilon^{2}=0.1$. To this end, we carry out a standard benchmark problem using initial condition

[TABLE]

In what follows, if not explicitly specified, the default values of parameters are given as $M=1,\beta=4,\mathbf{x}\in(0,2\pi)^{2},\rho=0.9,tol=10^{-3}$, and the space is discretized by $128\times{128}$ meshes. To quantify the deviation of the height function, define the roughness function $W(t)$:

[TABLE]

where $\bar{\phi}(t)=\frac{1}{\left|\Omega\right|}\int_{\Omega}\phi\left(\mathbf{x},t\right)\,\mathrm{d}{\mathbf{x}}$ is the average. The function $W(t)$ characterizes the mean size of the network cell.

For a fixed fractional index $\alpha=0.7$ and final time $T=30$, we first apply a constant time step $\tau=10^{-3}$, i.e., $N=30000$, to compute the solution. Recall that the intrinsically initial singularity of solution that presented in Figure 1 early, the numerical results suggest the time mesh should be refined near the initial time. As a consequence, we could obtain the reference solution where the parameter values $\gamma=3,T_{0}=0.01,N_{0}=30$ are applied in the cell $[0,T_{0}]$ and the uniform mesh is used over the remainder with the total numbers $N=30000$ as before. For the adaptive time-stepping technique, taking the analogous numerical strategy in the initial time and choosing parameter $\tau_{\min}=\tau_{N_{0}}=10^{-3}$ and $\tau_{\max}=10^{-1}$ for the remainder.

Figure 2 shows the plots of the energies and roughnesses for the Slope-Model and No-Slope-Model using uniform, graded and adaptive time steps, respectively. In our computation, we see that using uniform time step may produce incorrect steady-state solution while the energies still decay monotonically and the roughnesses show analogous plots compared with the reference solution. However, a comparison between the numerical results of using graded and adaptive time steps, we see that the curves are indistinguishable. For the adaptive strategy, the density of circle indicates the size of the time step, we then also observe that small time steps are used at the early stage of the computation because the quick transition of solution. Subsequent, large time steps are allowed due to the solution changes slowly. Again, small time steps are employed when capturing the steep structural transition from one stage to the next one. As a result, the total numbers of adaptive time steps are 4443 and 3586 for the Slope-Model and No-Slope-Model, respectively, while it takes 30000 constant time steps. The above observations show that the effect of the adaptive time approach on efficiency is significant dramatic.

Below we make a comparison among different fractional orders $\alpha$ to exploit how fractional index affects the evolution dynamics. Always, the third time mesh strategy is employed to solve the problem (2.2) with initial condition (4.5) in what follows. Figures 3 and 4 are the plots of the energies and roughnesses of the Slope-Model and No-Slope-Model for three different values of $\alpha$. We observe that in all cases the original energy decays rapidly and the smaller fractional order $\alpha$ the faster the energy dissipates, later it decays slower as smaller fractional index $\alpha$, and they reach the analogous steady-state in the end. The above observations may indicate that the time-fractional operator could affects the time scaling of the evolution dynamics, while the steady-state may not be affected.

Example 4.4

We investigate here the coarsening dynamics using the CN-SAV schemes (3.12)-(3.13) and (3.17)-(3.18). The initial condition is a random state by assigning a random number varying from $-0.001$ to $0.001$ to each grid points with parameter value $\epsilon=0.03$.

To discover the scaling of effective free energy and roughness in the time-fractional MBE model during coarsening, define the absolute value of the Slope-Model for each linear line of energy and roughness, $\beta(\alpha),R(\alpha)$, as

[TABLE]

and $E(\alpha,t)=\beta^{0}(\alpha)-\beta(\alpha)\log_{10}(t)$ for the No-Slope-Model, where $E(\alpha,t)$ and $W(\alpha,t)$ correspond to the energy $E(t)$ and roughness $W(t)$ with the fractional index $\alpha$, respectively.

Figures 5 and 6 show the time evolutions of energy and roughness with different fractional orders $\alpha$ for the Slope-Model during $t\in[1,500]$ where the adaptive time-stepping parameters $\tau_{\min}=\tau_{N_{0}}=1.25\times 10^{-4}$ and $\tau_{\max}=10^{-1}$. We observe that the energy dissipation approximately as $O(t^{-\frac{\alpha}{3}})$ and the growth rate of the roughness is approximately as $O(t^{\frac{\alpha}{3}})$ which are consistent with consistent with $O(t^{-\frac{1}{3}})$ and $O(t^{\frac{1}{3}})$, respectively, as recorded in [6, 7] for the integer Slope-Model. The time evolutions of energy and roughness for the No-Slope-Model using adaptive parameters $\tau_{\min}=\tau_{N_{0}}=3.32\times 10^{-5}$ and $\tau_{\max}=10^{-1}$ that are depicted in Figures 7 and 8, respectively. We see that the energy dissipation approximately as $O(-\log_{10}(t))$ and the growth rate of the roughness is approximately as $O(t^{\frac{\alpha}{2}})$ that are consistent with consistent with $O(-\log_{10}(t))$ and $O(t^{\frac{1}{2}})$ as reported in [6, 7] for the integer No-Slope-Model, respectively.

Furthermore, the drawings of the Figure LABEL:Slope-Dynamics-alpha-04-07-10 displays the numerical solutions of the height function $\phi$ and its Laplacian $\Delta\phi$ for the Slope-Model with different fractional orders $\alpha$. Based on Figure LABEL:Slope-Dynamics-alpha-04-07-10 and additional results not shown here for brevity, we observe that the edges of the pyramids generate a random distributed network over the domain and the pyramids become large when time increases. Additionally, coarsening dynamics, at beginning, appear to be faster as smaller $\alpha$ while it would be much slower as time evolves. Also, the observed phenomena are in good agreement with the published results [12].

5 Concluding remarks

In simulating the time-fractional phase field equations including the Molecular Beam Epitaxial model considered in this paper, the initial singularity should be treated properly because it always destroys the time accuracy of numerical algorithms especially near the initial time.

However, it seems challenging to build time-stepping approaches maintaining the discrete energy dissipation law based on the conventional L1 formula, especially on general nonuniform time meshes. Nonetheless, the energy stable schemes permitting adaptive time-stepping strategies are very attractive because they could be applicable for time-fractional phase-field models and for long-time simulations approaching the steady state. As an interesting remedy, the novel L1*+* formula, is proposed to approximate the fractional Caputo’s derivative. As a consequence, coupled with the scalar auxiliary variable, we suggest two linearized second-order energy stable CN-SAV schemes for the MBE model with and without slope selection, respectively, by virtue of the naturally positive semi-definite property of a discrete quadratic form. Furthermore, for the long-time simulations approaching the steady state, the fast L1*+* version incorporated with adaptive time-stepping strategy is developed for time-fractional phase field equations. Ample numerical examples are presented to validate the effectiveness of CN-SAV schemes.

The L1*+* formula would be superior to some widespread approximations, such as the L1, Alikhanov and BDF2-like formulas, because it is second-order accurate for both smooth and non-smooth solutions, and the convergence order is independent of the fractional order $\alpha\in(0,1)$. Thus it becomes critical to establish the rigorous theory on consistency, stability and convergence of nonuniform L1*+* formula. These issues will be addressed in a forthcoming report.