1. Introduction

Fractional models, in which a standard time or space differential operator are replaced by a corresponding fractional operator, have gained considerable popularity and importance during the last few decades, although fractional calculus is an old topic in mathematics, see ^[1] for historical notes. Fractional calculus is now used to describe a broad range of non-classical phenomena in the applied sciences, engineering, and finance due to the intrinsic non-local property of fractional derivatives, for example, the filtration of solutes in porous soils ^[2], diffusion of water molecules in brain tissues ^[3], electrical charge transport in polymer networks ^[4], the relationship between certain option pricing and heavy-tailed stochastic process ^[5], anomalous diffusion process for continuous time random walk models ^[6].

It is well known that the derivation of the analytical solutions to the fractional differential equations is generally difficult and computation of them is very expensive due to infinite series in the analytical solutions. On the other hand, the implementation of numerical approaches to solve the fractional differential equations also has essential difficulties and challenges due to the non-local nature of the fractional operators (space fractional) and the dependence on the full history (time fractional). However, in recent years, a number of successful numerical approaches for fractional differential equations have been considered such as finite difference methods ^{[7,8,9,10,11]}, spectral methods ^[12,13], finite element methods ^{[14,15,16,17]}, and discontinuous Galerkin methods ^[18,19]. Many of these approaches have limitations in terms of computational efficiency when two and three spatial dimensions are considered. Recently, Yang et al. in ^[20] proposed a new approach using a matrix transfer technique with finite difference and finite element methods to solve the time-space fractional diffusion equation in two spatial dimensions with homogeneous Dirichlet boundary conditions. The solution is advanced in time by computing the function of a matrix times a vector by the preconditioned Lanczos method. This concept was also considered in ^[21] using the finite element method in space and a semi-implicit Euler approximation in time. The computation of the fractional power of a matrix times a vector was done by the contour integral method, the extended Krylov subspace method, and the preassigned poles and interpolation nodes method.

We here concern ourselves with the following fractional-in-space phase-field models:

● Allen-Cahn type

$u_t + \big(-\Delta\big)^{\alpha} u + \frac{1}{\epsilon} F^{\;'}(u) = 0 \qquad \quad \;\;\;\; \hbox{in } \Omega \times (0, T),$

(1.1a)

$u = g^D \qquad \quad \;\; \hbox{in } \Gamma^D \times (0, T),$

(1.1b)

$\frac{\partial u}{\partial n} = g^N \qquad \quad \; \; \hbox{in } \Gamma^N \times (0, T),$

(1.1c)

$u(x, 0) = u_0(x) \qquad \; \; \hbox{in } \Omega.$

(1.1d)

● Cahn-Hilliard type

$u_t = \Delta w \qquad \qquad \qquad \qquad \quad \; \hbox{ in } \Omega \times (0, T),$

(1.2a)

$w = \big(-\Delta\big)^{\alpha} u + \frac{1}{\epsilon} F^{\;'}(u) \qquad \quad \hbox{in } \Omega \times (0, T),$

(1.2b)

$u = g^D \qquad \qquad \qquad \qquad \quad \;\; \hbox{ in } \Gamma^D \times (0, T),$

(1.2c)

$\frac{\partial u}{\partial n} = \frac{\partial w}{\partial n} = 0 \qquad \qquad \qquad \quad \;\; \hbox{ in } \Gamma^N \times (0, T),$

(1.2d)

$u(x, 0) = u_0(x) \qquad \qquad \qquad \qquad \;\; \hbox{ in } \Omega.$

(1.2e)

The function $F$ in (1.1) and (1.2) represents a configuration potential which may have two (or more) wells. The general form of $F$ is given by

$F(v) = r(v) - \frac{\lambda}{2} v^2,$

(1.3)

where $r(\cdot)$ is a function and $\lambda \geq 0$ is a constant. Upon the smoothness of the function $r(v)$, several types of significant choices have been proposed for $F$. For instance, the non-smooth logarithmic potential is ^{[22,23,24,25]}

$F(v) = \frac{\theta}{2} \big[(1+v)\ln(1+v) + (1-v)\ln(1-v) \big] - \frac{\theta_c}{2} v^2$

(1.4)

with $0 < \theta \leq \theta_c$, where $\theta$ and $\theta_c$ are the absolute temperature and the transition temperature, respectively. The smooth and convex one is the standard double-well potential ^[26,27], namely,

$F(v) = \frac{1}{4} (v^2 -1)^2.$

(1.5)

It is an approximation of the logarithmic potential in case the absolute temperature $\theta$ is close to the transition temperature $\theta_c$.

Then, $f(v): = F^{\; '}(v) = v(v^2-1)$ represents the bistable non-linearity for the double-well potential (1.5), whereas $f(v) = \frac{\theta}{2} \ln \left(\frac{1+v}{1-v} \right) -\theta_c v$ is for the logarithmic free energy (1.4). Further, the fractional-in-space phase-field models (1.1)-(1.2) can be viewed as the gradient flow of the energy

$\mathcal{E}(v) = \int_{\Omega} \left( \frac{1}{2} |\nabla v|^{(2\alpha)} + \frac{1}{\epsilon} F(v) \right) \; d \Omega.$

(1.6)

In problem (1.1)-(1.2), $u$ represents the concentration of one of the species of the alloy, $w$ is the chemical potential, the parameter $\epsilon$ represents the diffuse interface width parameter, $\Omega \subset \mathbb{R}^d (d = 1, 2, 3)$ is a bounded domain. The operator $\big(-\Delta\big)^{\alpha}$ denotes the fractional operators of order $\alpha \in (0.5, 1]$, see, e.g., ^[28,29,30].

The system (1.1) (or the system (1.2)) with $\alpha = 1$, is known as a scaled in time form of the Allen-Cahn (or Cahn-Hilliard) equation. The Cahn-Hilliard equation was originally introduced by Cahn and Hilliard in ^[31] to describe the phase separation and coarsening phenomena in a melted alloy, whereas the Allen-Cahn equation was introduced by Allen and Cahn in ^[32] to describe the motion of anti-phase boundaries in crystalline solids. The Allen-Cahn/Cahn-Hilliard equations are essential building blocks in the phase field methodology or the diffuse interface methodology for moving interface and free boundary problems, see, e.g., ^[33,34]. Both equations are particular cases of gradient flows, written as $u_t = -\frac{\delta E(u)}{\delta u}$, where $\frac{\delta E(u)}{\delta u}$ stands for the variational derivative of the free energy, either taken in the $L^2$-norm for Allen-Cahn or in $H^{-1}$-norm for Cahn-Hilliard. Since the Allen-Cahn/Cahn-Hilliard equations are gradient flows of the energy functional, the total energy is always non-increasing for both model equations. On the other hand, the Cahn-Hilliard model has a mass conservative property in contrast to the Allen-Cahn model.

There are several challenges for obtaining numerical approximations of these problems such as the existence of a nonlinear term and the presence of the small interfacial length parameter $\epsilon$. An appropriate numerical scheme requires a proper relation between physical and numerical scales, that is, the size of spatial mesh $h$ and time step $\tau$ have to properly be related to the interaction length $\epsilon$. In the literature, the classic Cahn-Hilliard equation has been studied with well known methods like spectral methods ^[35,36], collocation methods ^[37], finite element methods ^[22,38,39], discontinuous Galerkin methods ^{[39,40,41,42]}, whereas the classic Allen-Cahn equation has been investigated in ^[43] for finite element methods, and in ^[44,45] for discontinuous Galerkin methods. The resulting system obtained after the spatial discretization is an inherently stiff system due to the small positive parameter $\epsilon$. This is then handled by appropriate temporal discrimination methods, such as the implicit approximations (which are not energy-stable) ^[26,27], the implicit-explicit (IMEX) techniques (which do not introduce any numerical dissipation) ^{[46,47,48,49]}, and the average vector field (AVF) method ^[45,50]. In addition, in order to obtain unconditionally energy-stable schemes, splitting of the potential $F(v)$ into a convex and a non-convex part ^{[51,52,53,54,55]}, Taylor's formula ^[56] and a Lagrange multiplier formulation of the potential term ^[57,58] are considered.

Recently, there has been a fast increasing number of studies on front propagation of reaction diffusion systems with an anomalous diffusion as super diffusion, i.e., the fractional Allen-Cahn/Cahn-Hilliard equations. A reason for considering the fractional system (1.2) can be proved by observing that the Laplace operator (1.2b) in the Cahn-Hilliard equation ^[31] was actually replaced by a spatial convolution term in order to describe long-range interactions among particles. Due to analytical reasons, this nonlocal term has been substituted with the term $-\Delta u$. However, the use of the fractional Laplacian $\big(-\Delta\big)^{\alpha} u$ appears to be more adherent to the physical setting. The analysis of nonlocal Allen-Cahn and Cahn-Hilliard models has been studied in ^{[59,60,61,62,63]}. Especially, the fractional models offer insight that traditional approaches do not offer, such as in the case of diffusion in heterogeneous environments. Such super diffusion is related to Lïvy processes and can be modeled by a fractional operator $(-\Delta)^{\alpha}$ with $0.5 < \alpha \leq 1$. Ilić et al. in ^[66] showed that the fractional Laplace operator $(-\Delta)^{\alpha}$ has the same interpretation as $(-\Delta)$ in terms of its spectral decomposition for homogeneous boundary conditions. Further, the matrix transfer technique was introduced in ^[28,29] to compute the fractional Laplacian by first computing a matrix representation of the Laplace (independent of the discretization approach) and then raising it to the fractional order. It is noted that the fractional Laplacian $(-\Delta)^{\alpha}$ is a nonlocal operator and any approximation of it results in a large dense matrix. However, the sparse approximation of $(-\Delta)$ can be captured directly in numerical implementations, when the matrix transfer technique is applied.

We here solve the equation of the form (1.1) or (1.2) by computing the fractional power of a matrix times a vector. To compute the fractional power, the contour integral method proposed in ^[30] is applied. It is expected that the fractional reaction-diffusion models as (1.1) or (1.2) with smaller fractional order exhibit more heterogeneous environments. In addition, the sharp gradients and singularities emerge locally for small values of the parameter $\epsilon$. To handle these difficulties, we apply the symmetric interior penalty Galerkin (SIPG) method as a discontinuous Galerkin method for the spatial discretization. In contrast to continuous finite elements, the space of solutions or test functions in discontinuous Galerkin methods consist of piecewise discontinuous polynomials. That is, no continuity constraints are explicitly imposed on the state and test functions across the element interfaces. In this way, the SIPG approximation allows to capture singularities locally. They also allow the use of highly nonuniform and unstructured meshes, and have built-in parallelism which permits coarse-grain parallelization. In addition, the fact that the mass matrices are block diagonal is an attractive feature in the context of time dependent problems. Further, the fact that methods of this kind are locally conservative, which is a particularly relevant feature in the realm of numerical approximation of nonlinear hyperbolic conservation laws. On the other hand, the implicit-explicit (IMEX) methods are applied for the temporal discretization. In order to save computational cost we have addressed an adaptive-time stepping algorithm based on the difference between the first order IMEX method and the second order IMEX method.

The remainder of this paper is organized as follows: in the next section, we introduce the symmetric interior penalty Galerkin (SIPG) method as a discontinuous Galerkin discretization. In Section 3, we review the contour integral method, which allows us to approximate the fractional Laplacian by a fractional power of a matrix. The implicit-explicit methods are given in Section 4 for the temporal discretization. Also, an adaptive-time stepping algorithm is addressed to reduce the computational cost. In Section 6, several numerical examples are presented to show the effect of the fractional power. The paper ends with some conclusions and remarks.

---

2. Symmetric interior penalty Galerkin (SIPG) discretization

In this section, we introduce the symmetric interior penalty Galerkin (SIPG) discretization as a discontinuous Galerkin (DG) method. It is chosen due to the symmetric property of its bilinear form, i.e., $a_h(y, v) = a_h(v, y)$, see, e.g., ^[67].

We begin with the continuous weak formulation of the classical Allen-Cahn equation by choosing $\alpha = 1$

$u_t -\Delta u + \frac{1}{\epsilon}f(u) = 0 \qquad \quad \;\;\; \hbox{in } \Omega \times (0, T).$

(2.1)

Then, find $u(t) \in U$ such that

$(u_t, v) + a(u, v) = l(v) \quad \quad \quad \quad \forall v \in V, \;\; t \in (0, T],$

(2.2a)

$(u(\cdot, 0), v) = (u_0, v) \quad \quad \quad \quad \forall v \in V,$

(2.2b)

where the space of solutions $U$, and the space of test functions are defined by

$U = \{u \in H^{1}(\Omega): \; \; u|_{\Gamma^D} = g^D\}, \qquad V = \{v \in H^{1}(\Omega): \; \; v|_{\Gamma^D} = 0\},$

and the (bi)-linear forms are given by

$a(u, v) = \int_{\Omega} (\nabla u \cdot \nabla v ) \; dx, \quad \hbox{and} \quad l(v) = -\int_{\Omega} \frac{1}{\epsilon} f(u) v \; dx.$

We assume that the domain $\Omega$ is polygonal such that the boundary is exactly represented by boundaries of triangles. We denote $\{ \mathcal{T}_h\}_h$ as a family of shape-regular simplicial triangulations of $\Omega$. Each mesh $\mathcal{T}_h$ consists of closed triangles such that $\overline{\Omega} = \bigcup_{K \in \mathcal{T}_h} \overline{K}$ holds. We assume that the mesh is regular in the following sense: for different triangles $K_i, K_j \in \mathcal{T}_h$, $i \not = j$, the intersection $K_i \cap K_j$ is either empty or a vertex or an edge, i.e., hanging nodes are not allowed. The diameter of an element $K$ and the length of an edge $E$ are denoted by $h_{K}$ and $h_E$, respectively.

We split the set of all edges $\mathcal{E}_h$ into the set $\mathcal{E}^0_h$ of interior edges, the set $\mathcal{E}^{D}_h$ of Dirichlet boundary edges, and the set $\mathcal{E}^{N}_h$ of Neumann boundary edges so that $\mathcal{E}_h = \mathcal{E}^{B}_h \cup \mathcal{E}^{0}_h$ with $\mathcal{E}^B_h = \mathcal{E}^{D}_h \cup \mathcal{E}^{N}_h$. Let the edge $E$ be a common edge for two elements $K$ and $K^e$. For a piecewise continuous scalar function $u$, there are two traces of $u$ along $E$, denoted by $u|_E$ from inside $K$ and $u^e|_E$ from inside $K^e$. The jump and average of $u$ across the edge $E$ are defined by:

$\left[\kern-0.15em\left[u \right]\kern-0.15em\right] = u|_E\mathbf{n}_{K}+u^e|_E\mathbf{n}_{K^e}, \quad \left\{ {\left\{ u \right\}} \right\} = \frac{1}{2}\big( u|_E+u^e|_E \big),$

(2.3)

where $\mathbf{n}_K$ (resp. $\mathbf{n}_{K^e}$) denotes the unit outward normal to $\partial K$ (resp. $\partial K^e$).

Similarly, for a piecewise continuous vector field $\nabla u$, the jump and average across an edge $E$ are given by

$\left[\kern-0.15em\left[{\nabla u} \right]\kern-0.15em\right] = \nabla u|_E \cdot \mathbf{n}_{K}+\nabla u^e|_E \cdot \mathbf{n}_{K^e}, \quad \left\{ {\left\{ {\nabla u} \right\}} \right\} = \frac{1}{2}\big(\nabla u|_E+\nabla u^e|_E \big).$

(2.4)

For a boundary edge $E \in K \cap \Gamma$, we set $\left\{ {\left\{ {\nabla u} \right\}} \right\} = \nabla u$ and $\left[\kern-0.15em\left[u \right]\kern-0.15em\right] = u\mathbf{n}$, where $\mathbf{n}$ is the outward normal unit vector on $\Gamma$.

For continuous finite element methods (FEMs), the idea is to approximate (2.1) using a conforming, finite dimensional space $V_h \subset V$. On the other hand, we point out that in discontinuous Galerkin methods the space of solutions or test functions consist of piecewise discontinuous polynomials. That is, no continuity constraints are explicitly imposed on the state and test functions across the element interfaces. As a consequence, weak formulations must include jump terms across interfaces, and typically penalty terms are added to control the jump terms. Then, we define the spaces of test functions, of the solutions by

$V_h = U_h = \{{u \in L^2(\Omega)}:{ u\mid_{K}\in \mathbb{P}^r(K) \quad \forall K \in \mathcal{T}_h}\},$

(2.5)

where $\mathbb{P}^r(K)$ is the set of polynomials of degree at most $r$ in $K$. Note that the space $U_h$ of discrete solutions and the space of test functions $V_h$ are identical due to the weak treatment of boundary conditions in DG methods. Note also that the space $V_h$ is non-conforming such that $V_h \not\subset V.$

Now, we are ready to set up the SIPG discretization of the continuous weak formulation (2.1). Multiply (2.1) by a test function $v \in V_h$, and then integrate over each element $K \in \mathcal{T}_h$

$\sum \limits_{K \in \mathcal{T}_h} \int_{K} \Big( u_t v - \Delta u \, v \Big) \; dx = -\sum \limits_{K \in \mathcal{T}_h } \int_{K} \frac{1}{\epsilon} f(u) v \; dx.$

An application of integration by parts on each element integral and the definition of the jump operator give us

$\begin{align*} \sum \limits_{K \in \mathcal{T}_h} \int_{K} \Big( u_t v + \nabla u \cdot \nabla v \Big) \; dx &- \sum \limits_{E \in \mathcal{E}^0_h \cup \mathcal{E}^D_h} \int_{E} \left[\kern-0.15em\left[{v\nabla u} \right]\kern-0.15em\right] \; ds \\ = & -\sum \limits_{K \in \mathcal{T}_h } \int_{K} \frac{1}{\epsilon} f(u) v \; dx + \sum \limits_{E \in \mathcal{E}^N_h} \int_{E} g^N v \; ds. \end{align*}$

The following equality

$\left[\kern-0.15em\left[{v\nabla u} \right]\kern-0.15em\right] = \left\{ {\left\{ {\nabla u} \right\}} \right\} \cdot \left[\kern-0.15em\left[v \right]\kern-0.15em\right] + \left[\kern-0.15em\left[{\nabla u} \right]\kern-0.15em\right] \cdot \left\{ {\left\{ v \right\}} \right\},$

which one can verify easily and the fact that $\left[\kern-0.15em\left[{\nabla u} \right]\kern-0.15em\right] = 0$ ($u$ is assumed to be smooth) yield

$\begin{align*} \sum \limits_{K \in \mathcal{T}_h} \int_{K} \Big( u_t v + \nabla u \cdot \nabla v \Big) \; dx &- \sum \limits_{E \in \mathcal{E}^0_h \cup \mathcal{E}^D_h} \int_{E} \left\{ {\left\{ {\nabla u} \right\}} \right\} \cdot \left[\kern-0.15em\left[v \right]\kern-0.15em\right] \; ds \\ = & -\sum \limits_{K \in \mathcal{T}_h } \int_{K} \frac{1}{\epsilon} f(u) v \; dx + \sum \limits_{E \in \mathcal{E}^N_h} \int_{E} g^N v \; ds. \end{align*}$

To handle the coercivity of the left hand side and control the jump terms, we add the following equalities via $\left[\kern-0.15em\left[u \right]\kern-0.15em\right] = 0$ on the interior edges $E \in \mathcal{E}^0_h$

$\begin{eqnarray*} \sum \limits_{E \in \mathcal{E}^0_h \cup \mathcal{E}^D_h } \int_{E} \left\{ {\left\{ {\nabla v} \right\}} \right\} \cdot \left[\kern-0.15em\left[u \right]\kern-0.15em\right] \; ds & = & \sum \limits_{E \in \mathcal{E}^{D}_h} \int_{E} g^D \left\{ {\left\{ {\nabla v} \right\}} \right\} \; ds, \\ \sum \limits_{\mathcal{E}^0_h \cup \mathcal{E}^D_h } \frac{\sigma }{h_E} \int_{E} \left[\kern-0.15em\left[u \right]\kern-0.15em\right] \cdot \left[\kern-0.15em\left[v \right]\kern-0.15em\right] \; ds & = & \sum \limits_{E \in \mathcal{E}^{D}_h} \frac{\sigma }{h_E} \int_{E} g^D \, \left[\kern-0.15em\left[v \right]\kern-0.15em\right] \; ds, \end{eqnarray*}$

where $\sigma$ is the penalty parameter, which should be chosen sufficiently large to ensure the stability of the SIPG scheme, see, e.g., ^[67].

Then, the weak formulation of the Allen-Cahn equation (2.1), discretized by the SIPG method reads as: find $u_h \in U_h$ such that

$(\frac{\partial u_h}{\partial t}, v) + a_h(u_h, v) = l_h(v) \quad \quad \quad \quad \forall v \in V_h, \;\; t \in (0, T],$

(2.6a)

$(u_h(\cdot, 0), v) = (u_0, v) \quad \quad \quad \quad \forall v \in V_h,$

(2.6b)

where the (bi)-linear forms are given by

$\begin{align} a_h(u, v) = & \sum \limits_{K \in \mathcal{T}_h} \int \limits_{K} \big(\nabla u \cdot \nabla v )\; dx - \sum \limits_{ E \in \mathcal{E}^0_h \cup \mathcal{E}^D_h} \int \limits_E \Big( \left\{ {\left\{ {\nabla u} \right\}} \right\} \cdot \left[\kern-0.15em\left[v \right]\kern-0.15em\right] + \left\{ {\left\{ {\nabla v} \right\}} \right\} \cdot \left[\kern-0.15em\left[u \right]\kern-0.15em\right] \Big) \; ds \nonumber \\ & + \sum \limits_{ E \in \mathcal{E}^0_h \cup \mathcal{E}^D_h} \frac{\sigma }{h_E} \int \limits_E \left[\kern-0.15em\left[u \right]\kern-0.15em\right] \cdot \left[\kern-0.15em\left[v \right]\kern-0.15em\right] \; ds, \end{align}$

(2.7a)

$l_h(v) = -\sum \limits_{K \in \mathcal{T}_h} \int \limits_{K} \frac{1}{\epsilon} f(u)v \; dx + \sum \limits_{E \in \mathcal{E}_h^{D}} \int \limits_E g^D \big(\frac{\sigma }{h_E} \left[\kern-0.15em\left[v \right]\kern-0.15em\right] -\left\{ {\left\{ {\nabla v} \right\}} \right\} \big) \; ds + \sum \limits_{E \in \mathcal{E}^N_h} \int_{E} g^N v \; ds,$

(2.7b)

and $u_h(\cdot, 0)$ is a $L_2$-orthogonal projection of the initial condition $u_0$ onto $U_h$.

Analogously, the weak formulation of the classical Cahn-Hilliard equation ($\alpha = 1$ in (1.2)), discretized by the SIPG method reads as: find $u_h(t), w_h(t) \in U_h$ such that

$(\frac{\partial u_h}{\partial t}, v) + a_h(w_h, v) = 0 \quad \quad \quad \quad \forall v \in V_h, \;\; t \in (0, T],$

(2.8a)

$-(w_h, v_h) + a_h(u_h, v) = l_h(v) \quad \quad \quad \quad \forall v \in V_h, \;\; t \in (0, T],$

(2.8b)

$(u_h(\cdot, 0), v) = (u_0, v) \quad \quad \quad \quad \forall v \in V_h.$

(2.8c)

For each time step, we can expand the discrete solutions as

$u_h(t) = \sum \limits_{i = 1}^{N} \sum \limits_{j = 1}^{n_{p}} U_{j}^{i} \phi_{j}^{i}, \qquad w_h(t) = \sum \limits_{i = 1}^{N} \sum \limits_{j = 1}^{n_{p}} W_{j}^{i} \phi_{j}^{i},$

(2.9)

where $U_{j}^{i}, W_{j}^{i}$, and $\phi_{j}^{i}$ are the unknown coefficients and the basis functions, respectively, for $j = 1, 2, \cdots, n_{p}$ and $i = 1, 2, \cdots, N$. The number $N$ denotes the number dG elements and $n_{p}$ is the local dimension of each dG element with

$n_{p} = \frac{(p+1)(p+2)}{2},$

where $p$ is the degree of the polynomial order.

Inserting (2.9) into (2.6) and (2.8), we obtain the semi-discrete formulations of the Allen-Cahn equation

$\frac{dU}{dt} + \mathbf{M}^{-1} \mathbf{L} U = \mathbf{M}^{-1}\mathbf{B}(U),$

(2.10)

and of the Cahn-Hilliard equation

$\frac{dU}{dt} + \mathbf{M}^{-1} \mathbf{L} W = 0,$

(2.11a)

$- W + \mathbf{M}^{-1} \mathbf{L} U = \mathbf{M}^{-1}\mathbf{B}(U),$

(2.11b)

where $U, W$ are the unknown coefficient vectors $U = (U_{1}^{1}, \cdots, U^{1}_{n_p}, \cdots, U_{N}^{1}, \cdots, U^{N}_{n_p})$, $W = (W_{1}^{1}, \cdots, W^{1}_{n_p}, \cdots, W_{N}^{1}, \cdots, W^{N}_{n_p})$, $\mathbf{M}$ is the mass matrix, $\mathbf{L}$ is the stiffness matrix corresponding to $a_h(u, v)$, and $\mathbf{B}(\cdot)$ is the nonlinear vector of the unknown coefficient vector $U$ corresponding to $l_h(v)$.

We are now ready to employ the matrix transfer technique introduced in ^[28], which states that the error introduced by approximating the fractional Laplacian by a fractional power of the matrix $\mathbf{A} = \mathbf{M}^{-1}\mathbf{L}$ converges at the same rate as the underlying discretization method. Then, the fractional Laplacian operator is approximated as [20,p. 1162,eqn. (2.2)]

$(-\Delta)^{\alpha} u \approx \mathbf{A}^{\alpha} U.$

(2.12)

In the following section, we employ the contour integral method introduced in ^[30] to compute the fractional power of $\mathbf{A}$ times a vector.

---

3. Contour integration method (CIM)

In this section, we review the contour integral method, which allows us to approximate the fractional Laplacian by a fractional power of a matrix. An analytic function $h$ of a square matrix $\mathbf{A}$ can be represented as a contour integral in the complex plane [68,Definition 1.11]

$h(\mathbf{A}) = \frac{1}{2 \pi i} \int_{\Gamma} h(z) (z \mathbf{I}-\mathbf{A})^{-1} dz,$

(3.1)

where $i\equiv \sqrt{-1}$, and $\Gamma$ is a closed contour lying in the region of analyticity of $h$ and enclosing the spectrum of $\mathbf{A}$. Then, a numerical quadrature method is applied to the integration (3.1) to approximate $h(\mathbf{A})$.

We here compute the vector $h(\mathbf{A})\mathbf{b}$ for a given vector $\mathbf{b}$ by using the definition (3.1) with the technique proposed in ^[30]. The basic principle is based on an application of the midpoint rule over a circle contained within an annulus whose outer boundary maps to the interval $(-\infty, 0]$ and whose inner boundary maps to the interval $[\lambda_1, \lambda_N]$, which are the eigenvalues of $\mathbf{A}$, see Figure 1 (taken from ^[21]).

Then, the vector $h(\mathbf{A})\mathbf{b}$ is computed via the following quadrature formula

$h(\mathbf{A})\mathbf{b} \approx \hbox{Im} \sum \limits_{i = 1}^{n_q} w_j (\eta_i \mathbf{I}-\mathbf{A})^{-1} \mathbf{b},$

(3.2)

where the weights and shifts are denoted by $w$ and $\eta$, respectively, and $n_q$ is the number of quadrature points.

The SIPG discretization of the Dirichlet problem provides a nonsingular and real-symmetric matrices $\mathbf{L}$ and $\mathbf{M}$. Then, by using the symmetry property of these matrices, we can integrate over only the upper half the contour. The algorithm based on the method in ^[20] is given in Algorithm 1. In the algorithm, we use the routines ellipkkp and ellipjc, which are described in ^[69] to compute complex arguments.

Algorithm 1 CIM for computing $\mathbf{A}^{\alpha}$ for the Dirichlet problem
1: l = eigs(L, M, 1, 'SM'); l1 = l(1);	% min. eigenvalue of A
2: l = eigs(L, M, 1, 'LM'); lN = l(1);	% max. eigenvalue of A
3: k = (sqrt(lN/l1)-1)/(sqrt(lN/l1)+1);	% a convenient constant
4: [K Kp] = ellipkkp(-log(k)/pi);	% elliptic integrals
5: t =.5i*Kp-K+(nq-.5:-1:0)*2*K/nq;	% midpoint rule points
6: [sn cn dn] = ellipjc(t, -log(k)/pi);	% jacobi elliptic functions
7: xi = sqrt(l1*lN)*(1/k+sn)./(1/k-sn);	% quadrature nodes
8: dxidt = cn.*dn./(1/k-sn). ^ 2;	% derivative wrt t
9: wts = h(xi).*dxidt;	% quadrature weights
10: v = zeros(length(b), 1);	% initialize output
11: for $i=1:n_q$ do
12:     y = (xi(j)*M-L)$\backslash$ (M*b);
13:     v = v + wtj(j)*y;	% update solution vector
14: end for
15: v = -4*K*sqrt(l1*lN)*imag(v)/(k*pi*nq);	% scale the solution

However, many applications require Neumann-type boundary conditions. As done in the Dirichlet problem, a contour $\Gamma$ surrounding the eigenvalues of $\mathbf{A}$ cannot be found. Therefore, the Algorithm 1 should be modified to compute (3.2). Burrage et al. in [21,Sec. 4] handle this problem by adding a correction term. It is shown in Algorithm 2. The first line of the Algorithm 2 yields the first non-zero eigenvalue of the matrix $\mathbf{A}$.

Algorithm 2 CIM for computing $\mathbf{A}^{\alpha}$ for the Neumann problem
1: l = eigs(L, M, 3, 'SM'); l1 = l(1);	% min. eigenvalue of A
2: l = eigs(L, M, 1, 'LM'); lN = l(1);	% max. eigenvalue of A
3:                 $\vdots$
4: v = -4*K*sqrt(l1*lN)*imag(v)/(k*pi*nq);	% scale the solution
5: e = ones(length(b), 1);
6: v = v + (e'*(M*(b-v)))/(e'*M*e)*e;	% corrector term

---

4. Implicit-explicit schemes

After spatial discretization of the phase-field models, the leading system is typically stiff for small values of the parameter $\epsilon$. Explicit methods are not suitable for stiff systems, whereas implicit methods require the solution of nonlinear equations at each time step. Therefore, the implicit-explicit (IMEX) method can play an important role for such problems, see, e.g., ^[70,71]. In such a procedure, the Laplacian term is discretized implicitly in time and the nonlinear terms are discretized explicitly. This can also be recognized and analyzed as a splitting technique. In addition, it typically allows for a larger time step than explicit methods, while avoiding the use of nonlinear solvers.

We first divide the time interval $[0, T]$ as follows

$0 = t_0 < t_1 < \cdots < t_{N_T} = T$

with the time step size $\tau_n = t_{n}-t_{n-1}, \; \; n = 1, 2, \cdots, N_T$. Then, we can consider the first-and second-order IMEX approximations of the following system of ordinary differential equations (ODEs) for the Allen-Cahn model

$\frac{dU}{dt} + \mathbf{A}^{\alpha} U = \mathbf{M}^{-1} \mathbf{B}(U)$

(4.1)

and for the Cahn-Hilliard model

$\frac{dU}{dt} + \mathbf{A} W = 0,$

(4.2a)

$- W + \mathbf{A}^{\alpha} U = \mathbf{M}^{-1} \mathbf{B}(U).$

(4.2b)

---

4.1. First-order implicit-explicit scheme

We here consider the semi-implicit scheme as a first-order IMEX scheme:

● The Allen-Cahn model (2.10)

$\frac{ U^{n+1} -U^n}{\tau_n} + \mathbf{A}^{\alpha} U^{n+1} = \mathbf{M}^{-1} \mathbf{B}(U^n).$

(4.3)

● The Cahn-Hilliard model (2.11)

$\frac{ U^{n+1} -U^n}{\tau_n} + \mathbf{A} W^{n+1} = 0,$

(4.4a)

$\mathbf{A}^{\alpha} U^{n+1} - W^{n+1} = \mathbf{M}^{-1} \mathbf{B}(U^n).$

(4.4b)

---

4.2. Second-order implicit-explicit schemes

The modified Crank-Nicolson/Adams-Bashforth scheme is applied for the Allen-Cahn model:

$\frac{U^{n+1}-U^n}{\tau_n} + \mathbf{A}^{\alpha} \Big( \frac{9}{16} U^{n+1} + \frac{3}{8} U^n + \frac{1}{8} U^{n-1} \Big) = \frac{3}{2} \mathbf{M}^{-1}\mathbf{B}(U^n)- \frac{1}{2} \mathbf{M}^{-1}\mathbf{B}(U^{n-1}),$

(4.5)

while the Crank-Nicolson/Adams-Bashforth scheme is applied for the Cahn-Hilliard model:

$\frac{ U^{n+1} -U^n}{\tau_n} + \mathbf{A} W^{n+1} = 0,$

(4.6a)

$\frac{1}{2}\mathbf{A}^{\alpha} \big(U^{n+1} + U^{n} \big) - W^{n+1} = \frac{3}{2} \mathbf{M}^{-1}\mathbf{B}(U^n)- \frac{1}{2} \mathbf{M}^{-1}\mathbf{B}(U^{n-1}).$

(4.6b)

Remark 4.1. In our numerical simulations, we use two different matrix functions $h(z)$ to compute the fractional matrix $\mathbf{A}^{\alpha}$, which is formed on both left-and right-hand side of the IMEX schemes. When we first apply the Laplace transform and then Laplace inversion to the ODE system (2.10) or (2.11), the matrix function $h(z)$ is defined in terms of an exponential function, see ^[72]. Then, we have $h(z)$ as

$h(z) = \frac{1}{\exp(\tau z^{\alpha})}$

for the left-hand side. Note that the fraction is due to inversion in the IMEX schemes. On the other hand, for the right-hand side, we choose as $h(z) = z^{\alpha}$ as was done in ^[30].

Remark 4.2. The phase-field models such as Allen-Cahn equations and Cahn-Hilliard are obtained as gradient flows of an energy functional. Therefore, the total energy is always non-increasing for both model equations in the continuous setting. In our numerical simulations, we show the non-increasing property of the energy. However, it can not be justified at a theoretical level due to the explicit treatment of the nonlinear terms.

---

4.3. Time-step size adaptivity

For small values of the parameter $\epsilon$, the transition layer moves slowly and then an inordinate number of time steps is required to resolve the dynamics response of the fractional-in-space phase-field system. To reduce the amount of work, adaptivity in time should be used. Our time-step size adaptivity is based on the ideas presented in ^[73,74]. To update the time-step size, we use the difference between two solutions, which are a predictor and a corrector. The first-order IMEX schemes are chosen as a predictor, whereas the second-order IMEX schemes are chosen as a corrector. The time-step adaptive algorithm is presented in Algorithm 3. To update the time-step size, we use the following controller

$\tau^{*}_{n+1} = \rho \left( \frac{Tol}{e_n} \right)^{1/2} \tau_n,$

(4.7)

where $\rho$ is a safety coefficient, which is introduced to reduce the probability of rejecting $\tau^{*}_{n+1}$. In numerical examples, we take $\rho = 0.9$ as suggested in ^[75]. The parameter $Tol$ determines the required accuracy of the numerical solution. The impact of $Tol$ on the number of times steps will be studied in Section 6. Finally, to avoid a strong increase or decrease of subsequent time steps, we use the following formula as proposed in the deterministic framework ^[76]

Algorithm 3 Time-step adaptive algorithm

1: Given $U_0$, $\tau_0$, Tol

2: for $n=1, 2, \cdots N_T$ do

3:     Compute $P_n$ using a first-order implicit-explicit scheme.

4:     Compute $C_n$ using a second-order implicit-explicit scheme.

5:     Calculate $e_n = \frac{\|C_n -P_n\|}{\|C_n\|}$.

6:     Set $reject=0$.

7:     if $e_n > Tol$ then

8:         Recalculate time-step size $\tau_n \leftarrow \Lambda(e_n, \tau_n)$.

9:         Update $reject=reject +1$.

10:        goto step 3.

11:     else

12:         Update time-step size $\tau_{n+1}=\Lambda(e_n, \tau_n)$.

13:         continue

14:     end if

15: end for

$\Lambda(e_n, \tau_n) = \min \{s_{max} \tau_n, \max \{s_{min} \tau_n, \tau^{*}_{n+1} \} \}.$

(4.8)

In numerical simulations, $s_{min} = 0.1$ and $s_{max} = 2$ are used. A step size is accepted if $e_n \leq Tol$, otherwise it is rejected.

The time-step size adaptivity allows us to reduce the computation time by factors of hundreds compared to the uniform step size.

---

5. The discrete energy stability

In this section, we show the discrete energy stability property of the fractional Allen-Cahn equation (1.1) with the double-well potential (1.5) as done in ^[81]. After being semi-discretized in space, the following ODE system is obtained

$\frac{dU}{dt} + \mathbf{A}^{\alpha} U + \frac{1}{\epsilon} \left( U^3-U \right) = 0.$

(81)

If we define the semi-discrete energy as:

$\mathcal{E}_h(U) = \frac{1}{4\epsilon} \sum \limits_{i} (U_i^2-1)^2 + \frac{1}{2}U^T \mathbf{A}^{\alpha} U,$

(52)

then the ODE system (5.1) can be viewed as the gradient flow of the energy $\mathcal{E}_h(U)$, i.e.,

$\frac{dU}{dt} = - \nabla_U \mathcal{E}_h(U).$

(5.3)

Taking the $L^2$ inner product of (5.3) with $-\frac{dU}{dt}$ gives us

$\frac{d \mathcal{E}_h(U)}{dt} = - \| \frac{dU}{dt} \|^2_2 \leq 0,$

(5.4)

which implies that the energy is non-increasing. Now, we show that this energy stability can be inherited by the first-order scheme (4.3).

Theorem 5.1. Under the assumption $\tau \leq \epsilon/2$, $\max \limits_{x \in \overline{\Omega}} |u_0(x)| \leq 1$, and $\|\big(I + \tau \mathcal{A}^{\alpha} \big)^{-1}\|_{\infty} \leq 1$, the numerical solutions obtained by the scheme (4.3) can guarantee the discrete energy property, i.e.,

$\mathcal{E}_{h, n}(U^{n+1}) \leq \mathcal{E}_{h, n}(U^{n}) \qquad \forall n = 0, 1, \ldots,$

(5.5)

where

$\mathcal{E}_{h, n}(U) = \frac{\tau}{4\epsilon} \sum \limits_{i} (U_i^2-1)^2 + \frac{\tau}{2}U^T \mathbf{A}^{\alpha} U.$

(5.6)

Proof. We show the maximum principle property of the first-order scheme (4.3) by induction. First, we have $\|U^0\|_{\infty} \leq 1$ from the assumption. Now, we assume that the result holds for $n = m$, i.e., $\|U^m\|_{\infty} \leq 1$. Next, we show that it is true for $n = m+1$. It follows from the numerical scheme (4.3) that

$\big( I + \tau \mathcal{A}^{\alpha} \big) U^{m+1} = U^m + \frac{\tau}{\epsilon} \big( U^m - (U^m)^3 \big).$

(5.7)

Each element of $U^m + \frac{\tau}{\epsilon} \big(U^m -(U^m)^3 \big)$ is of the form

$g(x) = x + \frac{\tau}{\epsilon} (x-x^3).$

Then, the following statement holds

$g^{\, '}(x) = 1 + \frac{\tau}{\epsilon} (1-3x^2) \geq 0, \qquad \forall x \in [-1, 1]$

provided that $0 \leq \tau \leq \epsilon/2$. It implies that $|g(x)| \leq 1$ for $|x| \leq 1$. As a result, we obtain $\|U^m + \frac{\tau}{\epsilon} \big(U^m -(U^m)^3 \big)\|_{\infty} \leq 1$ if $\|U^m\|_{\infty} \leq 1$. This, together with the assumption, completes the induction. Consequently, we have

$\|U^n\|_{\infty} \leq 1 \qquad \forall n = 0, 1, \ldots.$

(5.8)

Now, using the maximum principle we show the discrete energy decay property of the scheme (4.3). Taking the difference of the discrete energy between two time levels, we obtain

$\begin{eqnarray*} \mathcal{E}_{h, n}(U^{n+1}) - \mathcal{E}_{h, n}(U^{n}) & = & \frac{\tau}{4\epsilon} \sum \limits_{i} \left[\big((U_i^{n+1})^2-1\big)^2-\big((U_i^{n})^2-1\big)^2 \right] \\ && + \frac{\tau}{2} \left( (U^{n+1})^T \mathbf{A}^{\alpha} U^{n+1} - (U^{n})^T \mathbf{A}^{\alpha} U^{n} \right). \end{eqnarray*}$

The application of the $L^2$ inner product of (4.3) with $\big(U^{n+1}-U^n\big)^T$ yield

$\big(U^{n+1}-U^n\big)^T \big(U^{n+1}-U^n\big) + \tau \big( U^{n+1}-U^n\big)^T \mathbf{A}^{\alpha} U^{n+1} + \frac{\tau}{\epsilon} \big( U^{n+1}-U^n\big)^T \left( (U^n)^3 - U^n \right) = 0.$

The following equality holds

$\begin{eqnarray*} \big( U^{n+1}-U^n\big)^T \mathbf{A}^{\alpha} U^{n+1} & = & \frac{1}{2} \left( (U^{n+1})^T \mathbf{A}^{\alpha} U^{n+1} - (U^{n})^T \mathbf{A}^{\alpha} U^{n} \right) \\ &&+ \frac{1}{2} \left( (U^{n+1}-U^n)^T \mathbf{A}^{\alpha} (U^{n+1}-U^n) \right). \end{eqnarray*}$

We note that

$\frac{1}{4} \left( (a^2 -1)^2 - (b^2-1)^2 \right) \leq (b^3-b)(a-b) + (a-b)^2$

for all $a, b \in [-1, 1]$. An application of the above inequality and the property $\|U^n\|_{\infty} \leq 1 \; \; \forall n = 0, 1, \dots$ yields

$\begin{eqnarray*} \mathcal{E}_{h, n}(U^{n+1}) - \mathcal{E}_{h, n}(U^{n}) & \leq & \frac{\tau}{\epsilon} \sum \limits_i \left[\big((U_i^n)^3-U_i^n \big) \big(U_i^{n+1}-U_i^n \big) + \big( U_i^{n+1}-U_i^n \big)^2 \right] \\ && + \frac{\tau}{2} \left( (U^{n+1})^T \mathbf{A}^{\alpha} U^{n+1} - (U^{n})^T \mathbf{A}^{\alpha} U^{n} \right) \\ & = & \big( \frac{\tau}{\epsilon} - 1 \big) \sum \limits_i \big( U_i^{n+1}-U_i^n\big)^2 \\ && - \frac{\tau}{2} \left( (U^{n+1}-U^n)^T \mathbf{A}^{\alpha} (U^{n+1}-U^n) \right). \end{eqnarray*}$

We know that the matrices $\mathbf{L}$ and $\mathbf{M}$ are real and symmetric. Also, $\mathbf{L}$ is nonnegative definite, and $\mathbf{M}$ is positive definite. Moreover, the nonnegative definiteness of $\mathbf{L}$ implies $x^T \mathbf{M}^{-1/2} \mathbf{A} \mathbf{M}^{1/2} x \geq 0$ so that $\mathbf{A}$ is also nonnegative definite.

With the help of the nonnegative definiteness of $\mathbf{A}$ and the assumption $\tau \leq \epsilon/2$, the desired result is obtained.

---

6. Numerical results

In this section, we investigate the performance of our spatial and temporal discretization strategies for the fractional-in-space Allen-Cahn/Cahn-Hilliard equations. To achieve the required accuracy for all examples, $50$ quadrature points are used in the contour integral method described in Section 3. We use piecewise linear polynomials to form the SIPG discretization in space in all numerical experiments. The penalty parameter $\sigma$ in the SIPG discretization is chosen as $\sigma = 6$ on the interior edges $\mathcal{E}^0$ and $\sigma = 12$ on the boundary edges $\mathcal{E}^{\partial}$. If not specified, all examples are implemented on a mesh, constructed by first dividing $\Omega$ into $32 \times 32$ uniform squares and then dividing each square into two triangles.

We note that the reliable performance of the phase-field computations demands long simulation time, it is critical to understand the stability properties of the underlying numerical schemes. However, our simulations end with the small final time since our model problems (1.1) or (1.2) are scaled versions of the original Allen-Cahn/Cahn-Hilliard equations.

---

6.1. Allen-Cahn model

---

6.1.1. Motion of a circle with double-well potential

We first consider a motion of a circle on the domain $\Omega = [0, 1] \times [0, 1]$ with the double-well potential (1.5), taken from ^[77] in order to show the order of convergence. Suppose a radially symmetric initial condition is given as follows:

$u_0(x, y) = \tanh \frac{0.25 - \sqrt{ (x-0.5)^2 + (y-0.5)^2}}{ \sqrt{2 \epsilon}},$

which represents a circle centered at $(0.5, 0.5)$ with a radius $R_0 = 0.25$. It is well known that the solution with the initial condition $u_0(x, y)$ is radially symmetric and the radius of the interfacial circle shrinks by the rate of the curvature of the circle. The rest of problem data are

$\partial \Omega = \Gamma^N, \quad g^N = 0, \quad \epsilon = 0.01.$

Since the fractional-in-space Allen-Cahn equations (1.1) do not have exact solution, we choose a reference numerical solution $\hat u_{h, \tau}$ as the exact solution in order to compare with the corresponding coarse time-or space-stepping approximations. The reference solutions are computed on the mesh, which is constructed by first dividing $s \times s$ uniform squares and then dividing each square into triangles and at the final time $T = 0.0256$. The uniform time-stepping size is chosen as $\tau = T/s$. Here, the reference solution $\hat u_{h, \tau}$ is computed by taking $s$ to be equal to $64$.

Table 1 shows the $L_2$-norm error and convergence rate of various fractional orders, i.e., $\alpha = 1, 0.88, 0.64$, at $T = 0.0256$ by applying the first order IMEX scheme (4.3). The rate of convergence is independent of the fractional order $\alpha$. The results are similar to the observations in ^[59].

---

6.1.2. Dumbbell example with double-well potential

We now consider a dumbbell example, taken from ^[78], with the double-well potential (1.5). The data of problem are

$\Omega = [-1, 1]^2, \quad \partial \Omega = \Gamma^N, \quad g^N = 0, \quad \epsilon = 0.0025, \quad \tau_0 = 5 \times 10^{-5},$

with the following initial condition

$u_0(x, y) = \left\{ \begin{array}{ll} \tanh\Big(\frac{3}{\epsilon} \big((x-0.5)^2 + y^2 -(0.39)^2 \big) \Big),&\hbox{if } x >0.14, \\ \tanh\Big(\frac{3}{\epsilon} \big( y^2 -(0.15)^2 \big) \Big),&\hbox{if } -0.3 \leq x \leq 0.14, \\ \tanh\Big(\frac{3}{\epsilon} \big( (x+0.5)^2 + y^2 -(0.25)^22 \big) \Big),&\hbox{if } x < -0.3. \end{array} \right.$

Figure 2 displays the initial condition $u_0$, which is a dumbbell shape with unequal bells. The snapshots of the solution of the fractional-in-space Allen-Cahn equation in time are displayed for various fractional powers ($\alpha = 1, 0.9, 0.8$) in Figure 3. With standard diffusion, i.e., $\alpha = 1$, we see that the curvature drives towards a circle (constant curvature) in time. The motion of smaller fractional powers is similar, although the rate is slower.

The behaviour of the numerical value of the energy function (1.6) and the adaptive time-step size are displayed in Figure 4 with the tolerance parameter $Tol = 10^{-3}$. The energy function (1.6) decreases in time for all cases. Reducing the fractional power increases the required time to reach the metastable state. For parameters $Tol \in \{5.10^{-3}, 10^{-3}, 5.10^{-4} \}$, the evolution of the length of the time step is shown in Figure 5. The time-step size oscillates for smaller fractional powers, when the tolerance parameter is large. In addition, the number of time steps increases with decreasing the $Tol$ and the fractional power $\alpha$.

The number of time-steps is given in Table 2 with various tolerance parameters. It can be seen that the number of rejected time-steps is increasing for small fractional powers with large tolerance parameter. Table 3 also shows the performance of the adaptive time-step size with respect to the uniform time-step size, i.e., $\tau = 5 \times 10^{-5}$. As expected, the time-step size adaptivity allows us to reduce the computing time compared to the uniform time-step size.

Now, we investigate the coarsening rate of the Allen-Cahn equation with respect to changes of the fractional power. To measure this quantity, we consider the rate of change of the energy density, as done in ^[79,80]. The rate of change of the energy with respect to time is considered as the form of $at^{-b}$, where $t$ represents the time. The coefficients, i.e., $a$ and $b$, are computed by using the following commands in MATLAB^®:

$\mathrm{f = fittype('a*x^ \wedge (-b)'); }$

$\mathrm{cfun = fit(time, energy, f, 'Startpoint', [10, -1/3]); }$

Table 4 shows the values of the coefficients $a$ and $b$ obtained for the various values of the $\alpha$ in the adaptive time refinement. As the $\alpha$ decreases, the values of the coefficient $b$ also decreases, which means that the decay of the energy decreases.

---

6.1.3. Intersection of two dumbbells with double-well potential

We now investigate an intersection of two dumbbells on the Laplacian. The double-well potential (1.5) function is taken. The rest of problem data are

$\Omega = [-1, 1]^2, \quad \partial \Omega = \Gamma^N, \quad g^N = 0, \quad \epsilon = 0.01, \quad \tau_0 = 5 \times 10^{-5},$

with the following initial condition

$u_0(x, y) = u_0^1(x, y)u_0^2(x, y),$

where

$u_0^1(x, y) = \left\{ \begin{array}{ll} \tanh\Big(\frac{3}{\sqrt{\epsilon}} \big((x-0.5)^2 + (y-0.4)^2 -(0.25)^2 \big) \Big), & \hbox{if } x >0.3, \\ \tanh\Big(\frac{3}{\sqrt{\epsilon}} \big( (y-0.4)^2 -(0.15)^2 \big) \Big), & \hbox{if } -0.3 \leq x \leq 0.3, \\ \tanh\Big(\frac{3}{\sqrt{\epsilon}} \big( (x+0.5)^2 + (y-0.4)^2 -(0.25)^2 \big) \Big), & \hbox{if } x < -0.3, \end{array} \right.$

and

$u_0^2(x, y) = \left\{ \begin{array}{ll} \tanh\Big(\frac{3}{\sqrt{\epsilon}} \big(x^2 + (y-0.6)^2 -(0.25)^2 \big) \Big), & \hbox{if } y >0.4, \\ \tanh\Big(\frac{3}{\sqrt{\epsilon}} \big( x^2 -(0.15)^2 \big) \Big), & \hbox{if } -0.4 \leq y \leq 0.4, \\ \tanh\Big(\frac{3}{\sqrt{\epsilon}} \big( x^2 + (y-0.6)^2 -(0.25)^2 \big) \Big), & \hbox{if } y < -0.4. \end{array} \right.$

The initial function $u_0$, which is an intersection of two dumbbells, is shown in Figure 7. For various fractional powers $(\alpha = 1, 0.9, 0.8, 0.7, 0.6$), the snapshots of the solutions in time are displayed in Figure 6. As the previous example, reducing of fractional power decreases the rate of motion of initial curvature.

Figure 8 illustrates the behaviour of the numerical energy function (1.6) and the adaptive time-step size. The evolution of the length of the time-step for $\alpha = 0.8, 0.7, 0.6$ is shown in Figure 9 for parameters $Tol \in \{5.10^{-3}, 10^{-3}, 5.10^{-4} \}$. It is observed that decreasing the tolerance parameter makes the motion of time-step size smoother.

Table 5 shows the number of time-steps for $\alpha = 1, 0.9, 0.8, 0.7, 0.6$ with various tolerance parameters. The number of time steps increases with decreasing the parameter $Tol$ and the fractional power $\alpha$. Further, the number of the adaptive and uniform time-steps are displayed in Table 6. It can be seen that the time-step size adaptivity allows us to reduce the computation time by factors of hundreds compared to the uniform time-step size.

As the previous example, Table 7 shows the values of the coefficients $a$ and $b$ in the form $at^{-b}$. Similarly, the value of the $b$ decreases, when $\alpha$ decreases.

---

6.1.4. Spinodal decomposition with logarithmic free energy

We now consider a test example with the logarithmic free energy (1.4). The initial condition is a random state by randomly distributing numbers from $-0.01$ to $0.01$. The rest of problem data are

$\Omega = [0, 2 \pi]^2, \;\; \partial \Omega = \Gamma^D, \;\; g^D = 0, \;\; \tau_0 = 5 \times 10^{-5}, \;\; \theta = 0.1, \;\; \theta_c = 0.2.$

In this example, we investigate the effect of fractional power when a spinodal decomposition is considered. The initial state is well mixed, see Figure 10. The snapshots of the phase evolution for various values of fractional power ($\alpha = 1, 0.8, 0.6$) are illustrated in Figure 11 with $\epsilon = 10^{-3}$.

Early stages of phase transition yield a rapid movement to bulk regions for $\alpha = 1$. However, smaller fractional powers lead much more heterogeneous phase structures with smaller bulk regions. Figure 12 also shows the snapshots of phase evaluation at $t = 0.016$ with $\epsilon = 10^{-4}$.

The behaviour of the numerical energy function (1.6) and the adaptive time-step size versus time are displayed in Figure 13 for $\epsilon = 10^{-4}$. The numerical energy decrease is observed for all the cases. Lastly, Figure 14 shows the evolution of the time-step size for various tolerance parameters.

---

6.2. Cahn-Hilliard model

Finally, we consider a Cahn-Hilliard system with double-well potential and the following initial condition:

$\begin{align*} u_0(x, y)& = \tanh \left( \frac{1}{\sqrt{2}\epsilon} \left( \min \{\sqrt{(x+0.3)^2 + y^2}-0.2, \sqrt{(x-0.3)^2 + y^2}-0.2, \right. \right. \\ & \hspace{30mm}\left. \left. \sqrt{x^2 + (y+0.3)^2 }-0.2, \sqrt{x^2 + (y-0.3)^2 + x^2} \} \right) \right). \end{align*}$

The rest of problem data are

$\Omega = [-1, 1]^2, \quad \partial \Omega = \Gamma^N, \quad \epsilon = 0.1, \quad \tau_0 = 10^{-5}.$

Figure 15 plots the change of the discrete energy in time, which decreases as predicted, by using the semi-implicit method in time with the fixed time step. As the Allen-Cahn examples, the rate of the motion is slower for the small number of the fractional powers.

Further, Table 8 shows the values of the coefficients $a$ and $b$ in the form $at^{-b}$ with the fixed time step. As the previous examples, the smaller the value of $\alpha$, the smaller the energy decay does. The total mass $m = \int_{\Omega} u \; dx$ equals to a constant value $2.8754$.

---

7. Conclusions

In this paper we have investigated the numerical solutions of the fractional-in-space phase-field equations, discretized by the symmetric interior penalty Galerkin (SIPG) method in space and an implicit-explicit (IMEX) method in time. The contour integral method (CIM) has been used to compute the fractional power of a matrix times a vector. To reduce computation time, an adaptive time-step size method is proposed. The numerical results of the fractional-in-space phase-fields equations show that such a kind of modelling can aid to understand the effects of spatial heterogeneity. Although, the ideas expressed in this paper have applicability in this setting, the numerical approximations of the Riesz derivatives, discretized by discontinuous Galerkin methods, should also be considered in the more general framework.

---

8. Acknowledgements

This work was supported by the Forschungszentrum Dynamische Systeme (CDS): Biosystemtechnik, Otto-von-Guericke-Universität Magdeburg.

---