# Highly efficient and energy dissipative schemes for the time fractional
Allen-Cahn equation∗

Dianming Hou1 Chuanju Xu2,3

###### Abstract.

In this paper, we propose and analyze a time-stepping method for the time
fractional Allen-Cahn equation. The key property of the proposed method is its
unconditional stability for general meshes, including the graded mesh commonly
used for this type of equations. The unconditional stability is proved through
establishing a discrete nonlocal free energy dispassion law, which is also
true for the continuous problem. The main idea used in the analysis is to
split the time fractional derivative into two parts: a local part and a
history part, which are discretized by the well known L1, L1-CN, and
$L1^{+}$-CN schemes. Then an extended auxiliary variable approach is used to
deal with the nonlinear and history term. The main contributions of the paper
are: First, it is found that the time fractional Allen-Chan equation is a
dissipative system related to a nonlocal free energy. Second, we construct
efficient time stepping schemes satisfying the same dissipation law at the
discrete level. In particular, we prove that the proposed schemes are
unconditionally stable for quite general meshes. Finally, the efficiency of
the proposed method is verified by a series of numerical experiments.

###### Key words and phrases:

time fractional Allen-Cahn, time-stepping scheme, unconditional stability,
spectral method

###### 1991 Mathematics Subject Classification:

65N35, 65M70, 45K05, 41A05, 41A10, 41A25

∗ The work of D. Hou is supported by NSFC grant 12001248 and the Natural
Science Foundation of the Jiangsu Higher Education Institutions of China grant
BK20201020. The second author has received support from NSFC grant 11971408,
NNW2018-ZT4A06 project, and NSFC/ANR joint program
51661135011/ANR-16-CE40-0026-01.
1School of Mathematics and Statistics, Jiangsu Normal University, 221116
Xuzhou, China.
2School of Mathematical Sciences and Fujian Provincial Key Laboratory of
Mathematical Modeling and High Performance Scientific Computing, Xiamen
University, 361005 Xiamen, China.
3Corresponding author. Email<EMAIL_ADDRESS>(C. Xu)

##  1\. Introduction

As a class of mathematical models, gradient flows is partial differential
equations under the form:

$\displaystyle\displaystyle\frac{\partial\phi}{\partial
t}=-\mbox{grad}_{H}E(\phi),$ (1.1)

where $\phi$ is the state function (also called phase function in many cases),
$E(\cdot)$ is the free energy driving functional associated to the physical
problem, and $\mbox{grad}_{H}E(\cdot)$ is the functional derivative of $E$ in
the Sobolev space $H$. It has other names: it is often called variational
principle in mathematics and Onsager principle in physics. Obviously the
gradient flows satisfies the energy dissipation law:

$\frac{d}{dt}E(\phi)=\Big{(}\mbox{grad}_{H}E(\phi),\frac{\partial\phi}{\partial
t}\Big{)}=-\|\frac{\partial\phi}{\partial t}\|_{0}^{2},$ (1.2)

where $(\cdot,\cdot)$ and $\|\cdot\|_{0}$ stand for the $L^{2}(\Omega)$-inner
product and norm, respectivly. This means that the state function $\phi$
evolves in such a way that the energy functional $E$ dissipates in time, i.e.,
in the opposite direction to the gradient of $E$ at $\phi$. This makes the
models very useful in many fields of science and engineering, such as
interface dynamics [4, 5, 18, 40], thin films [17, 30], crystal growth [13,
12, 14], polymers [15, 16] and liquid crystals [9, 24, 25, 26].

In this paper, we are interested in the following model:

${}_{0}{}\\!D^{\alpha}_{t}\phi=-\mbox{grad}_{H}E(\phi),$ (1.3)

deriving from gradient flows having a modified dissipation mechanism. Here
$0<\alpha<1$, ${}_{0}{}\\!D^{\alpha}_{t}$ is the Caputo fractional derivative
defined by

${}_{0}{}\\!D^{\alpha}_{t}\phi(t)=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}(t-s)^{-\alpha}\phi^{\prime}(s)ds.$

In mathematics gradient flows involving fractional derivatives have been
extensively studied in recent years; see, .e.g, [6, 42, 10, 1, 2, 3, 41, 29,
36, 27, 21]. From the definition it is seen that the fractional derivative is
some kind of weighted average in the history of the traditional derivative.
This means that the change rate, i.e., the derivative, at the current time is
affected by the historical rates. In a larger field this property has been
found quite useful in describing the memory effect which can be present, for
example, in some materials such as viscoelastic materials or polymers.
Intuitively, the gradient flows model (1.3) can be used to describe the
systems in which dissipation of the associated free energy has memory effect
in some circumstances. One of typical examples of such models is the
fractional Allen-Cahn equation, which is also the focus of this paper. There
exist a number of studies for this equation. Tang et al. [37] proved that the
time-fractional phase field model admits an integral type energy dissipation
law. They investigated the L1 time stepping scheme on the uniform mesh, which
is of first order energy stable accuracy. Du et al. [11] developed several
time schemes based on the convex splitting and weighted stabilization, and
proved that the convergence rates of their schemes are of order-$\alpha$ in
the uniform mesh without regularity assumption on the solution. Recently, Liao
et al. [22] proposed an adaptive second-order Crank-Nicolson time-stepping
scheme using SAV approach for the time-fractional MBE model, and showed that
the proposed scheme are unconditional stable on the nonuniform mesh. Very
recently, Quan et al. [31] theoretically proved the time-fractional energy law
and the weighted dissipation law. Accordingly, they constructed a first order
numerical method on the uniform time mesh [32], which preserved the energy
laws. However, it seems not easy to construct higher order schemes for
nonuniform meshes satisfying the same energy laws.

The aim of the present paper is to propose easy-to-implement and
unconditionally stable schemes, which preserve a non-local energy dissipation
law to be specified. The main idea in constructing the schemes is to use
existing efficient approximations to discretize the local part and history
part of the fractional derivative respectively, and use auxiliary variable
approaches [33, 34, 19] to deal with the nonlinear potential in the free
energy. The contributions of the paper are threefold:

  * •

Finding of a non-local energy dissipation law of the time fractional gradient
flows.

  * •

Construction of several unconditionally stable schemes for the time fractional
Allen-Cahn equation, which satisfy a discrete version of the energy
dissipation law. It is proved that the stability and energy dissipation law
remain true on the graded mesh, which is useful in recovering the optimal
convergence order for typical solutions having low regularity at the initial
time.

  * •

The proposed schemes are very easy to implement. That is, only several
Poisson-type equations with constant coefficients need to be solved at each
time step. Furthermore, a fast evaluation technique based on the sum-of-
exponentials approach is used to accelerate the calculation and reduce the
storage.

The paper is organized as follows: In the next section, we derive the non-
local energy dissipation law for the time fractional gradient flows. In
Section 3, we construct and analyze the first order numerical scheme. A
discrete energy dissipation law of the scheme is established for general time
girds. In Section 4, we propose and analyze two higher order schemes: a
$2-\alpha$ order and a second order schemes based on Crank-Nicolson formula.
The unconditional stability of the both schemes are rigorously proved. The
numerical experiments are carried out in Section 5, not only to validate
stability and accuracy of the proposed methods, but also to numerically
investigate the coarsening dynamics. Finally, the paper ends with some
concluding remarks.

##  2\. Non-local energy dissipation law

We consider the time fractional gradient flows (1.3) in the bounded domain
$\Omega\in\mathbb{R}^{n}~{}(n=1,2,3)$. When $\alpha=1$, it follows from
integrating (1.2) from $t_{e}$ to $t_{l}$ for any $0\leq t_{e}<t_{l}$ that:

$\displaystyle
E(\phi(t_{l}))-E(\phi(t_{e}))=-\int_{t_{e}}^{t_{l}}\Big{\|}\frac{\partial\phi(\cdot,s)}{\partial
s}\Big{\|}_{0}^{2}ds\leq 0.$

That is, the free energy $E(\cdot)$ is decreasing in $t$. However, for
$0<\alpha<1$, this energy law does not hold no longer. Instead, the solution
of (1.3) satisfies an non-local energy law that we derive below. To see that,
we split the fractional derivative into two parts as follows:

${}_{0}{}\\!D^{\alpha}_{t}\phi(t)=D^{\alpha,\hat{t}}_{l,t}\phi+D^{\alpha,\hat{t}}_{h,t}\phi,$
(2.1)

where $0<\hat{t}<t$, and the local term $D^{\alpha,\hat{t}}_{l,t}\phi$ and the
history term $D^{\alpha,\hat{t}}_{h,t}\phi$ are respectively defined by

$\displaystyle
D^{\alpha,\hat{t}}_{l,t}\phi=\frac{1}{\Gamma(1-\alpha)}\int_{\hat{t}}^{t}(t-s)^{-\alpha}\phi^{\prime}(s)ds,\quad\quad
D^{\alpha,\hat{t}}_{h,t}\phi=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{\hat{t}}(t-s)^{-\alpha}\phi^{\prime}(s)ds.$

Now we introduce the non-local “energy”

$\overline{E}(\hat{t},t;\phi):=E(\phi)+F_{h}(\hat{t},t;\phi),\ \ 0<\hat{t}<t,$
(2.2)

where $F_{h}(\hat{t},t;\phi)$ is the non-local part of the energy, defined by

$\displaystyle
F_{h}(\hat{t},t;\phi):=\int_{\Omega}\int_{\hat{t}}^{t}D^{\alpha,\hat{t}}_{h,s}\phi\frac{\partial\phi}{\partial
s}dsd{\tenbi x},\ \ 0<\hat{t}<t.$

A direct calculation shows

$\displaystyle\displaystyle\frac{d}{dt}\overline{E}(\hat{t},t;\phi)=\frac{d}{dt}{E}(\phi)+\Big{(}D^{\alpha,\hat{t}}_{h,t}\phi,\frac{\partial\phi}{\partial
t}\Big{)}=\Big{(}\mbox{grad}_{H}E(\phi),\frac{\partial\phi}{\partial
t}\Big{)}+\Big{(}D^{\alpha,\hat{t}}_{h,t}\phi,\frac{\partial\phi}{\partial
t}\Big{)},\ \ 0<\hat{t}<t.$

Then using (1.3) and (2.1) gives

$\displaystyle\displaystyle\frac{d}{dt}\overline{E}(\hat{t},t;\phi)=(-_{0}{}\\!D^{\alpha}_{t}\phi(t),\frac{\partial\phi}{\partial
t}\Big{)}+\Big{(}D^{\alpha,\hat{t}}_{h,t}\phi,\frac{\partial\phi}{\partial
t}\Big{)}=-\Big{(}D^{\alpha,\hat{t}}_{l,t}\phi,\frac{\partial\phi}{\partial
t}\Big{)},\ \ 0<\hat{t}<t.$ (2.3)

This allows to establish the following “energy” decay property: for
$0<\hat{t}<t$, integrating (2.3) yields

$\displaystyle\overline{E}(\hat{t},t;\phi(t))-\overline{E}(\hat{t},\hat{t};\phi(\hat{t}))=-\int_{\Omega}\int^{t}_{\hat{t}}D^{\alpha,\hat{t}}_{l,s}\phi\frac{\partial\phi}{\partial
s}dsd{\tenbi x}.$ (2.4)

Let’s define the bilinear form
$\displaystyle\mathcal{A}_{\alpha}^{\hat{t},t}(\cdot,\cdot)$ with
$0<\hat{t}<t$ by: for the functions $\varphi$ and $\psi$,

$\displaystyle\mathcal{A}_{\alpha}^{\hat{t},t}(\varphi,\psi):=\frac{1}{\Gamma(1-\alpha)}\int_{\hat{t}}^{t}\int_{\hat{t}}^{s}(s-\sigma)^{-\alpha}\varphi(\sigma)\psi(s)d\sigma
ds.$

It has been known; see, e.g., [37, 20], that the bilinear form
$\displaystyle\mathcal{A}_{\alpha}^{\hat{t},t}(\cdot,\cdot)$ is positive for
any $0<\hat{t}<t$. That is, for any $\psi\in L^{2}(\hat{t},t)$ so that the
following expression makes sense, it holds:

$\displaystyle\mathcal{A}_{\alpha}^{\hat{t},t}(\psi,\psi)\geq 0.$ (2.5)

It then follows from (2.4) and (2.5) that

$\displaystyle\overline{E}(\hat{t},t;\phi(t))-\overline{E}(\hat{t},\hat{t};\phi(\hat{t}))=-\int_{\Omega}\mathcal{A}_{\alpha}^{\hat{t},t}(\partial_{t}\phi,\partial_{t}\phi)d{\tenbi
x}\leq 0.$

This can be regarded as an energy law associated to the fucntional
$\overline{E}(\hat{t},t;\phi)$ defined in (2.2). However, the inconvenience in
using $\overline{E}(\hat{t},t;\phi)$ is that it depends on $\hat{t}$, which
makes the fucntional discontinuous in $t$. For the numerical purpose it is
desirable to derive a continuous-in-time energy, which is dissipative in a
given time grid, so that we have a clear goal to construct our numerical
scheme satisfying the same dissipation law. To this end, for a given time
mesh, say $0=t_{0}<t_{1}<t_{2}<\cdots<t_{M}=T$, we define a new non-local
“energy” functional as follows:

$\overline{E}(\phi)=\\!\\!\\!\begin{array}[]{l}\begin{cases}E(\phi),\
t\in[t_{0},t_{1}],\\\\[9.0pt] E(\phi)+F_{h}(t_{1},t;\phi),\
t\in[t_{1},t_{2}],\\\\[9.0pt] \displaystyle
E(\phi)+F_{h}(t_{n},t;\phi)+\sum_{k=1}^{n-1}F_{h}(t_{k},t_{k+1};\phi),\
t\in[t_{n},t_{n+1}],n=2,\cdots,M-1.\end{cases}\end{array}$ (2.6)

It is readily seen that

$\displaystyle\displaystyle\frac{d}{dt}\overline{E}(\phi)=-\Big{(}D^{\alpha,t_{n}}_{l,t}\phi,\frac{\partial\phi}{\partial
t}\Big{)},\ \ t\in[t_{n},t_{n+1}],n=0,1,\cdots,M-1.$

Integrating the above equality in the interval
$[t_{n},t_{n+1}],n=0,1,\cdots,M-1$ gives

$\overline{E}(\phi(t_{n+1}))-\overline{E}(\phi(t_{n}))=-\int_{\Omega}\mathcal{A}_{\alpha}^{t_{n},t_{n+1}}(\phi_{t},\phi_{t})d{\tenbi
x}\leq 0,\mbox{ for all }n=0,1,\cdots M-1.$ (2.7)

We see that for any given time mesh $\\{t_{n}\\}_{n=0}^{M}$, the corresponding
non-local free energy, defined in (2.6) is dissipative at the grid points.  We
would like to point out that the functionals defined in (2.2) and (2.6) do not
necessarily have any physical meaning. The motivation for introducing such a
modified “energy” is purely mathematical. That is, we want to find suitable
functionals related to the equation, which decay in time. This provides
insight into how a stable scheme should look like.  Our aim in the next
section is to construct numerical schemes that satisfy the same dissipation
law.

##  3\. Numerical methods — a first order scheme

To simplify the presentation, we will only consider the time fractional Allen-
Cahn equation, i.e.,

${}_{0}{}\\!D^{\alpha}_{t}\phi=-\mbox{grad}_{H}E(\phi),$ (3.1)

subject to the periodic boundary condition or Neumann boundary condition,
where the free energy $E(\cdot)$ is defined by

$E(\phi):=\displaystyle\int_{\Omega}\Big{[}\frac{\varepsilon^{2}}{2}|\nabla\phi|^{2}+F(\phi)\Big{]}d{\tenbi
x},$ (3.2)

and $F(\cdot)$ is a nonlinear potential. $H:=L^{2}(\Omega)$.

Auxiliary variable approach. The proposed schemes are based on a reformulation
of the time fractional Allen-Cahn equation by introducing an auxiliary
variable — an approach intensively studied recently for gradient flows; see,
e.g., [35, 19] and the references therein. The key to is to rewrite the
original equation (3.1)-(3.2) into the following equivalent form:

${}_{0}{}\\!D^{\alpha}_{t}\phi-\varepsilon^{2}\Delta\phi+\Big{(}1-\frac{R(t)}{R(t)}\Big{)}\theta^{2}\Delta\phi+\frac{R(t)}{R(t)}F^{\prime}(\phi)=0,$
(3.3)

where

$\displaystyle R(t)=\sqrt{\overline{E}_{\theta}(\phi)+C_{0}},\ \ \
\theta^{2}\leq\varepsilon^{2},$ (3.4)

and, for the time grid $\\{t_{n}\\}_{0}^{M}$, $\overline{E}_{\theta}(\phi)$ is
defined by

$\overline{E}_{\theta}(\phi)=\\!\\!\\!\begin{array}[]{l}\begin{cases}E_{\theta}(\phi),\
t\in[t_{0},t_{1}],\\\\[9.0pt] E_{\theta}(\phi)+F_{h}(t_{1},t;\phi),\
t\in[t_{1},t_{2}],\\\\[9.0pt] \displaystyle
E_{\theta}(\phi)+F_{h}(t_{n},t;\phi)+\sum_{k=1}^{n-1}F_{h}(t_{k},t_{k+1};\phi),\
t\in[t_{n},t_{n+1}],n=2,\cdots,M-1,\end{cases}\end{array}$

with
$E_{\theta}(\phi):=\int_{\Omega}[\frac{\theta^{2}}{2}|\nabla\phi|^{2}+F(\phi)]d{\tenbi
x}$, and $C_{0}$ being a constant such that
$\overline{E}_{\theta}(\phi)+C_{0}>0$.

To find a suitable way to discretize the auxiliary variable $R(t)$, we take
the derivative of (3.4) with respect to $t$ to obtain the auxiliary equation:

$\displaystyle\displaystyle\frac{dR}{dt}=\frac{1}{2R(t)}\Big{(}-\theta^{2}\Delta\phi+F^{\prime}(\phi)+D^{\alpha,t_{n}}_{h,t}\phi,\frac{\partial\phi}{\partial
t}\Big{)},\ \ \forall t\in[t_{n},t_{n+1}],n=2,\cdots,M-1.$ (3.5)

Furthermore, we use the operator splitting (2.1) to rewrite the equation (3.3)
under the equivalent form: for $n=2,\cdots,M-1$,

$\displaystyle
D^{\alpha,t_{n}}_{l,t}\phi-\varepsilon^{2}\Delta\phi+\Big{(}1-\frac{R(t)}{R(t)}\Big{)}\theta^{2}\Delta\phi+\frac{R(t)}{R(t)}(F^{\prime}(\phi)+D^{\alpha,t_{n}}_{h,t}\phi)=0,\
\ \forall t\in[t_{n},t_{n+1}].$ (3.6)

Now we are led to discretize the equations (3.6) and (3.5). The great
advantage of this approach is that, although we have one more variable and one
more equation to discretize compared to the original equation, constructing
stable and efficient schemes with help of the auxiliary variable turns out to
be a much easier task.

Before describing our schemes, we first realize, by taking
$L^{2}(\Omega)$-inner products and integrating from $t_{n}$ to $t_{n+1}$ of
(3.6) and (3.5) with $\frac{\partial\phi}{\partial t}$ and $2R(t)$
respectively, that

$\begin{array}[]{r@{}l}\displaystyle\Big{[}R^{2}(t_{n+1})+\frac{\varepsilon^{2}-\theta^{2}}{2}\|\nabla\phi(\cdot,t_{n+1})\|^{2}_{0}\Big{]}-\Big{[}R^{2}(t_{n})+\frac{\varepsilon^{2}-\theta^{2}}{2}\|\nabla\phi(\cdot,t_{n})\|^{2}_{0}\Big{]}&\displaystyle=-\int_{\Omega}\mathcal{A}_{\alpha}^{t_{n},t_{n+1}}(\phi_{t},\phi_{t})d{\tenbi
x}\\\ &\leq 0.\end{array}$ (3.7)

Noticing

$\displaystyle
R^{2}+\frac{\varepsilon^{2}-\theta^{2}}{2}\|\nabla\phi\|^{2}_{0}=\overline{E}(\phi)+C_{0},$

with $\overline{E}(\cdot)$ being the non-local energy functional defined in
(2.6), we have

$\displaystyle\overline{E}(\phi(t_{n+1}))-\overline{E}(\phi(t_{n}))\leq 0.$
(3.8)

This is exactly the same dissipation law as (2.7), derived from the original
equation without the auxiliary variable. We will see that after
discretization, the discrete solution satisfies a discrete dissipation law
under the form (3.7) rather than (3.8).

Starting with the equivalent equations (3.5) and (3.6), we are now in a
position to construct various efficient time stepping schemes to calculate the
solution $\phi$.

A first order scheme. Let $\tau_{n}=t_{n}-t_{n-1},n=1,\cdots,M$ be the time
step size, and $\tau=\max\\{\tau_{n},n=1,\cdots,M\\}$ be the maximum step
size.

We propose to use the popular L1 approximation [28] to discretize the local
and history parts of the Caputo fractional derivative at $t=t_{n+1}$:

$\begin{array}[]{r@{}l}D^{\alpha,t_{n}}_{l,t_{n+1}}\phi=&\displaystyle\frac{1}{\Gamma(1-\alpha)}\int_{t_{n}}^{t_{n+1}}(t_{n+1}-s)^{-\alpha}\phi_{s}(s)ds\\\\[9.0pt]
=&\displaystyle\frac{1}{\Gamma(1-\alpha)}\frac{\phi(t_{n+1})-\phi(t_{n})}{\tau_{n+1}}\int_{t_{n}}^{t_{n+1}}(t_{n+1}-s)^{-\alpha}ds+e^{n+1}_{l,\tau}\\\\[12.0pt]
:=&L^{\alpha}_{l}\phi(t_{n+1})+e^{n+1}_{l,\tau},\end{array}$
$\begin{array}[]{r@{}l}D^{\alpha,t_{n}}_{h,t_{n+1}}\phi=&\displaystyle\frac{1}{\Gamma(1-\alpha)}\sum_{k=0}^{n-1}\int_{t_{k}}^{t_{k+1}}(t_{n+1}-s)^{-\alpha}\phi_{s}(s)ds\\\\[9.0pt]
=&\displaystyle\sum_{k=0}^{n-1}\frac{1}{\Gamma(1-\alpha)}\frac{\phi(t_{k+1})-\phi(t_{k})}{\tau_{k+1}}\int_{t_{k}}^{t_{k+1}}(t_{n+1}-s)^{-\alpha}ds+e^{n+1}_{h,\tau}\\\\[12.0pt]
:=&L^{\alpha}_{h}\phi(t_{n+1})+e^{n+1}_{h,\tau},\end{array}$

where the discrete fractional operators $L^{\alpha}_{l}$ and $L^{\alpha}_{h}$
are defined respectively by

$\displaystyle
L^{\alpha}_{l}\phi(t_{n+1})=b_{0}\frac{\phi(t_{n+1})-\phi(t_{n})}{\tau_{n+1}},$
$\displaystyle
L^{\alpha}_{h}\phi(t_{n+1})=\sum_{k=0}^{n-1}b_{n-k}\frac{\phi(t_{k+1})-\phi(t_{k})}{\tau_{k+1}},$

and the coefficients $b_{k}$ are given by

$\displaystyle
b_{n-k}=\frac{1}{\Gamma(1-\alpha)}\int_{t_{k}}^{t_{k+1}}(t_{n+1}-s)^{-\alpha}ds>0,k=0,1,\cdots,n.$

The truncation errors $e^{n+1}_{l,\tau}$ and $e^{n+1}_{h,\tau}$ are defined
respectively by

$\displaystyle
e^{n+1}_{l,\tau}=\frac{1}{\Gamma(1-\alpha)}\Big{[}\int_{t_{n}}^{t_{n+1}}(t_{n+1}-s)^{-\alpha}\phi_{s}(s)ds-\displaystyle\frac{\phi(t_{n+1})-\phi(t_{n})}{\tau_{n+1}}\int_{t_{n}}^{t_{n+1}}(t_{n+1}-s)^{-\alpha}ds\Big{]},$

and

$\displaystyle
e^{n+1}_{h,\tau}=\frac{1}{\Gamma(1-\alpha)}\Big{[}\int_{0}^{t_{n}}(t_{n+1}-s)^{-\alpha}\phi_{s}(s)ds-\displaystyle\sum_{k=0}^{n-1}\frac{\phi(t_{k+1})-\phi(t_{k})}{\tau_{k+1}}\int_{t_{k}}^{t_{k+1}}(t_{n+1}-s)^{-\alpha}ds\Big{]}.$

For the graded mesh, i.e., $t_{n}=\big{(}{n\over M}\big{)}^{r}T,r\geq
1,n=0,1,\cdots,M$, which is particularly interesting for this problem and also
the focus of this paper, a direct calculation gives

$b_{j}=\frac{T^{1-\alpha}}{\Gamma(2-\alpha)M^{(1-\alpha)r}}\Big{[}\big{(}(n+1)^{r}-(n-j)^{r}\big{)}^{1-\alpha}-\big{(}(n+1)^{r}-(n-j+1)^{r}\big{)}^{1-\alpha}\Big{]}.$

Noting that when $r=1$, it is the uniform mesh.

It can be proved [28, 20] that the truncation error $e^{n+1}_{l,\tau}$ and
$e^{n+1}_{h,\tau}$ can be bounded by $c_{\phi}\tau^{2-\alpha}$, where
$c_{\phi}$ is a positive constant depending on the regularity of $\phi$.

The above approximation motivates us to consider the following scheme for
(3.6) and (3.5):

$\displaystyle\displaystyle
L^{\alpha}_{l}\phi^{n+1}-\varepsilon^{2}\Delta\phi^{n+1}+\Big{(}1-\frac{R^{n+1}}{R^{n}}\Big{)}\theta^{2}\Delta\phi^{n}+\frac{R^{n+1}}{R^{n}}\big{(}F^{\prime}(\phi^{n})+L^{\alpha}_{h}\phi^{n+1}\big{)}=0,$
(3.9a)
$\displaystyle\displaystyle\frac{R^{n+1}-R^{n}}{\tau_{n+1}}=\frac{1}{2R^{n}}\Big{(}-\theta^{2}\Delta\phi^{n}+F^{\prime}(\phi^{n})+L^{\alpha}_{h}\phi^{n+1},\frac{\phi^{n+1}-\phi^{n}}{\tau_{n+1}}\Big{)},$
(3.9b)

where $\phi^{n}$ is an approximation to $\phi(t_{n})$. Intuitively, this is a
first order scheme since it is a combination of some approximations of first
order and $2-\alpha$ order to different terms of the equations. However a
rigorous proof of the convergence order remains an open question. Instead, we
will provide a proof for the stability of the scheme, and the convergence
order will be verified through the numerical experiments to be presented later
on.

Stability. The stability property of the first order scheme (3.9) is presented
and proved in the following theorem.

######  Theorem 3.1.

Without any restriction on the mesh, the scheme (3.9) is stable in the sense
that the following discrete energy dissipation law holds

$\displaystyle\displaystyle
E^{n+1}_{\varepsilon,\theta}-E^{n}_{\varepsilon,\theta}\leq 0,\ \
n=0,1,\cdots,$

where

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

###### Proof.

First, taking the $L^{2}-$inner products of (3.9a) and (3.9b) with
$\phi^{n+1}-\phi^{n}$ and $2R^{n+1}$ respectively, we obtain:

$\begin{array}[]{r@{}l}&\displaystyle\big{(}L^{\alpha}_{l}\phi^{n+1},\phi^{n+1}-\phi^{n}\big{)}+\displaystyle(\varepsilon^{2}-\theta^{2})\big{(}\nabla\phi^{n+1},\nabla(\phi^{n+1}-\phi^{n})\big{)}+\theta^{2}\|\nabla(\phi^{n+1}-\phi^{n})\|^{2}_{0}\\\\[9.0pt]
&\hskip
150.79968pt\displaystyle+\frac{R^{n+1}}{R^{n}}\big{(}-\theta^{2}\Delta\phi^{n}+F^{\prime}(\phi^{n}+L^{\alpha}_{h}\phi^{n+1}),\phi^{n+1}-\phi^{n}\big{)}=0,\\\\[9.0pt]
&\displaystyle
2R^{n+1}(R^{n+1}-R^{n})=\displaystyle\frac{R^{n+1}}{R^{n}}\big{(}-\theta^{2}\Delta\phi^{n}+F^{\prime}(\phi^{n})+L^{\alpha}_{h}\phi^{n+1},\phi^{n+1}-\phi^{n}\big{)}.\end{array}$

Then combining these two equalities, using the identity

$\displaystyle
2a^{n+1}(a^{n+1}-a^{n})=|a^{n+1}|^{2}-|a^{n}|^{2}+|a^{n+1}-a^{n}|^{2},$

and dropping some non-essential positive terms, we obtain

$\displaystyle
E^{n+1}_{\varepsilon,\theta}-E^{n}_{\varepsilon,\theta}\leq-\big{(}L^{\alpha}_{l}\phi^{n+1},\phi^{n+1}-\phi^{n}\big{)}=-\frac{b_{0}}{\tau_{n+1}}\|\phi^{n+1}-\phi^{n}\|^{2}_{0}\leq
0,~{}~{}n=0,1,2,\cdots,$

This ends the proof. ∎

Implementation. Beside of its unconditional stability proved in Theorem 3.1,
another notable property of the scheme (3.9) is that it can be efficiently
implemented. To see that, we first eliminate $R^{n+1}$ from (3.9a) by using
(3.9b) to obtain

$\displaystyle
b_{0}\frac{\phi^{n+1}-\phi^{n}}{\tau_{n+1}}-\varepsilon^{2}\Delta\phi^{n+1}+\theta^{2}\Delta\phi^{n}+\Big{[}1+\frac{1}{2|R^{n}|^{2}}(\gamma^{n},\phi^{n+1}-\phi^{n})\Big{]}\gamma^{n}=0,$
(3.10)

where

$\gamma^{n}:=-\theta^{2}\Delta\phi^{n}+F^{\prime}(\phi^{n})+L^{\alpha}_{h}\phi^{n+1}.$

Then reformulating (3.10) gives

$\displaystyle\displaystyle\Big{(}\frac{b_{0}}{\tau_{n+1}}I_{d}-\varepsilon^{2}\Delta\Big{)}\phi^{n+1}+(\gamma^{n},\phi^{n+1})\frac{\gamma^{n}}{2|R^{n}|^{2}}=\frac{b_{0}}{\tau_{n+1}}\phi^{n}-\theta^{2}\Delta\phi^{n}\\!-\\!\Big{[}\frac{R^{n}}{R^{n}}\\!-\\!\frac{1}{2|R^{n}|^{2}}(\gamma^{n},\phi^{n})\Big{]}\gamma^{n}.$
(3.11)

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

$\displaystyle\begin{cases}\begin{array}[]{r@{}l}&\displaystyle\big{(}\frac{b_{0}}{\tau_{n+1}}I_{d}-\varepsilon^{2}\Delta\big{)}\phi^{n+1}_{1}=-\frac{\gamma^{n}}{2|R^{n}|^{2}},\\\\[9.0pt]
&\displaystyle\mbox{Neumann boundary condition or periodic boundary condition
on $\phi^{n+1}_{1}$;}\\\ \end{array}\end{cases}$ (3.12a)
$\displaystyle\begin{cases}\begin{array}[]{r@{}l}&\displaystyle\big{(}\frac{b_{0}}{\tau_{n+1}}I_{d}-\varepsilon^{2}\Delta\big{)}\phi^{n+1}_{2}=g(\phi^{n}),\\\\[9.0pt]
&\displaystyle\mbox{Neumann boundary condition or periodic boundary condition
on $\phi^{n+1}_{2}$;}\\\ \end{array}\end{cases}$ (3.12b)
$\displaystyle\displaystyle\phi^{n+1}=(\gamma^{n},\phi^{n+1})\phi^{n+1}_{1}+\phi^{n+1}_{2}.$
(3.12c)

In a first look it seems that (3.12c) governing the unknown $\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^{n+1}_{2}),$
(3.13)

where

$\displaystyle\sigma^{n}=-(\gamma^{n},\phi^{n+1}_{1})=\Big{(}\gamma^{n},A^{-1}\frac{\gamma^{n}}{2|R^{n}|^{2}}\Big{)}\
\ \mbox{with }A=\frac{b_{0}}{\tau_{n+1}}I_{d}-\varepsilon^{2}\Delta.$

Note that $A^{-1}$ is a positive definite operator. Thus $\sigma^{n}\geq 0$.
Then it follows from (3.13) that

$(\gamma^{n},\phi^{n+1})=\frac{(\gamma^{n},\phi^{n+1}_{2})}{1+\sigma^{n}}.$
(3.14)

Using this expression, $\phi^{n+1}$ can be explicitly computed from (3.12c).

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

(i) Calculation of $\phi^{n+1}_{1}$ and $\phi^{n+1}_{2}$: 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.14), then $\phi^{n+1}$
using (3.12c).

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

##  4\. Higher order schemes

###  4.1. A $2-\alpha$ order scheme

Using L1-CN formula [20] to discrete both the local and history parts of the
fractional derivative at $t=t_{n+\frac{1}{2}}:=\frac{t_{n}+t_{n+1}}{2}$ gives

$\begin{array}[]{r@{}l}D^{\alpha,t_{n}}_{l,t_{n+\frac{1}{2}}}\phi=&\displaystyle\frac{1}{\Gamma(1-\alpha)}\int_{t_{n}}^{t_{n+\frac{1}{2}}}(t_{n+\frac{1}{2}}-s)^{-\alpha}\phi_{s}(s)ds\\\\[9.0pt]
=&\displaystyle\frac{1}{\Gamma(1-\alpha)}\frac{\phi(t_{n+1})-\phi(t_{n})}{\tau_{n+1}}\int_{t_{n}}^{t_{n+\frac{1}{2}}}(t_{n+1}-s)^{-\alpha}ds+e^{n+\frac{1}{2}}_{l,\tau}\\\\[12.0pt]
:=&\widetilde{L}^{\alpha}_{l}\phi(t_{n+\frac{1}{2}})+e^{n+\frac{1}{2}}_{l,\tau},\end{array}$
$\begin{array}[]{r@{}l}D^{\alpha,t_{n}}_{h,t_{n+\frac{1}{2}}}\phi=&\displaystyle\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}(s)ds\\\\[9.0pt]
=&\displaystyle\sum_{k=0}^{n-1}\frac{1}{\Gamma(1-\alpha)}\frac{\phi(t_{k+1})-\phi(t_{k})}{\tau_{k+1}}\int_{t_{k}}^{t_{k+1}}(t_{n+\frac{1}{2}}-s)^{-\alpha}ds+e^{n+\frac{1}{2}}_{h,\tau}\\\\[12.0pt]
:=&\displaystyle\widetilde{L}^{\alpha}_{h}\phi(t_{n+\frac{1}{2}})+e^{n+\frac{1}{2}}_{h,\tau},\end{array}$

where the L1-CN difference operators $\widetilde{L}^{\alpha}_{l}$ and
$\widetilde{L}^{\alpha}_{h}$ are defined respectively by

$\widetilde{L}^{\alpha}_{l}\phi(t_{n+\frac{1}{2}})=\widetilde{b}_{0}\frac{\phi(t_{n+1})-\phi(t_{n})}{\tau_{n+1}},\quad\widetilde{L}^{\alpha}_{h}\phi(t_{n+\frac{1}{2}})=\sum_{k=0}^{n-1}\widetilde{b}_{n-k}\frac{\phi(t_{k+1})-\phi(t_{k})}{\tau_{k+1}},$

with

$\displaystyle\widetilde{b}_{0}=\frac{\tau_{n+1}^{1-\alpha}}{\Gamma(2-\alpha)2^{1-\alpha}}>0,~{}~{}~{}~{}~{}\widetilde{b}_{n-k}=\frac{1}{\Gamma(1-\alpha)}\int_{t_{k}}^{t_{k+1}}(t_{n+\frac{1}{2}}-s)^{-\alpha}ds>0,k=0,1,\cdots,n-1.$

The corresponding truncation errors $e^{n+1}_{l,\tau}$ and $e^{n+1}_{h,\tau}$
are defined respectively by

$\displaystyle
e^{n+\frac{1}{2}}_{\l,\tau}=\frac{1}{\Gamma(1-\alpha)}\Big{[}\int_{t_{n}}^{t_{n+\frac{1}{2}}}(t_{n+\frac{1}{2}}-s)^{-\alpha}\phi_{s}(s)ds-\displaystyle\frac{\phi(t_{n+1})-\phi(t_{n})}{\tau_{n+1}}\int_{t_{n}}^{t_{n+\frac{1}{2}}}(t_{n+\frac{1}{2}}-s)^{-\alpha}ds\Big{]},$

and

$\displaystyle
e^{n+\frac{1}{2}}_{h,\tau}=\frac{1}{\Gamma(1-\alpha)}\Big{[}\int_{0}^{t_{n}}(t_{n+\frac{1}{2}}-s)^{-\alpha}\phi_{s}(s)ds-\displaystyle\sum_{k=0}^{n-1}\frac{\phi(t_{k+1})-\phi(t_{k})}{\tau_{k+1}}\int_{t_{k}}^{t_{k+1}}(t_{n+\frac{1}{2}}-s)^{-\alpha}ds\Big{]}.$

It has been proved [20] that the truncation error $e^{n+1}_{l,\tau}$ and
$e^{n+1}_{h,\tau}$ are both of $2-\alpha$ order with respect to $\tau$.

Applying the difference operators $\widetilde{L}^{\alpha}_{l}$ and
$\widetilde{L}^{\alpha}_{h}$ to discretize the fractional derivative and the
Crank-Nicolson scheme to the remaining terms of the system (3.5) and (3.6), we
obtain L1-CN scheme as follows:

$\displaystyle\displaystyle\widetilde{L}^{\alpha}_{l}\phi^{n+\frac{1}{2}}-\varepsilon^{2}\frac{\Delta(\phi^{n+1}+\phi^{n})}{2}+\Big{(}1-\frac{R^{n+1}+R^{n}}{2R^{n+\frac{1}{2}}}\Big{)}\theta^{2}\Delta\phi^{n+\frac{1}{2}}$
$\displaystyle\hskip
199.16928pt+\frac{R^{n+1}+R^{n}}{2R^{n+\frac{1}{2}}}\Big{(}F^{\prime}(\phi^{n+\frac{1}{2}})+\widetilde{L}^{\alpha}_{h}\phi^{n+\frac{1}{2}}\Big{)},$
(4.1a)
$\displaystyle\displaystyle\frac{R^{n+1}-R^{n}}{\tau_{n+1}}=\frac{1}{2R^{n+\frac{1}{2}}}\Big{(}-\theta^{2}\Delta\phi^{n+\frac{1}{2}}+F^{\prime}(\phi^{n+\frac{1}{2}})+\widetilde{L}^{\alpha}_{h}\phi^{n+\frac{1}{2}},\frac{\phi^{n+1}-\phi^{n}}{\tau_{n+1}}\Big{)},$
(4.1b)

where
$\phi^{n+\frac{1}{2}}:=\phi^{n}+\frac{\tau_{n+1}}{2\tau_{n}}[\phi^{n}-\phi^{n-1}]$
and $R^{n+\frac{1}{2}}:=R^{n}+\frac{\tau_{n+1}}{2\tau_{n}}[R^{n}-R^{n-1}]$ are
the explicit approximation to $\phi(t_{n+\frac{1}{2}})$ and
$R(t_{n+\frac{1}{2}}),$ respectively.

The scheme (4.1) is expected to have $2-\alpha$ order convergence, since
formally the approximation to the fractional derivative
${}_{0}{}\\!D^{\alpha}_{t_{n+{1\over 2}}}\phi$ is of $2-\alpha$ order, and the
remaining approximations are based on the Crank-Nicolson formula, which is
second-order accurate.

The unconditional stability of the L1-CN scheme (4.1) is proved in the
following theorem.

######  Theorem 4.1.

The L1-CN scheme (4.1) satisfies the energy law:

$\displaystyle\displaystyle\frac{\varepsilon^{2}-\theta^{2}}{2}\|\nabla\phi^{n+1}\|^{2}_{0}+\frac{\theta^{2}}{4}\|\nabla(\phi^{n+1}-\phi^{n})\|_{0}^{2}+|R^{n+1}|^{2}$
$\displaystyle\displaystyle-\Big{[}\frac{\varepsilon^{2}-\theta^{2}}{2}\|\nabla\phi^{n}\|^{2}_{0}+\frac{\theta^{2}}{4}\Big{(}\frac{\tau_{n+1}}{\tau_{n}}\Big{)}^{2}\|\nabla(\phi^{n}-\phi^{n-1})\|_{0}^{2}+|R^{n}|^{2}\Big{]}\leq
0,~{}~{}n=0,1,\cdots.$

Therefore, the scheme (4.1) is unconditionally stable when i) the mesh is
uniform, i.e., $\frac{\tau_{n+1}}{\tau_{n}}=1$; ii) $\theta=0$. In the former
case, the discrete energy
$\frac{\varepsilon^{2}-\theta^{2}}{2}\|\nabla\phi^{n+1}\|^{2}_{0}+\frac{\theta^{2}}{4}\|\nabla(\phi^{n+1}-\phi^{n})\|_{0}^{2}+|R^{n+1}|^{2}$
decreases during the time stepping, while in the latter case, the energy
$\widetilde{E}^{n+1}_{\varepsilon}:=\frac{\varepsilon^{2}}{2}\|\nabla\phi^{n+1}\|^{2}_{0}+|R^{n+1}|^{2}$
decays in time.

###### Proof.

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

$\begin{array}[]{r@{}l}&\displaystyle(\widetilde{L}^{\alpha}_{l}\phi^{n+\frac{1}{2}},\phi^{n+1}-\phi^{n})+\frac{\varepsilon^{2}-\theta^{2}}{2}\big{(}\|\nabla\phi^{n+1}\|^{2}-\|\nabla\phi^{n}\|^{2}\big{)}\\\\[9.0pt]
&\hskip
113.81102pt\displaystyle+\frac{\theta^{2}}{2}\big{(}\nabla(\phi^{n+1}-\phi^{n}-\frac{\tau_{n+1}}{\tau_{n}}(\phi^{n}-\phi^{n-1})),\nabla(\phi^{n+1}-\phi^{n})\big{)}\\\\[9.0pt]
&\hskip
113.81102pt\displaystyle+\frac{R^{n+1}+R^{n}}{2R^{n+\frac{1}{2}}}\Big{(}-\theta^{2}\Delta\phi^{n+\frac{1}{2}}+F^{\prime}(\phi^{n+\frac{1}{2}})+\widetilde{L}^{\alpha}_{h}\phi^{n+\frac{1}{2}},\frac{\phi^{n+1}-\phi^{n}}{\tau_{n+1}}\Big{)}=0,\\\\[20.0pt]
&\displaystyle\frac{|R^{n+1}|^{2}-|R^{n}|^{2}}{\tau_{n+1}}=\frac{R^{n+1}+R^{n}}{2R^{n+\frac{1}{2}}}\Big{(}-\theta^{2}\Delta\phi^{n+\frac{1}{2}}+F^{\prime}(\phi^{n+\frac{1}{2}})+\widetilde{L}^{\alpha}_{h}\phi^{n+\frac{1}{2}},\frac{\phi^{n+1}-\phi^{n}}{\tau_{n+1}}\Big{)}.\end{array}$

Summing up the above two equalities, applying the identity

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

and dropping some non-essential positive terms, we obtain

$\displaystyle\displaystyle\frac{\varepsilon^{2}-\theta^{2}}{2}\|\nabla\phi^{n+1}\|^{2}_{0}+\frac{\theta^{2}}{4}\|\nabla(\phi^{n+1}-\phi^{n})\|_{0}^{2}+|R^{n+1}|^{2}$
$\displaystyle\displaystyle-\Big{[}\frac{\varepsilon^{2}-\theta^{2}}{2}\|\nabla\phi^{n}\|^{2}_{0}+\frac{\theta^{2}}{4}\Big{(}\frac{\tau_{n+1}}{\tau_{n}}\Big{)}^{2}\|\nabla(\phi^{n}-\phi^{n-1})\|_{0}^{2}+|R^{n}|^{2}\Big{]}$
$\displaystyle\leq-\big{(}\widetilde{L}^{\alpha}_{l}\phi^{n+\frac{1}{2}},\phi^{n+1}-\phi^{n})=-\frac{\widetilde{b}_{0}}{\tau_{n+1}}\|\phi^{n+1}-\phi^{n}\|^{2}_{0}\leq
0.$

This completes the proof. ∎

The L1-CN scheme (4.1) can be efficiently implemented by following the lines
similar to the first order scheme (3.9) described in the previous section.

###  4.2. A second order scheme

We first integrate the equation (3.6) and (3.5) from $t_{n}$ to $t_{n+1}$, and
multiply by $\frac{1}{\tau_{n+1}}$ to give:

$\displaystyle\displaystyle\frac{1}{\tau_{n+1}}\int_{t_{n}}^{t_{n+1}}D^{\alpha,t_{n}}_{l,t}\phi
dt-\frac{1}{\tau_{n+1}}\int_{t_{n}}^{t_{n+1}}\Big{[}\varepsilon^{2}\Delta\phi+\Big{(}1-\frac{R(t)}{R(t)}\Big{)}\theta^{2}\Delta\phi+\frac{R(t)}{R(t)}\big{(}F^{\prime}(\phi)+D^{\alpha,t_{n}}_{h,t}\phi\big{)}\Big{]}dt=0,$
$\displaystyle\displaystyle\frac{1}{\tau_{n+1}}\int_{t_{n}}^{t_{n+1}}\frac{dR}{dt}dt=\frac{1}{\tau_{n+1}}\int_{t_{n}}^{t_{n+1}}\frac{1}{2R(t)}\Big{(}-\theta^{2}\Delta\phi+F^{\prime}(\phi)+D^{\alpha,t_{n}}_{h,t}\phi,\frac{\partial\phi}{\partial
t}\Big{)}dt,\ n=0,1,\cdots.$ (4.2)

The idea is to construct second order approximations for the necessary terms
involved in the above equations. We define the finite difference operators
$\widehat{L}^{\alpha}_{l}$ and $\widehat{L}^{\alpha}_{h}$ by:

$\begin{array}[]{r@{}l}\widehat{L}^{\alpha}_{l}\phi(t_{n+\frac{1}{2}}):=&\displaystyle\frac{1}{\tau_{n+1}}\int_{t_{n}}^{t_{n+1}}\Big{[}\frac{1}{\Gamma(1-\alpha)}\int_{t_{n}}^{t}(t-s)^{-\alpha}\frac{\phi(t_{n+1})-\phi(t_{n})}{\tau_{n+1}}ds\Big{]}dt\\\\[12.0pt]
=&\displaystyle\widehat{b}_{0}\frac{\phi(t_{n+1})-\phi(t_{n})}{\tau_{n+1}},\end{array}$
$\begin{array}[]{r@{}l}\widehat{L}^{\alpha}_{h}\phi(t_{n+\frac{1}{2}})=&\displaystyle\frac{1}{\tau_{n+1})}\int_{t_{n}}^{t_{n+1}}\Big{[}\sum_{k=1}^{n}\frac{1}{\Gamma(1-\alpha)}\int_{t_{k}}^{t_{k+1}}(t_{n+\frac{1}{2}}-s)^{-\alpha}\frac{\phi(t_{k+1})-\phi(t_{k})}{\tau_{k+1}}ds\Big{]}dt\\\\[9.0pt]
=&\displaystyle\sum_{k=0}^{n-1}\widehat{b}_{n-k}\frac{\phi(t_{k+1})-\phi(t_{k})}{\tau_{k+1}},\end{array}$

where

$\displaystyle\widehat{b}_{0}=\frac{\tau_{n+1}^{1-\alpha}}{\Gamma(3-\alpha)},~{}~{}~{}~{}~{}\widehat{b}_{n-k}=\frac{1}{\Gamma(1-\alpha)\tau_{n+1}}\int_{t_{n}}^{t_{n+1}}\int_{t_{k}}^{t_{k+1}}(t-s)^{-\alpha}dsdt,k=0,1,\cdots,n-1.$

We want to use these two operators to approximate the local term
$\frac{1}{\tau_{n+1}}\int_{t_{n}}^{t_{n+1}}D^{\alpha,t_{n}}_{l,t}dt$ and the
history term
$\frac{1}{\tau_{n+1}}\int_{t_{n}}^{t_{n+1}}D^{\alpha,t_{n}}_{h,t}dt$
respectively. Note that a similar operator without splitting, called $L1^{+}$
formula, has been used in [22] to approximate the Caputo fractional
derivative, and the resulting scheme has been numerically found to be second
order accurate. However, there is no available analysis for the truncation
error, nor for the stability.

Applying $\widehat{L}^{\alpha}_{l}\phi(t_{n+\frac{1}{2}})$ and
$\widehat{L}^{\alpha}_{h}\phi(t_{n+\frac{1}{2}})$ to approximate
$\frac{1}{\tau_{n+1}}\int_{t_{n}}^{t_{n+1}}D^{\alpha,t_{n}}_{l,t}dt$ and
$\frac{1}{\tau_{n+1}}\int_{t_{n}}^{t_{n+1}}D^{\alpha,t_{n}}_{h,t}dt$
respectively, and trapezoidal formula to approximate the remaining terms in
(4.2), we arrive at the following $L1^{+}$-CN scheme:

$\displaystyle\displaystyle\widehat{L}^{\alpha}_{l}\phi^{n+\frac{1}{2}}-\varepsilon^{2}\frac{\Delta(\phi^{n+1}+\phi^{n})}{2}+\Big{(}1-\frac{R^{n+1}+R^{n}}{2R^{n+\frac{1}{2}}}\Big{)}\theta^{2}\Delta\phi^{n+\frac{1}{2}}$
$\displaystyle\hskip
199.16928pt+\frac{R^{n+1}+R^{n}}{2R^{n+\frac{1}{2}}}\Big{(}F^{\prime}(\phi^{n+\frac{1}{2}})+\widehat{L}^{\alpha}_{h}\phi^{n+\frac{1}{2}}\Big{)},$
(4.3a)
$\displaystyle\displaystyle\frac{R^{n+1}-R^{n}}{\tau_{n+1}}=\frac{1}{2R^{n+\frac{1}{2}}}\Big{(}-\theta^{2}\Delta\phi^{n+\frac{1}{2}}+F^{\prime}(\phi^{n+\frac{1}{2}})+\widehat{L}^{\alpha}_{h}\phi^{n+\frac{1}{2}},\frac{\phi^{n+1}-\phi^{n}}{\tau_{n+1}}\Big{)},$
(4.3b)

where
$\phi^{n+\frac{1}{2}}:=\phi^{n}+\frac{\tau_{n+1}}{2\tau_{n}}[\phi^{n}-\phi^{n-1}]$
and $R^{n+\frac{1}{2}}:=R^{n}+\frac{\tau_{n+1}}{2\tau_{n}}[R^{n}-R^{n-1}]$ are
the explicit approximation to $\phi(t_{n+\frac{1}{2}})$ and
$R(t_{n+\frac{1}{2}}),$ respectively.

We leave the error estimation as an open question, but present the stability
result for the scheme in the following theorem.

######  Theorem 4.2.

For the $L1^{+}$-CN scheme (4.3), it holds

$\displaystyle\displaystyle\frac{\varepsilon^{2}-\theta^{2}}{2}\|\nabla\phi^{n+1}\|^{2}_{0}+\frac{\theta^{2}}{4}\|\nabla(\phi^{n+1}-\phi^{n})\|_{0}^{2}+|R^{n+1}|^{2}$
$\displaystyle\displaystyle-\Big{[}\frac{\varepsilon^{2}-\theta^{2}}{2}\|\nabla\phi^{n}\|^{2}_{0}+\frac{\theta^{2}}{4}\Big{(}\frac{\tau_{n+1}}{\tau_{n}}\Big{)}^{2}\|\nabla(\phi^{n}-\phi^{n-1})\|_{0}^{2}+|R^{n}|^{2}\Big{]}$
$\displaystyle\displaystyle\leq-\frac{\widehat{b}_{0}}{\tau_{n+1}}\|\phi^{n+1}-\phi^{n}\|^{2}_{0}\leq
0,~{}~{}n=0,1,\cdots.$

This implies that the scheme is unconditionally stable in the cases of the
uniform mesh or $\theta=0$.

###### Proof.

The proof is very similar to Theorem 4.1, we leave it to interested readers. ∎

######  Remark 4.1.

One can also reformulate the original equation (3.1) into the following
equivalent form:

$\displaystyle\displaystyle
D^{\alpha,t_{n}}_{l,t}\phi-\varepsilon^{2}\Delta\phi+\Big{(}1-\frac{R(t)}{\sqrt{\overline{E}_{\theta}(\phi)+C_{0}}}\Big{)}\theta^{2}\Delta\phi+\frac{R(t)}{\sqrt{\overline{E}_{\theta}(\phi)+C_{0}}}\Big{(}F^{\prime}(\phi)+D^{\alpha,t_{n}}_{h,t}\phi\Big{)}=0,$
(4.4a)
$\displaystyle\displaystyle\frac{dR}{dt}=\frac{1}{2\sqrt{\overline{E}_{\theta}(\phi)+C_{0}}}\Big{(}-\theta^{2}\Delta\phi+F^{\prime}(\phi)+D^{\alpha,t_{n}}_{h,t}\phi,\frac{\partial\phi}{\partial
t}\Big{)},\ \ t\in[t_{n},t_{n+1}],n=0,1,\cdots.$ (4.4b)

Starting with this equivalent system and following the discussion in the above
sections, it is also possible to construct unconditionally stable schemes
based on L1, L1-CN and $L1^{+}-$CN formula for the discretization of the
fractional derivatives. However, compared to the schemes constructed for the
reformulation (3.5)-(3.6), a drawback using (4.4) is that one has to compute
$\overline{E}^{n}_{\theta}$ or $\overline{E}^{n+\frac{1}{2}}_{\theta}$, which
is more expensive than computing $R^{n}$ or $R^{n+\frac{1}{2}}$. Remember that
$\overline{E}^{n+\frac{1}{2}}_{\theta}$ is an explicit approximation to
$\overline{E}_{\theta}(\phi(t_{n+1/2}))$ involving the nonlocal terms
$F_{h}(t_{n},t_{n+\frac{1}{2}};\phi)+\sum_{k=1}^{n-1}F_{h}(t_{k},t_{k+1};\phi)$.
It is notable that the new SAV approach developed recently in [7, 8] may be
applied to deal with the nonlinear term and the history part of the fractional
derivative. Then energy stable schemes can be constructed based on the
following reformulation:

$\begin{array}[]{r@{}l}&\displaystyle
D^{\alpha,t_{n}}_{l,t}\phi-\varepsilon^{2}\Delta\phi+\big{(}1-\eta(t)\big{)}\theta^{2}\Delta\phi+\eta(t)\big{(}F^{\prime}(\phi)+D^{\alpha,t_{n}}_{h,t}\phi\big{)}=0,\\\\[9.0pt]
&\displaystyle\frac{d\overline{E}_{\theta}(\phi)}{dt}=\eta(t)\Big{(}-\theta^{2}\Delta\phi+F^{\prime}(\phi)+D^{\alpha,t_{n}}_{h,t}\phi,\frac{\partial\phi}{\partial
t}\Big{)},\end{array}$

where the scalar auxiliary function $\eta(t)$ is a Lagrange multiplier with
$\eta(0)=1.$

##  5\. Numerical results

This section is devoted to numerical investigation of the proposed schemes in
terms of the accuracy and stability. For the comparison purpose, we will
repeat most of the numerical examples in our previous work [20]. In the
following examples, we always set $\theta=0$ and $C_{0}=0$ in the schemes
unless specified otherwise. The spatial discretization is the Fourier method
or Legendre Galerkin spectral method using numerical quadratures. In order to
test the accuracy of the proposed schemes, the error is measured by the
maximum norm, i.e., $\displaystyle\max_{1\leq n\leq
M}\|\phi^{n}-\phi(t_{n})\|_{\infty}$ or $\displaystyle\max_{1\leq n\leq
M}\|\phi^{n}_{M}-\phi^{2n}_{2M}\|_{\infty}$, the latter norm will be used when
the exact solution is not available. In order to reduce the computational
complexity, a fast evaluation technique based on the sum-of-exponentials
approach [23, 38] is used to calculate the history part
$D^{\alpha,t_{n}}_{h,t}$ of the time fractional derivative.

###  5.1. Convergence order test

######  Example 5.1.

Consider the following fractional Allen-Cahn equation:

$\displaystyle{}_{0}{}\\!D^{\alpha}_{t}\phi-\varepsilon^{2}\Delta\phi-\phi(1-\phi^{2})=s,\quad({\tenbi
x},t)\in(0,2\pi)^{2}\times(0,T],$

subject to the periodic boundary condition, where $s({\tenbi x},t)$ is a
fabricated source term such that the exact solution is

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

The Fourier spectral method with $128\times 128$ modes is used to discretize
the equations in space. 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 as functions
of the time step sizes in log-log scale with $T=1$. It is observed, as
expected, that the L1 scheme (3.12), L1-CN scheme (4.1), and L$1^{+}$-CN
scheme (4.3) achieve respectively the first order, $2-\alpha$ order, and
second order convergence for all tested $\alpha$.

(a) first order scheme

(b) L1-CN scheme

(c) L$1^{+}$-CN scheme

Figure 1.  (Example 5.1) Error decay versus the time step sizes for the first
order scheme, L1-CN scheme, and L$1^{+}$-CN scheme with several different
$\alpha$.

######  Example 5.2.

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

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

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

In this test, a Legendre Galerkin spectral method with polynomials of degree
32 in each spatial direction is used for the spatial discretization. The
purpose of this test is to not only verify the convergence rate of the
schemes, but also investigate the impact of the regularity on the accuracy. In
particular, we are interested in studying the impact of the graded mesh
parameter $r$ on the convergence rate. It will help us to choose the optimal
value of $r$ to recover the convergence rate of the proposed schemes for low
regular solutions. The calculation is performed by using the L1-CN scheme
(4.1) and L$1^{+}$-CN scheme (4.3) with $M=64\times 2^{k},k=1,2,\cdots,9$. In
Figure 2, we plot the $L^{\infty}$ errors in log-log scale with respect to the
maximum time step size, i.e., $\tau=t_{M}-t_{M-1}$. The presented results are
in a perfect agreement with the expected convergence rates, i.e., $\min\\{\mu
r,2-\alpha\\}$ order for the L1-CN scheme, and $\min\\{\mu r,2\\}$ order for
the L$1^{+}$-CN scheme. This suggests use of the graded mesh with
$r=\frac{2-\alpha}{\mu}$ for the L1-CN scheme and $r=\frac{2}{\mu}$ for the
L$1^{+1}$-CN scheme. Doing so the schemes reach the optimal convergence rates
for solutions of this class.

(a) L1-CN scheme

(b) L$1^{+}$-CN scheme

Figure 2. (Example 5.2) Error history for the L1-CN scheme and L$1^{+}$-CN
scheme for different values of the mesh parameter $r$ and fractional
derivative order $\alpha$: the first line corresponds to L1-CN scheme; the
second line is for L$1^{+}$-CN scheme.

######  Example 5.3.

Consider the fractional Allen-Cahn Neumann problem in the domain
$(-1,1)\times(-1,1)$ with the initial condition $\phi({\tenbi x},0)=\cos(4\pi
x)\cos(4\pi y)$ without a source term. The exact solution is unavailable.

The spatial discretization uses the Legendre spectral method with high enough
mode number to avoid possible spatial error contamination. Since the exact
solution is unknown, the error of the numerical solution is defined as:

$\displaystyle e(\tau)=\displaystyle\max_{1\leq n\leq
M}\|\phi^{n}_{M}-\phi^{2n}_{2M}\|_{\infty}$

It is known that the time fractional operator creates some kind of singularity
at the initial time, and the solution behaves like $t^{\alpha}$ with respect
to the time variable; see, e.g., [11] for a numerical confirmation of this
behavior. In this example, we will take a closer look at the initial error to
investigate how the initial low regularity affects the accuracy of the
computed solutions. We also study the impact of the mesh ratio $r$ on the
convergence order. The calculation is run up $T=0.01$ using graded mesh with
grid points $M$ ranging from 20 to $10\times 2^{10}$. In Figure 3, we present
the errors in log-log scale with respect to the maximum mesh size, i.e.,
$\tau=\displaystyle\max_{1\leq n\leq M}(t_{n}-t_{n-1})=t_{M}-t_{M-1}$. It is
observed in the figure that the L1-CN scheme (4.1) and L$1^{+}$-CN scheme
(4.3) attain the convergence rate $\min\\{\alpha r,2-\alpha\\}$ and
$\min\\{\alpha r,2\\}$ respectively for all tested $\alpha$ and $r$. Clearly
the optimal $r$ is $\frac{2-\alpha}{\alpha}$ for the scheme (4.1) and
$\frac{2}{\alpha}$ for the scheme (4.3). In these cases, both schemes reach
the theoretical convergence order, i.e., $2-\alpha$ order and second order
respectively. For the uniform mesh, i.e., $r=1$, the schemes lose the optimal
convergence order. This is indicative that the regularity of the solution is
lower than what the theoretical convergence order of the schemes requests.

The stability of the proposed schemes is investigated through running the
calculation with $\varepsilon^{2}=0.001$ for long time, i.e., $T=50$, using a
time step size as large as possible. In view of the singularity feature at the
beginning time, we split the interval $[0,T]$ into two subintervals $[0,1]$
and $(1,T]$. We compute the solution using the graded mesh with optimal $r$ in
$[0,1]$, and using the uniform mesh with the time step size ${\triangle t}$ in
$(1,T]$. The computed modified energies and original energy are presented in
Figure 4. The modified energies shown in the left figures exhibits dissipative
features during the running time even if a large time step size ${\triangle
t}=1$ is used. This demonstrates good stability and modified energy
dissipation properties of the schemes proved in Theorems 4.1 and 4.2. However,
as shown in the right figures, the original energy is dissipative only for
relatively small time step size. Precisely, the original energy fails to keep
dissipation during some time period for the solution computed with ${\triangle
t}=1$. It is noteworthy that this failure is not due to the instability of the
schemes, but due to possible large error caused by the use of the large time
step size. Another notable fact is that all the original energy curves
coincide with each other for the time step sizes ranging from $0.0001$ to
$0.1$. In fact the computed original energy is a key indication of the
efficiency of the numerical methods for phase field models. The observed
dissipation feature of the original energy signifies that the solution evolves
in the right way, which is important for long time simulation. The last point
we want to emphasize is that although the time step size ${\triangle t}=1$ was
not able to produce correct original energy dissipation, use of larger
${\triangle t}$ is still possible through an adaptive strategy. That is,
adaptively utilize larger ${\triangle t}$ during the time period when the
phase transition is slow. Doing so will benefit the most from the
unconditional stability of the proposed schemes.

(a) L1-CN scheme

(b) L$1^{+}$-CN scheme

Figure 3. (Example 5.3) Effect of the graded mesh parameter $r$ on the
convergence rate of the L1-CN scheme (the first line) and L$1^{+}$-CN scheme
(the second line) with $\varepsilon^{2}=0.01$.

(a) L1-CN scheme with $\alpha=0.6$

(b) L$1^{+}$-CN scheme with $\alpha=0.6$

Figure 4. (Example 5.3) Evolution in time of the modified energy and original
energy computed by using different time step sizes.

###  5.2. Order sensibility of a benchmark problem

In this test case, we intend to apply the proposed schemes to the interface
moving problem governed by the fractional Allen-Cahn equation in the domain
$(-32,32)\times(-32,32)$. The initial state of the interface is the circle of
the radius $R_{0}=8$. It is known that the circle interface will shrink and
eventually disappear due to the driving force. This problem was also studied
in [20] to verify the performance of the scheme proposed in that paper. In the
classical case, i.e., $\alpha=1$, it was shown [39, 5] that the radius $R(t)$
of the circle at the given time $t$ evolves as $R(t)=\sqrt{R_{0}^{2}-2t}$.
That is, the ratio is monotonously decreasing and vanishes at $T=32$. This
problem has been frequently served as a benchmark to test the efficiency of
the numerical methods.

The problem is discretized by the L1-CN time stepping scheme using the graded
mesh with $r=\frac{2-\alpha}{\alpha}$ in the subinterval $[0,1]$ and the
uniform mesh in the subinterval $(1,T]$. The Legendre spectral method in space
uses polynomials of degree $128$ in each direction. For comparison purposes
the computation is performed with the same meshes for the same fractional
orders $\alpha$ as in [20]. The computed interface evolution is shown in
Figure 5. The computed $R^{2}$ and the free energy $E(\phi^{n})$ versus the
time are also plotted in Figure 6. The interface movement shown in Figure 5 is
almost the same compared to the results reported in [20]. The agreement on
$R^{2}$ and $E(\phi)$ between the current scheme and the one in [20] can be
observed equally. This demonstrates the efficiency of the both methods
proposed in [20] and in the present paper. However, as we have already
emphasized, the novelty of present work is the rigorous proof of an energy
dissipation law, both in the continuous and discrete cases, not only in the
uniform mesh but also in the graded mesh.

(a) $\alpha=1$ and ${\triangle t}=0.01$

(b) $\alpha=0.9,M=100$, and ${\triangle t}=0.01$

(c) $\alpha=0.4,M=1000$, and ${\triangle t}=0.01$

Figure 5. Snapshots of the interface evolution simulated by using the L1-CN
scheme with $C_{0}=1000$ for $\alpha=1,0.9$, and $0.4$.

(a) Interface radius as functions of the time.

(b) Free energy dissipation.

Figure 6. Comparison of the computed $R^{2}$ and total free energy computed
between the L1-CN scheme and the scheme in [20].

###  5.3. Coarsening dynamics

Finally we test the L$1^{+}$-CN scheme (4.3) to the two-phase coarsening
problem, by solving the fractional Allen-Cahn equation with
$\varepsilon^{2}=0.001$ in $(-1,1)^{2}$. The initial condition is a random
data, same as used in [20]. The simulation is performed in the graded mesh
with $r=\frac{2}{\alpha},M=100$ in $[0,1]$ and the uniform mesh with time step
size ${\triangle t}=0.01$ in $(1,T]$. The spatial spectral approximation uses
$128\times 128$ basis functions. Figure 7 presents some snapshots of the
simulated phase function and the computed free energy $E(\phi)$ versus time
for a number of fractional orders $\alpha=1,0.9,0.7$, and $0.5$. It is
observed that the results are not sensitive to the fractional orders at the
early stage as there is no distinguishable difference on the isoline plots
among different fractional orders before $t=5$. After that, the solutions
apparently start to deviate, and develop into long time phase coarsening. Note
that similar results and interpretation of these results have been presented
in [20]. The computed free energy $E(\phi^{n})$ shown in the last row figures
is also in a good agreement with the result reported in [20].

| $\alpha=1$ | $\alpha=0.9$ | $\alpha=0.7$ | $\alpha=0.5$
---|---|---|---|---
$t=0$ |  |  |  |
$t=5$ |  |  |  |
$t=20$ |  |  |  |
$t=50$ |  |  |  |
$t=100$ |  |  |  |
Energy |  |  |  |
Figure 7. Snapshots of the simulated phase field evolution starting with a
random initial data for $\alpha=1,0.9,0.7,0.5$ (the first five rows);
Comparisons of the computed energy between the L$1^{+}$-CN scheme and the
scheme in [20].

##  6\. Concluding remarks

We have proposed a class of energy dissipative schemes for the time fractional
Allen-Cahn equation. The construction of the schemes made use of a new idea to
reformulate the original equation. By splitting the time fractional derivative
into a local part and a history part, and adding the history part to the new
defined energy, a dissipation law for the new energy can be established in any
given time grid. Based on this splitting and an auxiliary variable approach,
several schemes of different convergence orders were constructed by combining
the L1 and L$1^{+}$ discretizations to the time fractional derivative and
Crank-Nicolson formula to other necessary terms. The main property of the
proposed schemes is its unconditional stability for general meshes. The proved
stability of the schemes built on the graded mesh is of particular interesting
because this type of mesh has been found very useful to treat with the
starting point singularity of the time fractional differential equations.
Moreover the splitting-based approach allows use of the sum-of-exponentials
techinique to fast evaluate the history part of the fractional derivative
without affecting the stability property of the schemes. The efficiency of the
proposed method was verified by a series of numerical experiments. the authors
proved that the fractional derivative of the traditional free energy is always
nonpositive.

It is notable that there exist some similar dissipation laws for memorized or
mean energy, and different energy laws lead to SAV-based schemes having quite
different stability properties. It seems to us that the energy dissipation law
established in the current paper facilitates construction of high order stable
schemes. It is also worth to mention that the idea of the present work is most
likely extendable to some other gradient flows, such as the Cahn-Hilliard
equation and molecular beam epitaxial growth models.

## 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] D. M. Anderson, G. B. Mcfadden, and A. A. Wheeler.  Diffuse-interface methods in fluid mechanics.  Annu. Rev. Fluid Mech., 30:139–165, 1997.
  * [5] 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.
  * [6] 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.
  * [7] Q. Cheng, C. Liu, and J. Shen.  A new Lagrange multiplier approach for gradient flows.  Comput. Methods Appl. Mech. Engrg., 367:113070, 2020.
  * [8] Q. Cheng and J. Shen.  Global constraints preserving scalar auxiliary variable schemes for gradient flows.  SIAM J. Sci. Comput., 42(4):A2489–A2513, 2020.
  * [9] M. Doi and S. F. Edwards.  THE THEORY OF POLYMER DYNAMICS.  Clarendon Press, 1986.
  * [10] 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.
  * [11] Q. Du, J. Yang, and Z. Zhou.  Time-fractional Allen-Cahn equations: analysis and numerical methods.  arXiv:1906.06584v1, pages 1–24, 2019.
  * [12] K. R. Elder and M. Grant.  Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals.  Phys. Rev. E, 70(5):051605, 2004.
  * [13] K. R. Elder, M. Katakowski, M. Haataja, and M. Grant.  Modeling elasticity in crystal growth.  Phys. Rev. Lett., 88(24):245701, 2002.
  * [14] K. R. Elder, N. Provatas, J. Berry, P. Stefanovic, and M. Grant.  Phase-field crystal modeling and classical density functional theory of freezing.  Phys. Rev. B, 75(6):794–802, 2007.
  * [15] J. Fraaije.  Dynamic density functional theory for microphase separation kinetics of block copolymer melts.  J. Chem. Phys., 99(11):9202–9212, 1993.
  * [16] J. Fraaije and G. Sevink.  Model for pattern formation in polymer surfactant nanodroplets.  Macromolecules, 36(21):7891–7893, 2003.
  * [17] L. Giacomelli and F. Otto.  Variatonal formulation for the lubrication approximation of the hele-shaw flow.  Calc. Var. Partial Dif., 2001.
  * [18] 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.
  * [19] 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.
  * [20] D. Hou, H. Zhu, and C. Xu.  Highly efficient schemes for time fractional Allen-Cahn equation using extended SAV approach.  arXiv:1910.09087, pages 1–28, 2019.
  * [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.  Commun. Nonlinear Sci. Numer. Simulat., 90:105366, 2020.
  * [22] B. Ji, H. Liao, Y. Gong, and L. Zhang.  Adaptive second-order Crank-Nicolson time-stepping schemes for time-fractional molecular beam epitaxial growth models.  SIAM J. Sci. Comput., 42(3):B738–B760, 2020.
  * [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] R. Larson.  Arrested tumbling in shearing flows of liquid-crystal polymers.  Macromolecules, 23(17):3983–3992, 1990.
  * [25] R. Larson and H. $\ddot{O}$ttinger.  Effect of molecular elasticity on out-of-plane orientations in shearing flows of liquid-crystalline polymers.  Macromolecules, 24(23):6270–6282, 1991.
  * [26] F. Leslie.  Theory of Flow Phenomena in Liquid Crystals.  Advances in Liquid Crystals, 4:1–81, 1979.
  * [27] 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.
  * [28] Y. Lin and C. Xu.  Finite difference/spectral approximations for the time-fractional diffusion equation.  J. Comput. Phys., 225(2):1533–1552, 2007.
  * [29] 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.
  * [30] F. Otto.  Lubrication approximation with prescribed nonzero contact anggle.  Commun. Part. Diff. Eq., 1998.
  * [31] C. Quan, T. Tang, and J. Yang.  How to define dissipation-preserving energy for time-fractional phase-field equations.  CSIAM Trans. Appl. Math., 1(3):478–490, 2020.
  * [32] C. Quan, T. Tang, and J. Yang.  Numerical energy dissipation for time-fractional phase-field equations.  arXiv:2009.06178, pages 1–22, 2020.
  * [33] J. Shen, J. Xu, and J. Yang.  The scalar auxiliary variable (SAV) approach for gradient flows.  J. Comput. Phys., 353:407–416, 2018.
  * [34] 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.
  * [35] Jie Shen, Jie Xu, and Jiang Yang.  The scalar auxiliary variable (SAV) approach for gradient flows.  Journal of Computational Physics, 353:407–416, 2018.
  * [36] 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.
  * [37] 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.
  * [38] 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.
  * [39] 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.
  * [40] 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.
  * [41] 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.
  * [42] 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.


# Exploring the Impact of Source Code Linearity on the Programmers’
Comprehension of API Code Examples

Seham Alharbi  0000-0002-1436-7992<EMAIL_ADDRESS>University of
YorkYorkUnited Kingdom and  Dimitris Kolovos  0000-0002-1724-6563
<EMAIL_ADDRESS>University of YorkYorkUnited Kingdom

(2024)

###### Abstract.

Context: Application Programming Interface (API) code examples are an
essential knowledge resource for learning APIs. However, a few user studies
have explored how the structural characteristics of the source code in code
examples impact their comprehensibility and reusability.
Objectives: We investigated whether the (a) linearity and (b) length of the
source code in API code examples affect users’ performance in terms of
correctness and time spent. We also collected subjective ratings.
Methods: We conducted an online controlled code comprehension experiment with
61 Java developers. As a case study, we used the API code examples from the
Joda-Time Java library. We had participants perform code comprehension and
reuse tasks on variants of the example with different lengths and degrees of
linearity.
Findings: Participants demonstrated faster reaction times when exposed to
linear code examples. However, no substantial differences in correctness or
subjective ratings were observed.
Implications: Our findings suggest that the linear presentation of a source
code may enhance initial example understanding and reusability. This, in turn,
may provide API developers with some insights into the effective structuring
of their API code examples. However, we highlight the need for further
investigation.

API comprehension, API code examples, source code linearity, human factors in
software engineering, controlled experiment

††journalyear: 2024††copyright: acmlicensed††conference: 32nd IEEE/ACM
International Conference on Program Comprehension; April 15–16, 2024; Lisbon,
Portugal††booktitle: 32nd IEEE/ACM International Conference on Program
Comprehension (ICPC ’24), April 15–16, 2024, Lisbon, Portugal††doi:
10.1145/3643916.3644395††isbn: 979-8-4007-0586-1/24/04††ccs: Software and its
engineering Documentation††ccs: Software and its engineering Programming by
example††ccs: Human-centered computing Empirical studies in HCI

##  1\. Introduction

Working through Application Programming Interface (API) code examples has been
proven to be the most preferred learning strategy for both beginner and
experienced API users (Meng et al., 2018). Surprisingly, little is known about
how the different source code structures in these examples affect their
comprehensibility and reusability. Furthermore, existing work (discussed in
Section 5) appears to be focused on examining the impact of several source
code characteristics on comprehension only in the context of generic software.
Therefore, to fill this gap, the study presented in this paper focuses
specifically on API code examples. Moreover, unlike existing studies, we
examine different source code constructs that illustrate the same API usage
and functionality, narrowing the evaluation focus to the constructs’ impact on
comprehension. We also assess an additional concept, which is the impact of
the examined source code structures on the reusability of API code examples.

In this study, we are particularly interested in exploring the impact of two
source code aspects: the degree of linearity and length. Linear source code
refers to code that can be read primarily in a sequential order without
interference from interdependent methods or classes. Considering the absence
of jumps between method definitions in such a code, we hypothesise it may be
easier to comprehend. In addition, linear API code examples may be easier to
reuse and adapt into one’s codebase, as they typically contain a single self-
contained method that can be copied and edited, as opposed to non-linear code
examples that involve multiple methods.

Through our study, we aim to help API developers understand how to structure
their API code examples more effectively, thereby enhancing the examples’
comprehensibility and reusability. This, in turn, would promote the
learnability of their APIs.

This work is a part of larger research (Alharbi et al., 2022) in which we
develop tools and techniques to enhance the maintainability and
comprehensibility of API code examples. All data collected or used in this
study is available in our replication package.111Replication package:
https://figshare.com/s/52e11ece2f39bac64bcb

##  2\. Methodology

###  2.1. Research Questions

We intended to answer the following research questions:

RQ1::

In terms of correctness and time spent, how does the linearity of an API code
example impact a programmer’s performance in tasks that require code
comprehension?

RQ2::

What effects does the length of a linear API code example222Please refer to
our (replication package $\rightarrow$ examples) for some sample code
illustrating linear and non-linear API code examples. have on its
comprehensibility and reusability?

RQ3::

Does the degree of linearity in a non-linear API example affect its
comprehensibility and reusability?

###  2.2. Study Design

This study was conducted online to allow access to a large and more diverse
pool of participants. We utilised Gorilla (gor, 2023), a widely used online
experiment builder, which provided all the features we needed in our study
(randomisation, counterbalancing presentation of formatted source code,
accurate reaction times and integration with participant recruitment
platforms). We recruited Java developers from
Prolific333https://www.prolific.com – a participant recruitment platform for
online research – and assigned them to two main groups: linear vs non-linear
(between-subjects). Each group consisted of two sub-groups that correspond to
the treatment categories shown in Table 1. Participant assignment to groups
was fully randomised and balanced, with a 1:1 ratio. Each participant
completed two code comprehension and reuse tasks from the same treatment
category (highlighted in the same colour in Table 1). The order in which each
participant received tasks was also randomised to eliminate any potential
order effects. Ethical approval was obtained before the study was conducted.

###  2.3. Independent Variables

We considered a single independent variable: the source code structure of API
code examples. Two primary source code factors were systematically varied:
code linearity, which is manipulated using the source code linearity metric
(i) proposed by Peitek et al. (Peitek et al., 2020); and code length, which is
varied by adjusting the number of lines of code (LOC). As shown in Table 1, we
combined these two factors and generated four treatment categories: linear-
short, linear-long, and non-linear with varying levels of linearity $(i)$,
ranging from $(10.00<i\leq 15.00)$ to $(15.00<i\leq 20.00)$. These selected
values reflect a diverse spectrum of code linearity.

###  2.4. Dependent Variables

We measured three dependent variables: reaction time, correctness and
subjective rating.

Time Duration Marking. For the comprehension phase, reaction time was defined
as the amount of time that elapsed between a participant’s initial view of an
API code example and submission of their overall comprehension rating
Similarly, for the code-reuse phase, reaction time was defined as the time
between a participant’s first view of the required code-reuse task and the
submission of their solution code.

Judging Correctness. We ensured marking consistency by defining a set of
correctness categories and criteria444Detailed criteria are available in our
replication package.: correct (A), almost correct (B), partially correct (C),
incorrect (D) and absent (F).

Table 1. API code examples and their variants, metric values and study results. Variants sharing the same colour belong to the same treatment category. The imbalance in the number of responses (N) between two variants of the same example is due to our exclusion of responses with inaccurately reported break times. API Code Example | Variant | Metrics | N | Correctness | Comprehension | Reuse
---|---|---|---|---|---|---
|  | LOC | Complexity |  Linearity (i) |  # Method Calls |  |  |  Median Reaction Time |  p-value∗ |  Median Reaction Time |  p-value∗
Date Example | Linear | 31 | 7 | 0.00 | 0 | 17 | 15 (88%) | 44s | 0.01 | 6m 11s | 0.36
| Non-linear | 14 | 1 | 12.95 | 5 | 13 | 8 (62%) | 1m 38s | 7m 14s
Chronology Example | Linear | 25 | 4 | 0.00 | 0 | 18 | 10 (56%) | 43s | 0.03 | 3m 37s | 0.10
| Non-linear | 11 | 1 | 11.41 | 4 | 14 | 13 (93%) | 1m 49s | 6m 15s
Duration Example | Linear | 35 | 7 | 0.00 | 0 | 16 | 9 (56%) | 3m 31s | 0.60 | 16m 18s | 0.40
| Non-linear | 21 | 2 | 19.47 | 6 | 12 | 5 (42%) | 2m 17s | 24m 59s
Interval Example | Linear | 46 | 7 | 0.00 | 0 | 14 | 10 (71%) | 1m 14s | 0.33 | 9m 21s | 0.06
| Non-linear | 33 | 5 | 19.43 | 6 | 12 | 9 (75%) | 2m 49s | 12m 56s
Overall | Linear |  |  |  | 65 | 44 (68%) | 58s |  | 7m 17s |
| Non-linear |  |  |  | 51 | 35 (69%) | 1m 50s |  | 10m 13s |
∗ Mann–Whitney U test |  |  |  |  |  |  |  |  |  |  |

###  2.5. Participants

Pilot. To validate the study design, we conducted a pilot with four
participants (average age 28.5 ± 9.5; average years of Java programming 3.2 ±
1.5). Based on the results, we reduced the number of tasks assigned to each
participant from four to two to minimise the experiment’s overall duration. We
also changed the online IDE555https://replit.com used due to its slow
execution time and enhanced the wording of the code optimisation question. The
data collected in the pilot study was not used in the final analysis.

Pre-screening survey. In addition to the pre-screeners provided by Prolific,
we created a separate programming knowledge survey to assess participants’
programming knowledge before they participated in the study. This was
essential since recent research revealed that, while recruitment platforms,
such as Prolific, greatly mitigate self-selection bias and some security
issues, their pre-screeners may not always be reliable (Danilova et al., 2021;
Russo, 2022; Reid et al., 2022). In this survey, we used the basic knowledge
questions and time limit recommended by Danilova et al. (Danilova et al.,
2021). Participants who correctly answered all the questions within the time
limit were manually invited to participate in the study.

Study. We recruited 61 Java developers from 14 countries with varying levels
of programming experience. Only 19 (31%) participants stated that they
previously used the Joda-Time Java library. Among them, only eight (42%) said
that they used it more than once, and none of them reported regular usage. The
participants’ prior programming experience was assessed using a validated
questionnaire that is based on self-estimates (Feigenspan et al., 2012;
Siegmund et al., 2014). Each participant was compensated £10 for their time
and effort. Additional information about the participant demographics is shown
in Table 2.

Table 2. Participant demographics. Category |  n=61
---|---
Student |  43 (70%)
Professional Developer |  18 (30%)
Programming Experience (in Years) |  5.9 ± 3.3
Java Programming Experience (in Years) |  3.2 ± 2.2
Familiarity with Joda-Time |  19 (31%)
Male |  54 (89%)
Female |  7 (11%)
Age (in Years) |  25.4 ± 6.1

###  2.6. Material

API code examples. We chose the Joda-Time666https://www.joda.org/joda-time/
Java library because it addresses a well-known concept (i.e. date and time
handling). Joda-Time met our selection criteria of: 1) not requiring prior
domain knowledge that would pose an unnecessary challenge to participants; 2)
being well-documented and 3) not being too popular so that an average Java
developer would not necessarily be familiar with it. Furthermore, we intended
to utilise the code examples available on the Joda-Time documentation page.6
However, we found that these examples were not complex enough. Thus, we
decided to create our own examples.

As shown in Table 1, we developed four API code examples, each of which
demonstrated a distinct usage of Joda-Time. This variation of examples was
important to minimise the risk of any potential learning effect arising from
within the examples. We then created linear and non-linear versions of each
example. We strove to make these examples look as natural as possible by: 1)
properly documenting them and 2) letting them depict real-world scenarios such
as data manipulation or meeting scheduling. For the non-linear versions of the
examples, we refactored the linear version by extracting some functionalities
into a set of utility methods and replacing the extracted code with method
calls using the extract method refactoring technique. Examples within the same
treatment category were relatively comparable in terms of length, complexity,
degree of linearity and the number of utility method calls.

Tasks. Oftentimes, when API users turn to code examples to learn a new API,
they typically have a specific problem in mind. They are hoping that the code
in the example they are reviewing will be reusable. If this is possible, they
copy and paste the example, then modify its source code by adding or deleting
statements to match their needs (Meng et al., 2018; Gao et al., 2020). In our
study, we wanted to simulate this behaviour. Therefore, each code-reuse task
had two parts: 1) code modification, in which participants were required to
make changes to address a specific problem; and 2) code optimisation, in which
we asked them to remove any unnecessary code that did not directly contribute
to their task solution. The tasks were generally easy and designed to be
solved with a few edits. Each API code example had a unique task that remained
the same for both versions of the example.

###  2.7. Experiment Procedure

After obtaining their consent, we asked the participants to complete a
demographics questionnaire. Subsequently, each of them was randomly assigned
two code examples from the same treatment category. This means that each
participant completed two distinct code-reuse tasks.

We designed each task to be completed in four sequential parts. The first was
the comprehension part, in which participants were asked to review the example
and rate their own understanding. The next part pertained to instructions, in
which participants were given a link to an online
IDE777https://www.jdoodle.com that contained a Java project of the example,
with Joda-Time imported and ready to use. Participants were also instructed on
how to download the example if they preferred using their own IDE. The third
part was the code-reuse task, in which participants answered a two-part
question (as explained in Section 2.6) and pasted their solution code in a
given text box. The final part involved post-task questions, in which
participants were asked to rate how difficult it was to reuse the code,
whether they employed the provided online IDE or their own, and report any
break time (if any was taken).

We only measured the time spent on two of the four parts: 1) comprehension
time and 3) reuse time. The rationale for separating the comprehension and
reuse of the same code example was to reduce participants’ use of the ‘as-
needed’ program comprehension strategy (Littman et al., 1987; Von Mayrhauser
and Vans, 1998; Feitelson, 2021).

###  2.8. Data Analysis

We manually analysed the correctness of responses for each task. First, we
converted the categories mentioned in Section 2.4 to numerical values (correct
(A) = 100%, almost correct (B) = 70%, partially correct (C) = 40%; both
incorrect (D) and absent (F) = 0%). We applied the same scale to both parts of
the code-reuse task (i.e. code modification and code optimisation). However,
when calculating the overall task score, we assigned more weight (90%) to the
first part of the task, as it required greater effort than the second part,
which accounted for only 10% of the total weight. Responses with an overall
score of 60% or higher were considered correct. This overall correctness
threshold ensures that participants achieve at least 70% in the first part of
the task.

When analysing reaction times automatically captured by Gorilla (gor, 2023),
we only considered correct responses. We used the Shapiro-Wilk test (Shapiro
and Wilk, 1965) to assess the normality of reaction times and correctness as
well as Levene’s test (Howard, 1960) to evaluate variance homogeneity. The
findings indicated non-normality and unequal variances for both correctness
and reaction times. Therefore, to test for statistically significant
differences, we used a non-parametric test, the Mann–Whitney U test (Wilcoxon
rank-sum test), with a significance level of $\alpha=0.05$.

##  3\. Results and Discussion

As shown in Table 1 and Figure 1, participants generally spent less time
comprehending and reusing linear code examples (both in mean and median
reaction times). This observation suggests that the source code linearity in
an API code example may affect a programmer’s performance. This influence has
a greater impact on comprehension and is statistically significant when the
linear API code example is also short (e.g. date and chronology examples).
Moreover, in terms of reusability, there appears to be a trend towards
significance in two of the API code examples (chronology and interval
examples), as reflected by their moderate p-values of $0.10$ and $0.06$,
respectively. However, the impact on correctness (Mann–Whitney U test,
$W=1690,p=0.428$) and subjective rating (as shown in Figure 2) was not
substantial (RQ1).

The Mann–Whitney U test revealed a significant difference in both
comprehension ($W=125,p=0.004$) and reusability ($W=99,p=0.000$) between the
groups that received linear-short and linear-long API code examples (RQ2).
Similarly, participants spent less time reusing the non-linear code examples
when the linearity value (i) was lower ($i<15.00$, Mann–Whitney U test:
$W=68,p=0.003$). Notably, unlike the comparison in RQ1, this comparison is
based on code examples illustrating different API usage; thus, the significant
differences in participants’ performance could be due to variations in the
implemented API functionality and required tasks (RQ3).

(a)

(b)

Figure 1. The time spent on (a) comprehending, and (b) reusing the API code
examples used in the study. Each boxplot represents the responses for one
version (linear or non-linear) of a single example.

(a)

(b)

Figure 2. Participants’ subjective ratings of API code example comprehension
(a) and reusability difficulty (b).

##  4\. Threats to Validity

Construct Validity. To mitigate construct validity threats, we used specific
metrics to manipulate our independent variables. We created code examples and
tasks that represented real-world scenarios and ensured a consistent
correctness evaluation. However, since the experiment was conducted online, we
had limited control over participants’ activities. To address this, we asked
them to self-report break times, which were subtracted from the time spent on
solving the task, and disclose whether they used their own IDEs. While this
approach provided insights into participants’ behaviour, it was not entirely
conclusive.

Internal Validity. We reduced internal validity threats by randomly assigning
participants to treatment groups and randomising the order in which they
viewed tasks. We also administered a validated programming experience
questionnaire and pre-screened participants for their programming knowledge.

External Validity. One potential threat to external validity was the study’s
limited scope. It focused solely on the API of one Java library and included
only a few code examples. Also, 70% of the participants were students, which
limited the generalisability of our findings. The study’s virtual setting and
the use of an online IDE, which lacks features such as auto-completion and
error checking, may not fully reflect the conditions of a traditional coding
environment.

##  5\. Related Work

A large number of studies investigated the impact of various source code
characteristics on programmers’ code comprehension. These characteristics
include the use of intermediate variables (Cates et al., 2021), certain
syntactic structures such as ifs and for loops (Ajami et al., 2019), different
identifier names (Hofmeister et al., 2019; Avidan and Feitelson, 2017) as well
as naming conventions (Sharif and Maletic, 2010; Binkley et al., 2013).
Moreover, some studies explored the effects of more global factors such as the
order of methods (Geffen and Maoz, 2016), code regularity (Jbara and
Feitelson, 2014), and the linearity of source code and reading order (Peitek
et al., 2020).

##  6\. Conclusion and Future Work

In this paper, we investigated the impact of source code linearity and length
on the comprehensibility and reusability of API code examples. We chose code
examples from the Joda-Time Java library and manipulated their structure.
Furthermore, we recruited 61 Java developers, assigned each one of them two
code examples, and asked them to complete code-reuse tasks. This study found
that participants demonstrated relatively faster reaction times when working
with linear API examples.

For future work, we intend to expand this study by incorporating a broader
range of APIs from diverse domains. This will involve utilising a larger set
of code examples with varying levels of linearity and increasing the number of
participants. Also, to better capture participants’ activities, we plan to
conduct this experiment in a laboratory setting. Additionally, we are
interested in determining whether the activities of participants, as they work
with code examples of different linearity, still align with the activities
reported in existing studies on the COIL888Collection and Organization of
Information for Learning. model (Gao et al., 2020; Sparman and Schulte, 2023).

## References

  * (1)
  * gor (2023) 2023\.  Gorilla Experiment Builder.  https://gorilla.sc
  * Ajami et al. (2019) Shulamyt Ajami, Yonatan Woodbridge, and Dror G. Feitelson. 2019.  Syntax, predicates, idioms — what really affects code complexity?  _Empirical Software Engineering_ 24, 1 (2 2019), 287–328.  https://doi.org/10.1007/s10664-018-9628-3
  * Alharbi et al. (2022) Seham Alharbi, Dimitris Kolovos, and Nicholas Matragkas. 2022.  Synthesising Linear API Usage Examples for API Documentation. In _2022 IEEE International Conference on Software Maintenance and Evolution (ICSME)_. 607–611.
  * Avidan and Feitelson (2017) Eran Avidan and Dror G. Feitelson. 2017.  Effects of Variable Names on Comprehension: An Empirical Study. In _IEEE 25th International Conference on Program Comprehension_. IEEE Computer Society, 55–65.  https://doi.org/10.1109/ICPC.2017.27
  * Binkley et al. (2013) Dave Binkley, Marcia Davis, Dawn Lawrie, Jonathan I. Maletic, Christopher Morrell, and Bonita Sharif. 2013.  The impact of identifier style on effort and comprehension.  _Empirical Software Engineering_ 18, 2 (4 2013), 219–276.  https://doi.org/10.1007/s10664-012-9201-4
  * Cates et al. (2021) Roee Cates, Nadav Yunik, and Dror G. Feitelson. 2021.  Does Code Structure Affect Comprehension? On Using and Naming Intermediate Variables. In _IEEE International Conference on Program Comprehension_ , Vol. 2021-May. IEEE Computer Society, 118–126.  https://doi.org/10.1109/ICPC52881.2021.00020
  * Danilova et al. (2021) Anastasia Danilova, Alena Naiakshina, Stefan Horstmann, and Matthew Smith. 2021.  Do you really code? Designing and evaluating screening questions for online surveys with programmers. In _Proceedings - International Conference on Software Engineering_. IEEE Computer Society, 537–548.
  * Feigenspan et al. (2012) Janet Feigenspan, Christian K¨astner, J¨org Liebig, Sven Apel, and Stefan Hanenberg. 2012.  Measuring Programming Experience. In _2012 20th IEEE International Conference on Program Comprehension_. IEEE, 73–82.  https://doi.org/10.1109/ICPC.2012.6240511
  * Feitelson (2021) Dror G. Feitelson. 2021.  Considerations and Pitfalls in Controlled Experiments on Code Comprehension. In _IEEE/ACM 29th International Conference on Program Comprehension (ICPC)_. 106–117.  https://doi.org/10.1109/ICPC52881.2021.00019
  * Gao et al. (2020) Gao Gao, Finn Vpichick, Michelle Ichinco, and Caitlin Kelleher. 2020.  Exploring Programmers’ API Learning Processes: Collecting Web Resources as External Memory. In _2020 IEEE Symposium on Visual Languages and Human-Centric Computing (VL/HCC)_. 1–10.
  * Geffen and Maoz (2016) Yorai Geffen and Shahar Maoz. 2016.  On method ordering. In _IEEE International Conference on Program Comprehension_ , Vol. 2016-July. IEEE Computer Society, 1–10.  https://doi.org/10.1109/ICPC.2016.7503711
  * Hofmeister et al. (2019) Johannes C. Hofmeister, Janet Siegmund, and Daniel V. Holt. 2019.  Shorter identifier names take longer to comprehend.  _Empirical Software Engineering_ 24, 1 (2 2019), 417–443.  https://doi.org/10.1007/s10664-018-9621-x
  * Howard (1960) Levene Howard. 1960.  Robust tests for equality of variances.  _Contributions to probability and statistics_ (1960), 278–292.
  * Jbara and Feitelson (2014) Ahmad Jbara and Dror G. Feitelson. 2014.  On the effect of code regularity on comprehension. In _22nd International Conference on Program Comprehension (ICPC 2014)_. Association for Computing Machinery, 189–200.  https://doi.org/10.1145/2597008.2597140
  * Littman et al. (1987) David Littman, Jeannine Pinto, Stanley Letovsky, and Elliot Soloway. 1987.  Mental models and software maintenance.  _Journal of Systems and Software_ 7, 4 (1987), 341–355.  https://www.sciencedirect.com/science/article/pii/0164121287900331
  * Meng et al. (2018) Michael Meng, Stephanie Steinhardt, and Andreas Schubert. 2018.  Application Programming Interface Documentation: What Do Software Developers Want?  _Journal of Technical Writing and Communication_ 48, 3 (2018), 295–330.  https://doi.org/10.1177/0047281617721853
  * Peitek et al. (2020) Norman Peitek, Janet Siegmund, and Sven Apel. 2020.  What drives the reading order of programmers? an eye tracking study. In _28th International Conference on Program Comprehension (ICPC ’20)_. ACM, 342–353.  https://doi.org/10.1145/3387904.3389279
  * Reid et al. (2022) Brittany Reid, Markus Wagner, Marcelo d’Amorim, and Christoph Treude. 2022.  Software Engineering User Study Recruitment on Prolific: An Experience Report. In _International Workshop on Recruiting Participants for Empirical Software Engineering (RoPES’22)_ , Vol. 5. Association for Computing Machinery.
  * Russo (2022) Daniel Russo. 2022.  Recruiting Software Engineers on Prolific. In _the1st International Workshop on Recruiting Participants for Empirical Software Engineering (RoPES 2022)_ , Vol. 5. Association for Computing Machinery.
  * Shapiro and Wilk (1965) Samuel Shapiro and Martin Wilk. 1965.  An analysis of variance test for normality (complete samples).  _Biometrika_ 52, 3/4 (1965), 591–611.
  * Sharif and Maletic (2010) Bonita Sharif and Jonathan I. Maletic. 2010.  An eye tracking study on camelcase and under-score identifier styles. In _IEEE 18th International Conference on Program Comprehension_. 196–205.  https://doi.org/10.1109/ICPC.2010.41
  * Siegmund et al. (2014) Janet Siegmund, Christian Kästner, Jörg Liebig, Sven Apel, and Stefan Hanenberg. 2014.  Measuring and modeling programming experience.  _Empirical Software Engineering_ 19, 5 (2014), 1299–1334.  https://doi.org/10.1007/s10664-013-9286-4
  * Sparman and Schulte (2023) Sören Sparman and Carsten Schulte. 2023.  Analysing the API learning process through the use of eye tracking. In _the 2023 Symposium on Eye Tracking Research and Applications (ETRA)_. 1–6.
  * Von Mayrhauser and Vans (1998) Anneliese Von Mayrhauser and A. Marie Vans. 1998.  Program Understanding Behavior During Adaptation of Large Scale Software. In _the 6th International Workshop on Program Comprehension (IWPC’98)_. 164–172.