Adaptive second-order Crank-Nicolson time-stepping schemes for time fractional molecular beam epitaxial growth models
Bingquan Ji Hong-lin Liao Yuezheng Gong Luming Zhang





Adaptive second-order Crank-Nicolson time-stepping schemes for
 time fractional molecular beam epitaxial growth models
    Bingquan JiDepartment of Mathematics, Nanjing University of Aeronautics and Astronautics,
211101, P. R. China. Bingquan Ji (jibingquanm@163.com).
  Hong-lin LiaoORCID 0000-0003-0777-6832; Corresponding author. Department of Mathematics,
Nanjing University of Aeronautics and Astronautics,
Nanjing 211106, P. R. China. Hong-lin Liao (liaohl@csrc.ac.cn)
is supported by a grant 1008-56SYAH18037
from NUAA Scientific Research Starting Fund of Introduced Talent.
  Yuezheng Gong Department of Mathematics, Nanjing University of Aeronautics and Astronautics,
Nanjing 210016, P. R. China;
Yuezheng Gong (gongyuezheng@nuaa.edu.cn) is partially supported by the NSFC grant No. 11801269,
 and the NSF grant No. BK20180413 of Jiangsu Province.
  Luming ZhangDepartment of Mathematics, Nanjing University of Aeronautics and Astronautics,
211101, P. R. China. Luming Zhang (zhanglm@nuaa.edu.cn)
is supported by the NSFC grant No. 11571181.
    
===================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




  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.

  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



§ 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
<cit.>;

continuum models in the form of partial differential systems
<cit.>;

and hybrid models
<cit.>

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[ϕ]
and associated with L^2 gradient flow, the height evolution model could be written as

    ∂_tϕ
    =-Mμ,

in which ϕ is a scaled height function of a thin film, positive constant M is the mobility coefficient
and μ=δE/δϕ is the variational derivative of the free energy E.












One intrinsic property of above system is the energy dissipation law, that is,

    E/t
    =⟨δE/δϕ,
    ∂ϕ/∂t|
    =-1/M∂_tϕ^2≤0,

where the L^2 inner product is defined by
⟨f,g|=∫_Ωfg𝐱, and the L^2 norm f=√(∫_Ωf^2𝐱), for all f,g∈L^2(Ω).
Notice that the periodic boundary condition or any other proper boundary that can satisfy the flux free condition at the boundary
∂_𝐧ϕ|_∂Ω=0 and
∂_𝐧Δϕ|_∂Ω=0 are chose to ensure that the boundary integrals resulted during the integration by parts vanish, where 𝐧 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
<cit.>

and
the stabilized semi-implicit method
<cit.>.

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 <cit.> and the scalar auxiliary variable (SAV) approach <cit.>.
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
<cit.>

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
<cit.>
to accurately describe anomalous diffusion problems.
Li et al. <cit.> 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. <cit.> 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. <cit.> 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. <cit.> 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 α.
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









    ∂_t^αϕ
    =-Mμ,


where the notation ∂_t^α:=_0^CD_t^α in (<ref>) denotes
the fractional Caputo derivative of order α with respect to t,

    (∂_t^αv)(t)
    :=(ℐ_t^1-αv')(t)=∫_0^tω_1-α(t-s)v'(s)s,   0<α<1,

involving the fractional  Riemann-Liouville integral ℐ_t^β of order β>0, that is,

    (ℐ_t^βv)(t)
    :=∫_0^tω_β(t-s)v(s)s,  where  ω_β(t):=t^β-1/Γ(β).

Specifically, in comparison to the decay property (<ref>) of the classical model, the energy dissipation law of the time-fractional MBE model (<ref>), see also <cit.>, is given by

    Eϕ(T)
    -Eϕ(0)
    =-1/M∫_Ω∫_0^T∂_t^αϕ∂_tϕt𝐱≤0.

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 (<ref>) in discrete sense.

We consider the nonuniform time levels 0=t_0<t_1<⋯<t_k-1<t_k<⋯<t_N=T with the time-step sizes τ_k:=t_k-t_k-1 for 1≤k≤N and the maximum time-step size τ:=max_1≤k≤Nτ_k. Also, let the local time-step ratio ρ_k:=τ_k/τ_k+1 and the maximum step ratio ρ:=max_k≥ 1ρ_k. Given a grid function {v^k}, put ▿_τv^k:=v^k-v^k-1, ∂_τv^k-1/2:=▿_τv^k/τ_k and
v^k-1/2:=(v^k+v^k-1)/2 for k≥1. Always,
let (Π_1,kv)(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

    Π_1v:=Π_1,kv  so that  
    (Π_1v)'(t)=∂_τv^k-1/2  for t_k-1<t≤ t_k and k≥1.



















As an essential mathematical feature of linear and nonlinear subdiffusion problems
including the time-fractional MBE model (<ref>),
the solution always lacks the smoothness near the initial time
although it would be smooth away from  t=0, see <cit.>.










Actually, in any numerical methods for solving
time-fractional diffusion equations, a key consideration is the initial singularity of solution, see the recent works <cit.>.
More directly, we apply a new L1-type  formula to the time-fractional problem (<ref>), see more details in subsection 3.1 and Example <ref>.
Figure <ref>  plots the discrete time derivative ∂_τu^k-1/2
near t=0 on the graded mesh t_k=(k/N)^γ when α=0.4.
They suggest that

    log|u_t(𝐱,t)|
    ≈(α-1)log(t)+C(𝐱)  as t→0.

It says that the solution possesses weakly singularity like u_t=O(t^α-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 (<ref>) 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 (<ref>) by using SAV technique, which provides an elegant platform for constructing energy stable schemes.
Although the L1 formula was applied in <cit.> 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 <ref>).
In section 3, this open problem guides us to develop a novel numerical Caputo derivative (<ref>), called L1^+ formula,
which preserves inherently the positive semi-definite property of integral kernel, see (<ref>).
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 α∈(0,1)
and substantially improves the accuracy of L1 formula, especially when the fractional order α→1.
Applying the new formula to the time-fractional MBE model (<ref>),
we obtain two Crank-Nicolson SAV schemes (<ref>)-(<ref>) and
(<ref>)-(<ref>), respectively.
Thus it is straightforward to confirm that both of them preserve the discrete energy dissipation law well,
see Theorems <ref>-<ref>,
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 (<ref>) are the following: suggest a simple L1^+ formula of Caputo derivative having second-order accuracy for any fractional order α∈(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.




§ EQUIVALENT PDE SYSTEMS

Let ε>0 is a positive constant and F⟨𝐯| is a nonlinear, smooth function of its argument 𝐯∈ℝ^d (d=1,2,3).
The Ehrlich-Schwoebel energy E[ϕ] is given by


    E[ϕ]=∫_Ω⟨ε^2/2Δϕ^2
    + F⟨∇ϕ||𝐱 .

There are two popular choices of the nonlinear bulk potential:
the double well potential, F⟨∇ϕ|=1/4∇ϕ^2-1^2,
for the case of slope selection model;
and the logarithmic potential, F⟨∇ϕ|=-1/2ln1+∇ϕ^2,
for the case without slope selection model.
Correspondingly, the governing equation of (<ref>) for the height function ϕ reads:

    ∂_t^αϕ
    =-MδE⟨ϕ|/δϕ
    =-M
    ⟨ε^2Δ^2ϕ
    +f⟨∇ϕ||,

in which the case of slope selection model

    f⟨∇ϕ|
    =-∇·∇ϕ^2-1∇ϕ,

while the case without slope selection model

    f⟨∇ϕ|
    =∇·1/1+∇ϕ^2∇ϕ.

In what follows, we denote the time fractional MBE model (<ref>) 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 <cit.>.



 §.§ The Slope-Model and its equivalent system

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

    u(t)
    =√(∫_Ω1/4⟨∇ϕ^2-1-β|^2𝐱+C_0),

where C_0>0 is constant that ensures the radicand positive.
Compared with those in <cit.>, the scalar auxiliary variable u(t) adds
an artificial parameter β to regularize the numerical approach.
As a consequence, the free energy of the Slope-Model is transformed into a quadratic form

    Eϕ,u
    =∫_Ω⟨ε^2/2Δϕ^2
    +β/2∇ϕ^2|𝐱
    +u^2
    -β/2+β^2/4Ω-C_0.

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

    ∂_t^αϕ   =-M
    ⟨ε^2Δ^2ϕ
    -βΔϕ
    -U⟨ϕ|u|,
    
    u_t   =-1/2∫_ΩU⟨ϕ|∂_tϕ𝐱,

where the expression of notation U⟨ϕ| given by

    U⟨ϕ|
    =∇·∇ϕ^2-1-β∇ϕ/√(∫_Ω1/4⟨∇ϕ^2-1-β|^2𝐱+C_0).


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

    ϕ⟨𝐱,0|
    =ϕ_0⟨𝐱|  and  
    u⟨0|
    =u⟨ϕ_0⟨𝐱||,

and the boundary conditions are same as the primitive model.
By taking the L^2 inner product of (<ref>) with ϕ, of (<ref>) with u, and then making time integrations on both sides, we see the equivalent system preserves the energy dissipation law

    Eϕ(T),u(T)
    -Eϕ(0),u(0)
    =-1/M∫_Ω∫_0^T∂_t^αϕ∂_tϕt𝐱≤0.

The non-positive of the right part of above equality is determined by <cit.>.



 §.§ 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 ϕ as follows

    v(t)
    =√(∫_Ω⟨1/2ln1+∇ϕ^2
    +β/2∇ϕ^2|𝐱+C_0),

where β and C_0 are similar to the previous ones. Therefore, the free energy of the No-Slope-Model could be rewritten into

    Eϕ,v=∫_Ω⟨ε^2/2Δϕ^2
    +β/2∇ϕ^2|𝐱
    -v^2+C_0.

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

    ∂_t^αϕ   =-M
    ⟨ε^2Δ^2ϕ
    -βΔϕ
    +V(ϕ)v|,
    
    v_t   =1/2∫_Ω
    V⟨ϕ|∂_tϕ𝐱,

in which

    V(ϕ)=
    ∇·⟨1/1+∇ϕ^2+β∇ϕ|/√(∫_Ω⟨1/2ln1+∇ϕ^2
    +β/2∇ϕ^2|𝐱+C_0),

with the following initial conditions

    ϕ⟨𝐱,0|
    =ϕ_0⟨𝐱|  and  
    v(0)
    =v⟨ϕ_0⟨𝐱||.

Similarly, the new system admits the following energy dissipation law

    Eϕ(T),v(T)
    -Eϕ(0),v(0)
    =-1/M∫_Ω∫_0^T∂_t^αϕ∂_tϕt𝐱≤0.





§ NOVEL L1^+ FORMULA AND CRANK-NICOLSON SAV SCHEMES

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

    (∂_τ^αv)^n
    :=∫_t_0^t_nω_1-α(t_n-s)(Π_1v)'(s)s
    =∑_k=1^na_n-k^(n)▿_τv^k,

where the corresponding discrete convolution kernels a_n-k^(n) are defined by

    a_n-k^(n)
    :=1/τ_k∫_t_k-1^t_kω_1-α(t_n-s)s  for 1≤k≤n.


Obviously, the discrete L1 kernels a_n-k^(n) are positive and decreasing, see also <cit.>,

    a_n-k^(n)>0   and   a_n-k-1^(n)>a_n-k^(n)  for 1≤k≤n-1.



Based on the above L1 formula on uniform time mesh, a linearized scheme by using the stabilized technique
via a stabilized term S(Δϕ^n-Δϕ^n-1) for a properly large scalar parameter S>0 is presented. We refer to <cit.> for more details.
The resulting stabilized semi-implicit scheme for the problem (<ref>) reads

    ⟨∂_τ^αϕ|^n
    =-M⟨ε^2Δ^2ϕ^n
    +f⟨∇ϕ^n-1||
    +S(Δϕ^n-Δϕ^n-1).

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

    a_n-k^(n)=a_n-k
    =1/τ^αω_2-α(n-k+1)-ω_2-α(n-k)  for 1≤k≤n,

the semi-implicit scheme (<ref>)
inherits a discrete energy dissipation law, see more details in <cit.>.
As seen, the proof of discrete version of energy dissipation law
(<ref>) relies on
the positive semi-definite property of a discrete quadratic form ∑_k=1^nw_k∑_j=1^ka_k-jw_j≥0.


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=▿_τv^k)

    ∑_k=1^n▿_τv^k(∂_τ^αv)^k
    =∑_k=1^nw_k∑_j=1^ka_k-j^(k)w_j≥0.





The open problem in Remark <ref> motivates us to design a novel discrete Caputo formula such that
it naturally possesses the energy decay law (<ref>) 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 (<ref>).

Fortunately, as pointed out early in <cit.>,
we know that the
weakly singular kernel ω_1-α is positive semi-definite, that is,

    ℐ_t^1(w ℐ_t^1-αw)(t)
       =∫_0^tw(μ)μ∫_0^μω_1-α(μ-s)w(s)s
       =1/2∫_0^t∫_0^tω_1-α(|μ-s|)w(μ)w(s)μs≥0

for t>0 and any w∈C[0,T]. Actually, one has (-1)^kω_1-α^(k)(t)≥0 for any integer k≥1
and ∫_0^∞ω_1-α(t)cos(θt)t≥0 for any θ>0, then the Plancherel's theorem implies
the positive semi-definite property (<ref>).
We will see that a discrete counterpart of (<ref>) yields a discrete Caputo approximation preserving
the desired energy dissipation property (<ref>) when it is applied to
time-fractional MBE model (<ref>).




 §.§ The L1^+ formula

The L1^+ formula for the Caputo derivative (<ref>) is defined at time t=t_n-1/2 as follows

    (∂_τ^αv)^n-1/2
    :=1/τ_n∫_t_n-1^t_n∫_0^tω_1-α(t-s)(Π_1v)'(s)st
    =∑_k=1^na̅_n-k^(n)▿_τv^k  for n≥1,

where the discrete convolution kernels a̅_n-k^(n) are defined by

    a̅_n-k^(n)
    :=1/τ_nτ_k∫_t_n-1^t_n∫_t_k-1^min{t,t_k}ω_1-α(t-s)st  for 1≤k≤n.

Obviously, the naturally nonuniform L1^+ approximation (∂_τ^αv)^n-1/2≈1/τ_n∫_t_n-1^t_n⟨∂_t^αv|(t)t 

ensures the positive semi-definite property (<ref>)
by taking w=Π_1v, that is,

    ∑_k=1^n▿_τv^k(∂_τ^αv)^k-1/2   =∑_k=1^nτ_k(Π_1,kv)'(∂_τ^αv)^k-1/2
       =∫_t_0^t_n(Π_1v)'(t)
    ∫_0^tω_1-α(t-s)(Π_1v)'(s)st
       =ℐ_t^1(Π_1v)'
    ℐ_t^1-α(Π_1v)'(t_n)≥0  for n≥1.

The definition (<ref>) and the arbitrariness of function v yield the following result.

The discrete convolution kernels a̅_n-k^(n)
in (<ref>) are positive, and for any real sequence {w_k}_k=1^n with n entries, it holds that


    ∑_k=1^nw_k∑_j=1^ka̅_k-j^(k)w_j≥0  for n≥1.







Before applying the L1^+ formula (<ref>) to the time-fractional MBE model (<ref>),
we show that it is a second-order approximation for the Caputo derivative (<ref>) numerically.
Consider a simple fractional ODE problem
∂_t^αu(t)=f(t) for 0<t<1, we run a Crank-Nicolson-type scheme
(∂_τ^αu)^n-1/2=f(t_n-1/2)
with uniform time-steps τ_k=τ,
by choosing a smooth solution u=ω_1+σ(t) with the regularity parameter σ=2.5.
The discrete l^∞ norm errors e(N)=max_1≤n≤N|u(t_n)-u^n| are listed in Table <ref>.
It seems that the numerical accuracy of (<ref>) is second-order accurate for any fractional order α∈(0,1).
If the solution has an initial singularity, the L1^+ formula (<ref>) can also achieve the second-order accuracy
by properly refining the mesh near t=0, see more tests in Example <ref>.



 §.§ 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 (<ref>)-(<ref>) from t=t_n-1 to t_n,
respectively, leads to

    1/τ_n∫_t_n-1^t_n∂_t^αϕt   =-M/τ_n∫_t_n-1^t_n⟨ε^2Δ^2ϕ
    -βΔϕ
    -U⟨ϕ|u|t,
    1/τ_n∫_t_n-1^t_nu_tt   =-1/2τ_n∫_t_n-1^t_n∫_ΩU⟨ϕ|∂_tϕ𝐱t.

Applying the L1^+ formula (<ref>), the trapezoidal formula, we have the following Crank-Nicolson SAV (CN-SAV) time-stepping  scheme for the Slope-Model

    ⟨∂_τ^αϕ|^n-1/2   =-M
    ⟨ε^2Δ^2ϕ^n-1/2
    -βΔϕ^n-1/2
    -U(ϕ̂^n-1/2)u^n-1/2|,
    ∂_τu^n-1/2   =-1/2∫_Ω
    U(ϕ̂^n-1/2)
    ∂_τϕ^n-1/2𝐱,

where ϕ̂^n-1/2:=ϕ^n-1+▿_τϕ^n-1/(2ρ_n-1) is the local extrapolation.

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


The CN-SAV scheme (<ref>)-(<ref>) preserves the energy dissipation law,

    Eϕ^n,u^n
    -Eϕ^0,u^0≤0,  for   1≤n≤N,

such that it is unconditionally stable, where

    Eϕ^n,u^n
    =∫_Ω⟨ε^2/2Δϕ^n^2
    +β/2∇ϕ^n^2|𝐱
    +⟨u^n|^2-C_0.




Taking the inner product of (<ref>)-(<ref>) with ▿_τϕ^n and 2τ_nu^n-1/2,
respectively, and adding the resulting two identities, we obtain

    -1/M⟨⟨∂_τ^αϕ|^n-1/2,▿_τϕ^n|   =
    ⟨ε^2Δ^2ϕ^n-1/2
    -βΔϕ^n-1/2,▿_τϕ^n|
    +⟨u^n|^2-⟨u^n-1|^2,

which implies that

    Eϕ^k,u^k
    -Eϕ^k-1,u^k-1
    =-1/M⟨(∂_τ^αϕ)^k-1/2,▿_τϕ^k|  for 1≤k≤n.

By summing the superscript k from 1 to n, we apply the property (<ref>) to derive that

    Eϕ^n,u^n
    -Eϕ^0,u^0
    =-1/M∫_Ωℐ_t^1(Π_1ϕ)^'ℐ_t^1-α(Π_1ϕ)^'(t_n)𝐱≤0  for 1≤n≤N.

It completes the proof.



Analogously, for the No-Slope-Model (<ref>)-(<ref>), we present the following CN-SAV time-stepping  scheme

    ⟨∂_τ^αϕ|^n-1/2   =-M
    ⟨ε^2Δ^2ϕ^n-1/2
    -βΔϕ^n-1/2
    +V(ϕ̂^n-1/2)v^n-1/2|,
    ∂_τv^n-1/2   =1/2∫_ΩV(ϕ̂^n-1/2)
    ∂_τϕ^n-1/2𝐱.


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


The CN-SAV scheme (<ref>)-(<ref>) satisfies the following discrete energy dissipation law,

    Eϕ^n,v^n
    -Eϕ^0,v^0≤0,

where

    Eϕ^n,v^n
    =∫_Ω⟨ε^2/2Δϕ^n^2
    +β/2∇ϕ^n^2|𝐱
    -⟨v^n|^2+C_0.


































 §.§ Further notes on L1^+ formula













































As seen in Table  <ref> and more tests for Example
presented in next section,
the accuracy of the L1^+ formula (<ref>) is dependent on the regularity of solution,
but would be independent of the fractional order α∈(0,1).
It is quite different from some exiting numerical approaches, such as L1, Alikhanov <cit.>
and BDF2-like formulas <cit.>, of which the consistency errors are dependent on the fractional order α.
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,

    ∑_i=1^mw_i ∂_t^α_iu=Δ u+f,

where 0<α_i<1 for 1≤ i≤ m, and  w_i is the corresponding weights.
One can construct the following second-order Crank-Nicolson-type time-stepping scheme

    ∑_i=1^mw_i(∂_τ^α_iv)^n-1/2=Δ u^n-1/2+f^n-1/2  for 1≤ n≤ N.

Obviously, the L1^+ formula (<ref>) 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 <cit.>.


The suggested L1^+ formula (<ref>) 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 <cit.>
for nonuniform L1 and L1-2_σ (Alikhanov)  formulas can not be applied to the L1^+ formula directly
because the corresponding discrete convolution kernels a̅_n-k^(n) in (<ref>) do not have the uniform monotonicity like (<ref>).
Actually, the definition (<ref>) and the integral mean-value theorem
yield the following result.


The positive discrete kernels a̅_n-k^(n) in (<ref>) fulfill

    a̅_0^(n)=1/Γ(3-α)τ_n^α,  a̅_1^(n)>a̅_2^(n)>⋯>a̅_n-1^(n)>0  for n≥2.




Notwithstanding, it ought to be emphasized that

    a̅_0^(n)-a̅_1^(n)
    =1/Γ(3-α)τ_n^αρ_n-11+ρ_n-1+ρ_n-1^2-α-(1+ρ_n-1)^2-α.

It is easily seen that a̅_0^(n)<a̅_1^(n) as α→0 and a̅_0^(n)>a̅_1^(n) as α→1. So the value of a̅_0^(n)-a̅_1^(n) may change the sign when the fractional order α varies over (0,1). At the same time, this situation is no worse than the case in the BDF2-like formulas <cit.>
in which the second kernel would be negative when α→1, see more details
and a potential remedy technique in <cit.>.
The theoretical investigations,
including the consistency, stability and convergence, of nonuniform L1^+ formula (<ref>)
 will be addressed in a separate technical report.






§ NUMERICAL ALGORITHMS AND EXAMPLES



 §.§ A fast version of L1^+ formula

It is evident that the approximations (<ref>) or (<ref>) 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 ω_1-α(t) efficiently on the interval [Δt, T],
see <cit.>.

For the given α∈(0, 1), an absolute tolerance error ϵ≪1, a cut-off time Δt>0 and a finial time T, there exists a positive integer N_q, positive quadrature nodes θ^ℓ and corresponding positive weights ϖ^ℓ (1≤ℓ≤N_q) such that

    |
    ω_1-α(t)
    -∑_ℓ=1^N_qϖ^ℓe^-θ^ℓt|≤ϵ,
      ∀ t∈[Δt, T].




The Caputo derivative (<ref>) 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 <cit.>. Then, the local part will be approximated by linear interpolation directly, the history part can be evaluated via the SOE technique, that is,

    (∂_τ^αv)(t_n-1/2)
       ≈1/τ_n∫_t_n-1^t_n∫_0^t_n-1v^'(s)
    ∑_ℓ=1^N_qϖ^ℓe^-θ^ℓ(t-s)st
       +1/τ_n∫_t_n-1^t_n∫_t_n-1^tω_1-α(t-s)(Π_1v)^'(s)st
       =a̅_0^(n)▿_τv^n
    +1/τ_n∑_ℓ=1^N_qϖ^ℓ∫_t_n-1^t_n∫_0^t_n-1v^'(s)
    e^-θ^ℓ(t-t_n-1)
    e^-θ^ℓ(t_n-1-s)st
       =a̅_0^(n)▿_τv^n
    +1/τ_n∑_ℓ=1^N_qϖ^ℓ
    b^(n,l)ℋ^ℓ(t_n-1),   n≥1,

where

    ℋ^ℓ(t_k)
    :=∫_0^t_ke^-θ^ℓ(t_k-s)v^'(s)s  with  ℋ^ℓ(t_0)=0,  
    b^(k,l):=∫_t_k-1^t_ke^-θ^ℓ(t-t_k-1)t.

By utilizing the linear interpolation and a recursive formula, we can approximate ℋ^ℓ(t_k) by

    ℋ^ℓ(t_k)
       ≈∫_0^t_k-1e^-θ^ℓ(t_k-s)v^'(s)s
    +∫_t_k-1^t_ke^-θ^ℓ(t_k-s)▿_τv^k/τ_ks
       =e^-θ^ℓτ_kℋ^ℓ(t_k-1)
    +c^(k,l)▿_τv^k,

where the positive coefficients

    c^(k,l):=1/τ_k∫_t_k-1^t_ke^-θ^ℓ(t_k-s)s,   k≥1.

Having taken this excursion through (<ref>)-(<ref>),
 we arrive at the fast L1^+ formula

    (∂_f^αv)^n-1/2
    =a̅_0^(n)▿_τv^n
    +1/τ_n∑_ℓ=1^N_qϖ^ℓ
    b^(n,l)ℋ^ℓ(t_n-1),

in which ℋ^ℓ(t_k) is computed by using the recursive relationship

    ℋ^ℓ(t_k)
    =e^-θ^ℓτ_kℋ^ℓ(t_k-1)
    +c^(k,l)▿_τv^k,   k≥1,   1≤ℓ≤N_q.




 §.§ 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 <ref>,
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 <cit.>, or
the time evolution of the total energy such as <cit.>.
We focus on the former and update the time step size by using the formula

    τ_ada⟨e,τ|
    =ρ⟨tol/e|^1/2τ,

where ρ 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 <ref>.
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.




