# LOWEST-ORDER NONSTANDARD FINITE ELEMENT METHODS FOR TIME-FRACTIONAL
BIHARMONIC PROBLEM

Shantiram Mahata1<EMAIL_ADDRESS>Neela Nataraj1
<EMAIL_ADDRESS>and Jean-Pierre Raymond2
1 Department of Mathematics, Indian Institute of Technology Bombay, Mumbai,
Maharashtra - 400076, India
2Université Paul Sabatier Toulouse III & CNRS, Institut de Mathématiques,
31062 Toulouse Cedex 9, France; Visiting Professor, Indian Institute of
Technology Bombay, Mumbai, Maharashtra - 400076, India<EMAIL_ADDRESS>toulouse.fr

###### Abstract

In this work, we consider an initial-boundary value problem for a time-
fractional biharmonic equation in a bounded convex polygonal domain in
$\mathbb{R}^{2}$ with clamped boundary conditions. After establishing the
well-posedness, we focus on some regularity results of the solution with
respect to the regularity of the problem data. The spatially semidiscrete
scheme covers several popular lowest-order piecewise-quadratic finite element
schemes, namely, Morley, discontinuous Galerkin, and $C^{0}$ interior penalty
methods, and includes both smooth and nonsmooth initial data. Optimal order
error bounds with respect to the regularity assumptions on the data are proved
for both homogeneous and nonhomogeneous problems. The numerical experiments
validate the theoretical convergence rate results.

Key words: Time-fractional equation, biharmonic, lowest-order FEM, smooth and
nonsmooth initial data, semidiscrete, error estimates

AMS subject classifications: 35R11, 65M15, 65M60

##  1 Introduction

Let $\Omega$ be a bounded convex polygonal domain in $\mathbb{R}^{2}$ with
boundary $\partial\Omega$ and $T$ be a fixed positive real number. For
$0<\alpha<1$, consider the following initial-boundary value problem for time-
fractional biharmonic equation that seeks $u(t)\in H_{0}^{2}(\Omega)$ such
that

$\displaystyle\begin{aligned} \partial_{t}^{\alpha}u(t)+Au(t)&=f(t)\text{ in
}\Omega,\;0<t\leq T,\\\ u(0)&=u_{0}\text{ in }\Omega,\\\ \end{aligned}$ (1.1)

where $A$ is the unbounded operator in $L^{2}(\Omega)$ corresponding to the
biharmonic operator $\Delta^{2}$ with clamped boundary conditions. For an
absolutely continuous real function $u$ over $[0,T]$,
$\partial_{t}^{\alpha}u(t)$ denotes the Caputo fractional derivative defined
by

$\displaystyle\partial_{t}^{\alpha}u(t)=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}(t-s)^{-\alpha}u^{\prime}(s)\,ds=\mathcal{I}^{1-\alpha}u^{\prime}(t),\;\;0<\alpha<1,$
(1.2)

where $u^{\prime}(t)\mathrel{\mathop{\mathchar 58\relax}}=\partial_{t}u(t)$
and $\mathcal{I}^{\beta}$, $0<\beta<\infty$, denotes the Riemann-Liouville
fractional integral defined by

$\displaystyle\mathcal{I}^{\beta}v(t)=\int_{0}^{t}\kappa_{\beta}(t-s)v(s)\,ds,\quad\kappa_{\beta}(t)\mathrel{\mathop{\mathchar
58\relax}}=t^{\beta-1}/\Gamma(\beta),$ (1.3)

with the gamma function
$\Gamma(z)=\int_{0}^{\infty}e^{-t}t^{z-1}\;dt,\;\text{Re}(z)>0$. In Definition
2.1, we will use (1.2) for absolutely continuous functions with values in a
Hilbert space $V^{*}$ introduced later.

It is well-known that fractional order models describe physical phenomena more
accurately compared to the usual integer-order models (cf. [podlubny]). The
equation in model (1.1) represents a particular case of the time-fractional
Cahn-Hilliard equation [huangstynes21, maskarikaraa22].

The main goal of this work is to study the convergence analysis of
semidiscrete approximation with lowest-order nonstandard finite element
methods (FEMs). Prior to the development of numerical methods, we first prove
the well-posedness of problem (1.1).

To present our contributions, we introduce below some standard notations. For
any real $r$, $H^{r}(\Omega)$ denotes the standard Sobolev space associated
with the Sobolev-Slobodeckii semi-norm $|\cdot|_{H^{r}(\Omega)}$ (cf.
[grisvard92]). Let $H^{r}_{0}(\Omega)$ denote the closure of $D(\Omega)$ in
$H^{r}(\Omega)$, $V\mathrel{\mathop{\mathchar 58\relax}}=H^{2}_{0}(\Omega)$,
and let $H^{-r}(\Omega)=(H^{r}_{0}(\Omega))^{*}$ be the dual of
$H^{r}_{0}(\Omega)$. The notations $(\cdot,\cdot)$ and $\|\cdot\|$ denote the
$L^{2}$ scalar product and $L^{2}$ norm in ${\Omega}$. For a Hilbert space $H$
with norm $\|\cdot\|_{H}$, let $L^{p}(0,T;H)$, $1\leq p<\infty$ denote the
standard Bochner space equipped with the norm
$\|g\|_{L^{p}(0,T;H)}=(\int_{0}^{T}\|g(t)\|_{H}^{p}\,dt)^{1/p}$. If
$p=\infty$, the space $L^{\infty}(0,T;H)$ has the norm
$\|g\|_{L^{\infty}(0,T;H)}=\inf\\{C\mathrel{\mathop{\mathchar
58\relax}}\|g(t)\|_{H}\leq C\text{ almost everywhere on }(0,T)\\}$. In
addition, for $m\geq 0$, $1\leq p\leq\infty$, $W^{m,p}(0,T;H)$
$\mathrel{\mathop{\mathchar 58\relax}}=\\{g\mathrel{\mathop{\mathchar
58\relax}}(0,T)\rightarrow H\text{ such that }\frac{d^{j}g}{dt^{j}}\in
L^{p}(0,T;H)\text{ for }0\leq j\leq m\\}$ is equipped with the norm
$\|g\|_{W^{m,p}(0,T;H)}=\sum_{j=0}^{m}\|\frac{d^{j}g}{dt^{j}}\|_{L^{p}(0,T;H)}$
(cf. [evans10]). Hereafter, $\gamma_{0}\in(1/2,1]$ is the elliptic regularity
index of the biharmonic problem introduced in Section 2.1. Throughout the
paper, $a\lesssim b$ denotes $a\leq Cb$, where $C$ is a positive generic
constant that may depend on fractional order $\alpha$ and the final time $T$
but is independent of the spatial mesh size.

The main contributions of this work are summarized below.

  * 1.

We study the existence, uniqueness, and stability of the solution of time-
fractional biharmonic problem (1.1) with clamped boundary conditions. In
Theorem 2.4, we prove that when initial data $u_{0}$ is simply an element of
$L^{2}(\Omega)$ and source function $f$ belongs to $W^{1,1}(0,T;V^{*})$,
problem (1.1) admits a unique weak solution with

$\displaystyle\|u\|_{L^{1}(0,T;V)}+\|u^{\prime}\|_{L^{1}(0,T;V^{*})}\lesssim\|u_{0}\|+\|f\|_{W^{1,1}(0,T;V^{*})}.$

In addition, we also derive regularity results for the solutions of both
homogeneous and nonhomogeneous problems useful for the error analysis.

  * 2.

In (3.4), we introduce a novel Ritz projection
$\mathcal{R}_{h}\mathrel{\mathop{\mathchar 58\relax}}V\rightarrow V_{h}$,
where $V_{h}$ denotes the space of lowest-order piecewise-quadratic
polynomials for (1.1) with clamped boundary conditions in the proposed
semidiscrete schemes and includes popular Morley, discontinuous Galerkin (dG),
and $C^{0}$ interior penalty ($C^{0}$IP) methods. We establish the quasi-
optimal approximation property displayed below in Lemma 3.4.

$\|{v-\mathcal{R}_{h}v}\|+h\mathinner{\lVert
v-\mathcal{R}_{h}v\rVert}_{h}\lesssim h^{1+\gamma}\mathinner{\lVert
v\rVert}_{H^{2+\gamma}(\Omega)},$

for all $\gamma\in[0,\gamma_{0}]$ and for all $v\in H^{2+\gamma}(\Omega)$.
These properties are crucial to obtaining optimal order error bounds for both
smooth and nonsmooth initial data.

  * 3.

We use an energy argument in a unified framework for the nonstandard FEM
analysis with lowest-order piecewise-quadratic polynomials for (1.1). Since
the solution of model (1.1) reflects a singular behavior near $t=0$ (see below
Theorems 2.6 and 2.8), a straightforward analysis of numerical scheme is
problematic. This is addressed by multiplication by weights of type $t^{j}$,
$j=1,2$ in the semidiscrete scheme. For initial data in $L^{2}(\Omega)$ and a
source term equal to zero, in Theorem 4.2 $(i)$ we obtain the following error
bound in the $L^{2}(\Omega)$ and energy norms for each fixed time $t\in(0,T]$,

$\displaystyle\mathinner{\lVert u(t)-u_{h}(t)\rVert}$
$\displaystyle\lesssim\big{(}h^{1+\gamma_{0}}t^{-\alpha}+h^{2}t^{-(1+\alpha)/2}\big{)}\mathinner{\lVert
u_{0}\rVert},$ $\displaystyle\mathinner{\lVert u(t)-u_{h}(t)\rVert}_{h}$
$\displaystyle\lesssim h^{\gamma_{0}}t^{-\alpha}\mathinner{\lVert
u_{0}\rVert}+\big{(}h^{1+\gamma_{0}}t^{-3\alpha/2}+h^{2}t^{-(\alpha+1/2)}\big{)}\mathinner{\lVert
u_{0}\rVert},$

where $h$ is the mesh size of the FEM. Similar error bounds are also proved
for the nonhomogeneous problem with an appropriate assumption on $f$ and
$u_{0}\in D(A)$ in Theorem 4.2 $(ii)$ for $t\in(0,T]$.

  * 4.

Numerical experiments that validate the theoretical results are presented for
examples with smooth and nonsmooth initial data.

There are technical difficulties that arise in the error analysis because (i)
the space discretization is not conforming and (ii) the error estimates are
established under lower regularity of the exact solution. To overcome these
challenges, we introduce a novel Ritz projection $\mathcal{R}_{h}$ using a
smoother (or companion) operator $Q$ in (3.4). This leads to an error equation
in (4.13) that is different from those obtained for the conforming second-
order problems [mustapha18, karaamustaphapani18]. The idea is also completely
new in comparison to the available works for time fractional fourth-order
problems available in the literature.

The article is organized as follows:

  * -

In Section 2, we introduce some notations, establish the well-posedness of
problem (1.1), and prove regularity results of the exact solution.

  * -

This is followed by a discussion on the semidiscrete schemes and error
analysis in Sections 3 and 4. Optimal order error bounds are proved for both
smooth and nonsmooth initial data.

  * -

Finally, Section 5 presents the results of the numerical implementations that
validate the theoretical estimates.

The semidiscrete convergence analysis in this work is inspired by the
technique in [karaamustaphapani18, mustapha18, smrksenergy] for time-
fractional parabolic problems with and without memory. These works discuss
only second-order problems and cover only the conforming case, where the
discrete trial space is a finite-dimensional subspace of $H_{0}^{1}(\Omega)$.
In contrast, the analysis of conforming FEMs for fourth-order problems demands
$C^{1}$ continuity and requires a high polynomial degree, which imposes many
conditions on the vertices and edges of an element, and hence, is
prohibitively expensive. Most of the available literature for (1.1) discuss
simply supported boundary conditions ($u=\Delta u=0$) that enables a mixed
formulation via an auxiliary variable $v=\Delta u$ and subsequently, the
reduction of the problem to a system of two second-order equations.

Below, we compare our findings with those already known from the literature
and review the existing results for (1.1).

  * -

The well-posedness result in Point (1) above is inspired by the work in [sy11]
for the second-order fractional diffusion-wave equation. For (1.1), the well-
posedness result is summarized in [lietal2022_cnc, Appendix]. However, in
[lietal2022_cnc, (A.1) and (A.2)], the well-posedness and regularity results
are stated for higher regularity assumptions on the problem data, for example,
initial data $u_{0}\in D(\Delta^{j+2})$ and source function $f$ such that
$\frac{\partial^{l}f(t)}{\partial t^{l}}\in D(\Delta^{j})$, where
$j\in\mathbb{N}$ and $l=0,1,2,3$. These studies to (1.1) in our work are done
comparatively less smoother assumptions on the problem data (cf. Theorems 2.4,
2.6, and 2.8).

  * -

The introduction of the Ritz projection for lowest-order piecewise-quadratic
based FEM on $H_{0}^{2}(\Omega)\cap H^{2+\gamma}(\Omega)$ and its quasi-
optimal approximation properties on $H^{2+\gamma}(\Omega)$ in $L^{2}(\Omega)$
and energy norms are new. The idea of this operator is borrowed from [cn2022]
for the analysis of FEMs for biharmonic plates. In [lietal2022_cnc], the
authors have introduced the Ritz projection on $H_{0}^{2}(\Omega)\cap
H^{j+1}(\Omega)$, $j\geq 2$ for the analysis of virtual element methods, and
have proved optimal projection errors which demand at least $H^{3}(\Omega)$
regularity of the continuous solution $u(t)$. We also have to mention that the
idea of using a Ritz projection is suggested, without proof, in
[danumjayaetal2021, Remark 3.1], for the Morley finite element approximation
of fourth-order nonlinear reaction-diffusion problems.

  * -

As far as we know, the semidiscrete error analysis with $u_{0}\in
L^{2}(\Omega)$ is a totally new result for the time-fractional biharmonic
problem with clamped boundary conditions.

  * -

The work in [cuifd2021] deals with finite difference method with clamped plate
problem.

Several numerical methods to (1.1) and its nonlinear variants with simply
supported  boundary conditions (that is, $u=\Delta u=0$) are studied in
literature with stability and convergence analysis. These methods include, for
example, Adomian decomposition method [jafarietal08], finite difference method
[zhang_pu2017], weak Galerkin method [yazdaniweak2022], mixed FEMs [liuetal14,
liuetal15, liuetal2015, huangstynes21, huangetal22], local discontinuous
Galerkin methods [weihe14, duetal2017dg], Petrov-Galerkin method [pg2020], and
spline-based methods [hadhoud19, zhang2020]. The FEM based works in
[liuetal14, liuetal15, liuetal2015, huangstynes21, huangetal22] for fourth-
order problems are studied for simply supported boundary conditions. An
auxiliary variable $v=\Delta u$ is introduced for the mixed formulation, and
subsequently, the problem reduces to a system of two second-order equations
with boundary conditions $u=v=0$.

##  2 Well-posedness and regularity

We introduce function spaces and norms, and present some useful properties of
the Mittag-Leffler functions in Subsection 2.1. The well-posedness results for
(1.1) are discussed in Subsection 2.2. In Subsection 2.3, regularity results
for the solutions of both homogeneous and nonhomogeneous problems are derived.

###  2.1 Preliminaries

### The elliptic operator $A$

To study the well-posedness results, we start with the solution
representation. The results are based on the eigenfunction expansion of the
corresponding elliptic operator associated with homogeneous boundary
conditions. Consider the following eigenvalue problem

$\displaystyle\begin{aligned} &\Delta^{2}\phi=\mu\phi\text{ in
}\Omega,\;\phi=\frac{\partial\phi}{\partial\nu}=0\text{ on
}\partial\Omega.\end{aligned}$ (2.1)

Here $\nu$ is the outward unit normal to $\partial\Omega$. It is well-known
[babuskaosborn91, p. 761] that there exists a family
$(\mu_{j},\phi_{j})_{j=1}^{\infty}$ such that $(\mu_{j},\phi_{j})$ is solution
of (2.1), $0<\mu_{1}\leq\mu_{2}\leq\ldots$, $\mu_{j}\rightarrow\infty$ as
$j\rightarrow\infty$, and the corresponding family of eigenfunctions
$\\{\phi_{j}\\}_{j=1}^{\infty}$ is an orthonormal basis of $L^{2}(\Omega)$.
For further development of this work, define the unbounded operator $(A,D(A))$
in $L^{2}(\Omega)$ by

$D(A)=\\{\phi\in H_{0}^{2}(\Omega)\mid\Delta^{2}\phi\in L^{2}(\Omega)\\}\mbox{
and }A\phi=\Delta^{2}\phi\text{ for all }\phi\in D(A).$

It is known that [blumrannacher80]

$D(A)\subset V\cap H^{2+\gamma_{0}}(\Omega)\mbox{ with
}V\mathrel{\mathop{\mathchar 58\relax}}=H_{0}^{2}(\Omega)\text{ and
}\gamma_{0}\in(1/2,1],$ (2.2)

where $\gamma_{0}$ is the so-called elliptic regularity index of the
biharmonic operator with clamped boundary conditions on $\partial\Omega$.
Define the domain of the fractional power $A^{r},\,r\in\mathbb{R}^{+}$
[pazybook, Chapter 2], of the operator $A$, by

$D(A^{r})=\Big{\\{}v=\sum_{j=1}^{\infty}v_{j}\phi_{j}\mid
v_{j}\in{\mathbb{R}},\
\sum_{j=1}^{\infty}\mu_{j}^{2r}|v_{j}|^{2}<\infty\Big{\\}},\;A^{r}v=\sum_{j=1}^{\infty}\mu_{j}^{r}v_{j}\phi_{j}=\sum_{j=1}^{\infty}\mu_{j}^{r}(v,\phi_{j})\phi_{j}.$
(2.3)

The space $D(A^{r})$ is equipped with the norm

$\displaystyle\mathinner{\\!\left\lVert
v\right\rVert}_{D(A^{r})}=\bigg{(}\sum_{j=1}^{\infty}\mu_{j}^{2r}|(v,\phi_{j})|^{2}\bigg{)}^{1/2}.$
(2.4)

The orthonormality of the basis $\\{\phi_{j}\\}_{j=1}^{\infty}$ is used above
and throughout in this section. In particular, for $r=0$, we obtain
$\mathinner{\lVert\cdot\rVert}_{D(A^{0})}=\mathinner{\lVert\cdot\rVert}_{L^{2}(\Omega)}=\mathinner{\lVert\cdot\rVert}$.
When $r=1/2$, we have $D(A^{1/2})=V$. From [Bensoussan-etal-2007, Part II,
Chapter 1, Proposition 6.1], it follows that

$D(A^{r})=[L^{2}(\Omega),D(A)]_{r}\text{ for all }r\in[0,1].$

Therefore, with (2.2) and interpolation, we can prove that

$D(A^{r})\subset V\cap H^{2+\gamma_{0}(2r-1)}(\Omega)\text{ for all
}r\in[1/2,1].$

Now, for $r>0$, let $D(A^{-r})\mathrel{\mathop{\mathchar
58\relax}}=(D(A^{r}))^{*}$ be the dual of $D(A^{r})$, consisting of all
bounded linear functionals on $D(A^{r})$. Identifying the dual of
$L^{2}(\Omega)$ with itself, we write $D(A^{r})\subset L^{2}(\Omega)\subset
D(A^{-r})$. The space $D(A^{-r})$ is a Hilbert space with the norm

$\displaystyle\mathinner{\\!\left\lVert
g\right\rVert}_{D(A^{-r})}=\bigg{(}\sum_{j=1}^{\infty}\frac{1}{\mu_{j}^{2r}}|\langle
g,\phi_{j}\rangle|^{2}\bigg{)}^{1/2},$ (2.5)

where for $g\in D(A^{-r})$ and $\phi\in D(A^{r})$, the symbol $\langle
g,\phi\rangle$ stands for the duality pairing of $g$ with $\phi$. Further,
with $g\in L^{2}(\Omega)$ and $\phi\in D(A^{r})$, we write $\langle
g,\phi\rangle=(g,\phi)$ (cf. [brezis, Chapter V]).

Let $a(\cdot,\cdot)\mathrel{\mathop{\mathchar 58\relax}}V\times
V\rightarrow\mathbb{R}$ be the bilinear form associated with $\Delta^{2}$
defined by

$\displaystyle\begin{aligned}
a(v,w)&=\int_{\Omega}D^{2}v\mathrel{\mathop{\mathchar
58\relax}}D^{2}w\,dx=\int_{\Omega}\sum_{i,j=1}^{2}\frac{\partial^{2}v}{\partial
x_{i}\partial x_{j}}\frac{\partial^{2}w}{\partial x_{i}\partial
x_{j}}\,dx,\end{aligned}$

where $D^{2}v$ denotes the Hessian matrix of $v$. Then, we can verify that

$\displaystyle(Av,w)=a(v,w)\text{ for all }v\in D(A)\text{ and for all }w\in
V.$

In the convergence analysis of the semidiscrete scheme, we need the following
version of Gronwall’s lemma (cf. [chenshih98, Lemma 2.1]).

######  Lemma 2.1 (Gronwall’s lemma).

Assume that $\phi,\psi$, and $\chi$ are three non-negative integrable
functions on $[0,T]$. If $\phi$ satisfies

$\displaystyle\phi(t)\leq\psi(t)+\int_{0}^{t}\chi(s)\phi(s)\,ds\text{ for
}t\in(0,T),$

then,

$\displaystyle\phi(t)\leq\psi(t)+\int_{0}^{t}\psi(s)\chi(s)e^{\int_{s}^{t}\chi(\tau)\,d\tau}\,ds\text{
for }t\in(0,T).$

### Mittag-Leffler functions and their properties

The Mittag-Leffler functions play a very important role in the analysis of
fractional differential equations. The two-parameter Mittag-Leffler function
$E_{\alpha,\beta}(\cdot)$ is defined [kilbasbook, p. 42] by

$\displaystyle
E_{\alpha,\beta}(z)=\sum_{j=0}^{\infty}\frac{z^{j}}{\Gamma(j\alpha+\beta)},\;\;z\in\mathbb{C},\;\;\alpha>0,\,\beta\in\mathbb{R}.$

It is a generalization of the exponential function $e^{z}$ in the sense that
$E_{1,1}(z)=e^{z}$. Further, we can directly verify that it is an entire
function in $z$. The following properties of the Mittag-Leffler function
$E_{\alpha,\beta}(\cdot)$ [kilbasbook, (1.8.28)] or [podlubny, Theorem 1.4],
and [sy11, Lemmas 3.2 and 3.3] are essential in Subsections 2.2 and 2.3.

######  Lemma 2.2 (Properties of Mittag-Leffler functions).

Let $r$ be a real number with $\frac{\alpha\pi}{2}<r<\min\\{\pi,\alpha\pi\\}$,
where $\alpha\in(0,2)$. Then for any real number $\beta$ and $r\leq|\arg
z|\leq\pi$,

$\displaystyle|E_{{\alpha},\beta}(z)|\lesssim\begin{cases}\displaystyle{\frac{1}{1+|z|^{2}}},\;\;(\beta-\alpha)\in\mathbb{Z}^{-}\cup\\{0\\},\\\
\displaystyle{\frac{1}{1+|z|}},\;\;\text{otherwise}.\end{cases}$ (2.6)

Further, for $\alpha,\mu>0$ and $k\in\mathbb{N}$,

$\displaystyle\begin{aligned} \frac{d^{k}}{dt^{k}}E_{\alpha,1}(-\mu
t^{\alpha})&=-\mu t^{\alpha-k}E_{\alpha,\alpha-k+1}(-\mu
t^{\alpha}),\;\;t\in(0,T],\\\ E_{\alpha,\alpha}(-\mu)&\geq
0,\;\;\alpha\in(0,1),\;\mu\geq 0.\end{aligned}$ (2.7)

###  2.2 Well-posedness

######  Definition 2.1 (Weak solution).

Let us assume that $u_{0}\in L^{2}(\Omega)$ and $f\in W^{1,1}(0,T;V^{*})$. A
function $u\in L^{1}(0,T;V)\cap W^{1,1}(0,T;V^{*})$ is called a weak solution
to (1.1) if the following identities hold.

$\displaystyle\partial_{t}^{\alpha}\langle u(t),v\rangle+a(u(t),v)=\langle
f(t),v\rangle\text{ in }L^{1}(0,T)\text{ for all }v\in V,$ (2.8) and
$\displaystyle\lim_{t\rightarrow 0}\mathinner{\lVert
u(t)-u_{0}\rVert}_{V^{*}}=0.$

######  Remark 2.3.

Since we have

$\displaystyle\|\partial_{t}^{\alpha}\langle u(t),v\rangle\|_{L^{1}(0,T)}$
$\displaystyle=\|\int_{0}^{t}\frac{(t-s)^{-\alpha}}{\Gamma(1-\alpha)}\langle
u^{\prime}(s),v\rangle\,ds\|_{L^{1}(0,T)}\lesssim\|\langle
u^{\prime}(\cdot),v\rangle\|_{L^{1}(0,T)}$
$\displaystyle\lesssim\|u^{\prime}\|_{L^{1}(0,T;V^{*})}\|v\|_{V}\text{ for all
}v\in V,$

we can define $\partial_{t}^{\alpha}u$ in ${L^{1}(0,T;V^{*})}$ by the equation

$\displaystyle\langle\partial_{t}^{\alpha}u(t),v\rangle\mathrel{\mathop{\mathchar
58\relax}}=\partial_{t}^{\alpha}\langle u(t),v\rangle\text{ for all }v\in
V\text{ and all }t\in(0,T].$

######  Theorem 2.4 (Well-posedness of (1.1)).

Assume that $u_{0}\in L^{2}(\Omega)$ and $f\in W^{1,1}(0,T;V^{*})$. Then,
(1.1) admits a unique weak solution determined by

$\displaystyle u(t)=E(t)u_{0}+\int_{0}^{t}F(t-s)f(s)\,ds,$ (2.9)

with the operators $E(t)\in{\mathcal{L}}(L^{2}(\Omega))$ and
$F(t)\in{\mathcal{L}}{(V^{*})}$ defined by

$\displaystyle\begin{aligned}
E(t)u_{0}&=\sum_{j=1}^{\infty}(u_{0},\phi_{j})E_{\alpha,1}(-\mu_{j}t^{\alpha})\phi_{j}(x),\\\
F(t)f(\cdot)&=\sum_{j=1}^{\infty}\langle f(\cdot),\phi_{j}\rangle
t^{\alpha-1}E_{\alpha,\alpha}(-\mu_{j}t^{\alpha})\phi_{j}(x),\;\;t\in(0,T].\end{aligned}$
(2.10)

Moreover, we have

$\displaystyle\|u\|_{L^{1}(0,T;V)}+\|u^{\prime}\|_{L^{1}(0,T;V^{*})}\lesssim\|u_{0}\|+\|f\|_{W^{1,1}(0,T;V^{*})}.$
(2.11)

The following stability properties of solution operator $F(t)$ of the
nonhomogeneous problem are essential in the proof of Theorems 2.4 and 2.8. For
related results of second-order time-fractional differential equations, we
refer to [jlpz15, Lemma 2.2].

######  Lemma 2.5 (Stability).

For $q\geq-1$, $t\in(0,T]$, and $v\in D(A^{q/2})$, we have

$\displaystyle\mathinner{\lVert F(t)v\rVert}_{D(A^{p/2})}\lesssim
t^{-1+\big{(}1+\frac{(q-p)}{2}\big{)}\alpha}\mathinner{\lVert
v\rVert}_{D(A^{q/2})},\;\;0\leq p-q\leq 4.$

###### Proof.

In view of (2.10), (2.4) (or (2.5)), and (2.6), we obtain

$\displaystyle\mathinner{\lVert F(t)v\rVert}_{D(A^{p/2})}^{2}$
$\displaystyle=\sum_{j=1}^{\infty}\mu_{j}^{p}t^{2\alpha-2}|E_{\alpha,\alpha}(-\mu_{j}t^{\alpha})|^{2}|\langle
v,\phi_{j}\rangle|^{2}\lesssim
t^{2\alpha-2}\sum_{j=1}^{\infty}\frac{\mu_{j}^{p}}{(1+(\mu_{j}t^{\alpha})^{2})^{2}}|\langle
v,\phi_{j}\rangle|^{2}$
$\displaystyle=t^{2\alpha-2+(q-p)\alpha}\sum_{j=1}^{\infty}\frac{(\mu_{j}t^{\alpha})^{p-q}}{(1+(\mu_{j}t^{\alpha})^{2})^{2}}\mu_{j}^{q}|\langle
v,\phi_{j}\rangle|^{2}\leq t^{2\alpha-2+(q-p)\alpha}\mathinner{\lVert
v\rVert}_{D(A^{q/2})}^{2},$

where we have used
$\frac{(\mu_{j}t^{\alpha})^{p-q}}{(1+(\mu_{j}t^{\alpha})^{2})^{2}}\leq 1$ for
$0\leq p-q\leq 4$ and $j\in\mathbb{N}$. This concludes the proof. ∎

Proof of Theorem 2.4. Step 1 (Representation for solution). If the function
$u$ is a solution to (1.1) in the sense of Definition 2.1, then for all
$j\in{\mathbb{N}}$, $(u(\cdot),\phi_{j})$ satisfies

$\displaystyle\begin{aligned} \partial_{t}^{\alpha}\langle
u(t),\phi_{j}\rangle+\mu_{j}(u(t),\phi_{j})&=\langle
f(t),\phi_{j}\rangle\text{ for all }t\in(0,T],\\\ \lim_{t\rightarrow
0}(u(t),\phi_{j})&=(u_{0},\phi_{j}).\end{aligned}$ (2.12)

It is well-known that for all $j\in{\mathbb{N}}$ [gorenflomainardi97, p. 242]
(see also [podlubny, kilbasbook, sy11]), (2.12) has a unique solution
$(u(\cdot),\phi_{j})\in C([0,T])$ defined by

$(u(t),\phi_{j})=(u_{0},\phi_{j})E_{\alpha,1}(-\mu_{j}t^{\alpha})+\int_{0}^{t}\langle
f(s),\phi_{j}\rangle(t-s)^{\alpha-1}E_{\alpha,\alpha}(-\mu_{j}(t-s)^{\alpha})\,ds.$

Let us set

$\displaystyle\displaystyle u(t)$
$\displaystyle=\sum_{j=1}^{\infty}(u(t),\phi_{j})\phi_{j}(x),\;\displaystyle
v(t)=\sum_{j=1}^{\infty}(u_{0},\phi_{j})E_{\alpha,1}(-\mu_{j}t^{\alpha})\phi_{j}(x),$
$\displaystyle\text{and }\displaystyle w(t)$
$\displaystyle=\sum_{j=1}^{\infty}\Big{(}\int_{0}^{t}\langle
f(s),\phi_{j}\rangle(t-s)^{\alpha-1}E_{\alpha,\alpha}(-\mu_{j}(t-s)^{\alpha})\,ds\Big{)}\phi_{j}(x).$

In the next steps, we show that $u(t)=v(t)+w(t)$ is indeed a weak solution.

Step 2 (Proofs of $v,w\in L^{1}(0,T;V)$). From (2.3), the definition of
$v(t)$, and the second inequality in (2.6) we obtain

$\displaystyle\mathinner{\lVert
A^{1/2}v(t)\rVert}^{2}=\sum_{j=1}^{\infty}|\mu_{j}^{1/2}(u_{0},\phi_{j})E_{\alpha,1}(-\mu_{j}t^{\alpha})|^{2}\lesssim
t^{-\alpha}\sum_{j=1}^{\infty}|(u_{0},\phi_{j})|^{2}\frac{\mu_{j}t^{\alpha}}{(1+\mu_{j}t^{\alpha})^{2}}\leq
t^{-\alpha}\mathinner{\lVert u_{0}\rVert}^{2},$

with $\frac{\mu_{j}t^{\alpha}}{(1+\mu_{j}t^{\alpha})^{2}}\leq 1$ and (2.4) in
the last step above. This and $D(A^{1/2})=V$ show that $\mathinner{\lVert
v(t)\rVert}_{V}\lesssim t^{-\alpha/2}\|u_{0}\|$ and hence $v\in L^{1}(0,T;V)$.

Next, to show $w\in L^{1}(0,T;V)$, we first invoke the last inequality and
then the first identity in (2.7) to arrive at

$\displaystyle\begin{aligned}
\int_{0}^{\eta}|t^{\alpha-1}E_{\alpha,\alpha}(-\mu_{j}t^{\alpha})|\,dt&=\int_{0}^{\eta}t^{\alpha-1}E_{\alpha,\alpha}(-\mu_{j}t^{\alpha})\,dt=-\frac{1}{\mu_{j}}\int_{0}^{\eta}\frac{d}{dt}E_{\alpha,1}(-\mu_{j}t^{\alpha})\,dt\\\
&=\frac{1}{\mu_{j}}(1-E_{\alpha,1}(-\mu_{j}\eta^{\alpha})),\;\;\eta>0.\end{aligned}$
(2.13)

Since $W^{1,1}(0,T;V^{*})\subset L^{\infty}(0,T;V^{*})$, (2.3), the definition
of $w(t)$, and (2.13) show

$\displaystyle\begin{aligned} &\mathinner{\lVert
A^{1/2}w(t)\rVert}^{2}=\sum_{j=1}^{\infty}\Big{|}\int_{0}^{t}\langle
f(s),\phi_{j}\rangle\mu_{j}^{1/2}(t-s)^{\alpha-1}E_{\alpha,\alpha}(-\mu_{j}(t-s)^{\alpha})\,ds\Big{|}^{2}\\\
&\lesssim\sum_{j=1}^{\infty}\frac{1}{\mu_{j}}\sup_{0\leq\tau\leq T}|\langle
f(\tau),\phi_{j}\rangle|^{2}\Big{|}\int_{0}^{t}\mu_{j}(t-s)^{\alpha-1}E_{\alpha,\alpha}(-\mu_{j}(t-s)^{\alpha})\,ds\Big{|}^{2}\\\
&=\sum_{j=1}^{\infty}\frac{1}{\mu_{j}}\sup_{0\leq\tau\leq T}|\langle
f(\tau),\phi_{j}\rangle|^{2}\big{(}1-E_{\alpha,1}(-\mu_{j}(t-s)^{\alpha})\big{)}^{2}\lesssim\|f\|_{L^{\infty}(0,T;V^{*})}^{2},\end{aligned}$

with the boundedness of $E_{\alpha,1}(\cdot)$ (cf. (2.6)) in the last step.
This shows $w\in L^{1}(0,T;V)$.

Step 3 (Proof of $v\in W^{1,1}(0,T;V^{*})$). The first equality in (2.7) leads
to

$\displaystyle
v^{\prime}(t)=-t^{\alpha-1}\sum_{j=1}^{\infty}\mu_{j}(u_{0},\phi_{j})E_{\alpha,\alpha}(-\mu_{j}t^{\alpha})\phi_{j}(x).$

This with (2.5) and (2.6) show

$\displaystyle\|v^{\prime}(t)\|_{V^{*}}^{2}=\sum_{j=1}^{\infty}\frac{t^{2\alpha-2}}{\mu_{j}}\mu_{j}^{2}|E_{\alpha,\alpha}(-\mu_{j}t^{\alpha})|^{2}|(u_{0},\phi_{j})|^{2}\lesssim
t^{\alpha-2}\sum_{j=1}^{\infty}\frac{\mu_{j}t^{\alpha}}{(1+(\mu_{j}t^{\alpha})^{2})^{2}}|(u_{0},\phi_{j})|^{2}.$

Since $\frac{\mu_{j}t^{\alpha}}{(1+(\mu_{j}t^{\alpha})^{2})}\leq 1$, the last
displayed inequality and (2.4) reveal $\|v^{\prime}(t)\|_{V^{*}}\lesssim
t^{\alpha/2-1}\|u_{0}\|$, and hence $v\in W^{1,1}(0,T;V^{*})$.

Step 4 (Proof of $\displaystyle w(t)=\int_{0}^{t}F(t-s)f(s)\>ds$). Set

$S_{N}(t,s)=\sum_{j=1}^{N}\langle
f(s),\phi_{j}\rangle(t-s)^{\alpha-1}E_{\alpha,\alpha}(-\mu_{j}(t-s)^{\alpha})\phi_{j}(x),\;t>s.$

An application of Lemma 2.5 with $p=-1=q$ yields

$\displaystyle\|S_{N}(t,s)\|_{V^{*}}^{2}=\sum_{j=1}^{N}\frac{1}{\mu_{j}}|\langle
f(s),\phi_{j}\rangle(t-s)^{\alpha-1}E_{\alpha,\alpha}(-\mu_{j}(t-s)^{\alpha})|^{2}\lesssim(t-s)^{2\alpha-2}\|f\|_{L^{\infty}(0,T;V^{*})}^{2},$

which implies
$\|S_{N}(t,s)\|_{V^{*}}\lesssim(t-s)^{\alpha-1}\|f\|_{L^{\infty}(0,T;V^{*})}$
for all $N\in\mathbb{N}$. In view of the Lebesgue-Böchner dominated
convergence theorem, we obtain

$\displaystyle\lim_{N\rightarrow\infty}\int_{0}^{t}S_{N}(t,s)\,ds=\int_{0}^{t}\lim_{N\rightarrow\infty}S_{N}(t,s)\,ds,$

and this concludes the proof.

Step 5 (Proof of $w\in W^{1,1}(0,T;V^{*})$). A differentiation of $w$ with
respect to $t$ yields

$\displaystyle\begin{aligned}
w^{\prime}(t)=\frac{d}{dt}\Big{(}\int_{0}^{t}F(s)f(t-s)\,ds\Big{)}&=\int_{0}^{t}F(s)f^{\prime}(t-s)\,ds+F(t)f(0).\end{aligned}$
(2.14)

That is,

$\displaystyle\begin{aligned}
\|w^{\prime}(t)\|_{V^{*}}&\leq\int_{0}^{t}\|F(s)f^{\prime}(t-s)\|_{V^{*}}\,ds+\|F(t)f(0)\|_{V^{*}}\\\
&\lesssim\int_{0}^{t}s^{\alpha-1}\|f^{\prime}(t-s)\|_{V^{*}}\,ds+t^{\alpha-1}\|f(0)\|_{V^{*}},\end{aligned}$

with a use of Lemma 2.5 for $p=-1=q$ in the last step above. Since $f\in
W^{1,1}(0,T;V^{*})$, from Young’s inequality for the convolution product to
bound the first term followed by an application of $t\leq T$ for both the
terms, we obtain

$\displaystyle\int_{0}^{T}\|w^{\prime}(t)\|_{V^{*}}\,dt\lesssim\|f^{\prime}\|_{L^{1}(0,T;V^{*})}+\|f(0)\|_{V^{*}}.$

This establishes $w\in W^{1,1}(0,T;V^{*})$.

Step 6 (Proof of $\lim_{t\rightarrow 0}\mathinner{\lVert
u(t)-u_{0}\rVert}_{V^{*}}=0$). Since $w\in W^{1,1}(0,T;V^{*})$, $w$ belongs to
$C([0,T];V^{*})$. Moreover, $w(0)=0$. Therefore,

$\displaystyle\lim_{t\rightarrow 0}\|w(t)\|_{V^{*}}=0.$

Next, we show that

$\displaystyle\lim_{t\rightarrow 0}\mathinner{\lVert v(t)-u_{0}\rVert}=0.$
(2.15)

Employ the definition of $v(t)$ and (2.3) for $r=0$ to obtain

$\displaystyle\mathinner{\lVert
v(t)-u_{0}\rVert}^{2}=\sum_{j=1}^{\infty}|(u_{0},\phi_{j})|^{2}\big{(}E_{\alpha,1}(-\mu_{j}t^{\alpha})-1\big{)}^{2}.$

Also, for each $j\in\mathbb{N}$, $\lim_{t\rightarrow
0}\big{(}E_{\alpha,1}(-\mu_{j}t^{\alpha})-1\big{)}=0$. Moreover, (2.6) applied
to the first inequality below shows, for $0\leq t\leq T$,

$\displaystyle\sum_{j=1}^{\infty}|(u_{0},\phi_{j})|^{2}\big{(}E_{\alpha,1}(-\mu_{j}t^{\alpha})-1\big{)}^{2}\lesssim\sum_{j=1}^{\infty}\Big{(}1+\frac{1}{(1+\mu_{j}t^{\alpha})^{2}}\Big{)}|(u_{0},\phi_{j})|^{2}\leq
2\|u_{0}\|^{2}<\infty.$ (2.16)

Using (2.6), we have

$\displaystyle|(u_{0},\phi_{j})|^{2}\big{(}E_{\alpha,1}(-\mu_{j}t^{\alpha})-1\big{)}^{2}\lesssim|(u_{0},\phi_{j})|^{2}\text{
for all }j\in\mathbb{N}.$

Since $u_{0}\in L^{2}(\Omega)$, the series on the left-hand side of (2.16) is
normally convergent, and so (2.15) is proved.

Since $u=v+w$, a combination of the limits obtained for $v$ and $w$ leads to

$\displaystyle\lim_{t\rightarrow
0}\|u(t)-u_{0}\|_{V^{*}}\leq\lim_{t\rightarrow
0}\|v(t)-u_{0}\|_{V^{*}}+\lim_{t\rightarrow 0}\|w(t)\|_{V^{*}}=0.$

Finally, utilizing (2.9) and (2.11), we can show (2.8). Thus, we have shown
that $u$ defined by (2.9) satisfies all the conditions of a weak solution.
$\Box$

###  2.3 Regularity results

The following regularity results for the solutions of both homogeneous and
nonhomogeneous problems are employed in the error analysis. The proofs of
these results are based on the solution representation in (2.9), the
properties of Mittag-Leffler functions in (2.6) and (2.7), and Lemma 2.5 .

######  Theorem 2.6 (Regularity for homogeneous case).

Let $u$ be the solution of (1.1) with $f=0$. Then, for $t\in(0,T]$, the
following results hold true:

  * (i)

(Nonsmooth initial data) For $u_{0}\in L^{2}(\Omega)$, we have

$\displaystyle\mathinner{\lVert u^{(i)}(t)\rVert}_{D(A^{p})}$
$\displaystyle\lesssim t^{-(i+\alpha p)}\mathinner{\lVert u_{0}\rVert}\text{
for }i\in\\{0,1\\},0\leq p\leq 1,$ (2.17) $\displaystyle\mathinner{\lVert
u^{\prime\prime}(t)\rVert}$ $\displaystyle\lesssim t^{-2}\mathinner{\lVert
u_{0}\rVert},$ (2.18) $\displaystyle\|\partial_{t}^{\alpha}u(t)\|$
$\displaystyle\lesssim t^{-\alpha}\|u_{0}\|.$ (2.19)

In addition, $\|u\|_{C([0,T];L^{2}(\Omega))}\lesssim\|u_{0}\|$.

  * (ii)

(Smooth initial data) For $u_{0}\in D(A)$,

$\displaystyle\mathinner{\lVert u(t)\rVert}_{D(A)}$
$\displaystyle\lesssim\mathinner{\lVert u_{0}\rVert}_{D(A)}$ (2.20)
$\displaystyle\mathinner{\lVert u^{\prime}(t)\rVert}_{D(A^{p})}$
$\displaystyle\lesssim t^{-(1+\alpha(p-1))}\mathinner{\lVert
u_{0}\rVert}_{D(A)}\text{ for }0\leq p\leq 1,$ (2.21)
$\displaystyle\mathinner{\lVert u^{\prime\prime}(t)\rVert}$
$\displaystyle\lesssim t^{\alpha-2}\mathinner{\lVert u_{0}\rVert}_{D(A)},$
(2.22) $\displaystyle\|\partial_{t}^{\alpha}u(t)\|$
$\displaystyle\lesssim\|u_{0}\|_{D(A)}.$ (2.23)

Furthermore, $\|u\|_{C([0,T];D(A))}\lesssim\|u_{0}\|_{D(A)}$.

Proof of case (i). In view of (2.3), (2.9) for $f=0$, and (2.6) it is easy to
see that for $p=0,1$,

$\displaystyle\mathinner{\lVert
A^{p}u(t)\rVert}^{2}=\sum_{j=1}^{\infty}|\mu_{j}^{p}(u_{0},\phi_{j})E_{\alpha,1}(-\mu_{j}t^{\alpha})|^{2}\lesssim\sum_{j=1}^{\infty}|(u_{0},\phi_{j})|^{2}\frac{\mu_{j}^{2p}}{(1+\mu_{j}t^{\alpha})^{2}}\leq
t^{-2p\alpha}\mathinner{\lVert u_{0}\rVert}^{2}.$

Thus for $i=0$, (2.17) follows for $0\leq p\leq 1$ by means of interpolation.

The first identity in (2.7) with $k=1$ and (2.9) with $f=0$ lead to

$\displaystyle
u^{\prime}(t)=-t^{\alpha-1}\sum_{j=1}^{\infty}\mu_{j}(u_{0},\phi_{j})E_{\alpha,\alpha}(-\mu_{j}t^{\alpha})\phi_{j}(x).$

In the next two bounds we utilize
$\frac{(\mu_{j}t^{\alpha})^{i-l}}{(1+(\mu_{j}t^{\alpha})^{2})^{2}}\leq 1$ for
$0\leq i-l\leq 4$ and $j\in\mathbb{N}$. The estimate in (2.6) and the last
displayed equality lead to

$\displaystyle\mathinner{\lVert u^{\prime}(t)\rVert}^{2}$
$\displaystyle\lesssim
t^{-2}\sum_{j=1}^{\infty}|(u_{0},\phi_{j})|^{2}\frac{(\mu_{j}t^{\alpha})^{2}}{(1+(\mu_{j}t^{\alpha})^{2})^{2}}\leq
t^{-2}\sum_{j=1}^{\infty}(u_{0},\phi_{j})^{2}=t^{-2}\mathinner{\lVert
u_{0}\rVert}^{2}.$

Analogous arguments for the case $i=p=1$ in (2.17) reveal

$\displaystyle\mathinner{\lVert Au^{\prime}(t)\rVert}^{2}$
$\displaystyle\lesssim
t^{2\alpha-2}\sum_{j=1}^{\infty}|(u_{0},\phi_{j})|^{2}\frac{\mu_{j}^{4}}{(1+(\mu_{j}t^{\alpha})^{2})^{2}}=t^{-2(1+\alpha)}\sum_{j=1}^{\infty}|(u_{0},\phi_{j})|^{2}\frac{(\mu_{j}t^{\alpha})^{4}}{(1+(\mu_{j}t^{\alpha})^{2})^{2}}$
$\displaystyle\leq
t^{-2(1+\alpha)}\sum_{j=1}^{\infty}|(u_{0},\phi_{j})|^{2}=t^{-2(1+\alpha)}\mathinner{\lVert
u_{0}\rVert}^{2}.$

The above displayed estimates show $\|u^{\prime}(t)\|\lesssim t^{-1}\|u_{0}\|$
and $\mathinner{\lVert Au^{\prime}(t)\rVert}\lesssim
t^{-(1+\alpha)}\mathinner{\lVert u_{0}\rVert}$. An interpolation argument
completes the proof of (2.17).

The first identity in (2.7) with $k=2$ and (2.9) with $f=0$ show

$\displaystyle
u^{\prime\prime}(t)=-t^{\alpha-2}\sum_{j=1}^{\infty}\mu_{j}(u_{0},\phi_{j})E_{\alpha,\alpha-1}(-\mu_{j}t^{\alpha})\phi_{j}(x).$

Now, proceed as in the case for $i=1,p=0$ in (2.17) to obtain (2.18). To prove
(2.19), use (1.1) and (2.17) with $i=0,p=1$ to conclude

$\displaystyle\|\partial_{t}^{\alpha}u(t)\|\leq\|Au(t)\|\lesssim
t^{-\alpha}\|u_{0}\|.$

Proceeding as in Step 6 in the proof of Theorem 2.4, we can show that $u\in
C([0,T];L^{2}(\Omega))$, and this completes the proof of $(i)$. ∎

Proof of (ii). Analogous arguments as in the proof of $(i)$, but now with
$u_{0}\in D(A)$, lead to

$\displaystyle\|u(t)\|^{2}_{D(A)}=\sum_{j=1}^{\infty}|\mu_{j}(u_{0},\phi_{j})E_{\alpha,1}(-\mu_{j}t^{\alpha})|^{2}\lesssim\sum_{j=1}^{\infty}\mu_{j}^{2}|(u_{0},\phi_{j})|^{2}\frac{1}{(1+\mu_{j}t^{\alpha})^{2}}\leq\mathinner{\lVert
u_{0}\rVert}_{D(A)}^{2},$ $\displaystyle\mathinner{\lVert
u^{\prime}(t)\rVert}^{2}\lesssim
t^{2\alpha-2}\sum_{j=1}^{\infty}\mu_{j}^{2}|(u_{0},\phi_{j})|^{2}\frac{1}{(1+(\mu_{j}t^{\alpha})^{2})^{2}}\leq
t^{2\alpha-2}\mathinner{\lVert u_{0}\rVert}_{D(A)}^{2},$ and
$\displaystyle\mathinner{\lVert Au^{\prime}(t)\rVert}^{2}\lesssim
t^{-2}\sum_{j=1}^{\infty}\mu_{j}^{2}|(u_{0},\phi_{j})|^{2}\frac{(\mu_{j}t^{\alpha})^{2}}{(1+(\mu_{j}t^{\alpha})^{2})^{2}}\leq
t^{-2}\mathinner{\lVert u_{0}\rVert}_{D(A)}^{2}.$

The last three displayed inequalities and an interpolation conclude the proof
of (2.20)-(2.21). The proof of (2.22) is similar to that of (2.18) and is
skipped. The proof of (2.23) is straightforward from (1.1) and (2.20).
Following the Step 6 of Theorem 2.4, we can show that $u\in C([0,T];D(A))$,
and this completes the proof of $(ii)$. $\Box$

######  Remark 2.7.

In view of the estimates in Theorem 2.6 $(i)$ and $(ii)$, by means of
interpolation we obtain for time $t\in(0,T]$,

$\displaystyle\|u(t)\|_{D(A^{p})}$ $\displaystyle\lesssim
t^{-\alpha(p-q)}\|u_{0}\|_{D(A^{q})}\text{ for }0\leq q\leq p\leq 1,$
$\displaystyle\|u^{\prime}(t)\|_{D(A^{p})}$ $\displaystyle\lesssim
t^{-(1+\alpha(p-q))}\|u_{0}\|_{D(A^{q})}\text{ for }0\leq p,q\leq 1,$
$\displaystyle\|u^{\prime\prime}(t)\|$ $\displaystyle\lesssim
t^{q\alpha-2}\|u_{0}\|_{D(A^{q})}\text{ for }0\leq q\leq 1,\text{ and},$
$\displaystyle\|\partial_{t}^{\alpha}u(t)\|$ $\displaystyle\lesssim
t^{-(1-q)\alpha}\|u_{0}\|_{D(A^{q})}\text{ for }0\leq q\leq 1.$

######  Theorem 2.8 (Regularity for nonhomogeneous case with zero initial
data).

Let $u$ be the solution of (1.1) with $u_{0}=0$ and $f\in
W^{1,\infty}([0,T];D(A^{q/2}))\cap W^{2,1}(0,T;L^{2}(\Omega))$, $q\in[-1,1]$.
Then, for all $\epsilon\in(0,1)$, we have

$\displaystyle\mathinner{\lVert
u(t)\rVert}_{D(A^{q/2+1-\epsilon/2})}\lesssim\epsilon^{-1}t^{\alpha\epsilon/2}\mathinner{\lVert
f\rVert}_{L^{\infty}(0,t;D(A^{q/2}))},\;\;t\in[0,T],$ (2.24)
$\displaystyle\mathinner{\lVert
u^{\prime}(t)\rVert}_{D(A^{q/2+1-\epsilon/2})}\lesssim\epsilon^{-1}t^{\alpha\epsilon/2}\mathinner{\lVert
f^{\prime}\rVert}_{L^{\infty}(0,t;D(A^{q/2}))}+t^{-1+\alpha\epsilon/2}\mathinner{\lVert
f(0)\rVert}_{D(A^{q/2})},\;\;t\in(0,T],$ (2.25)
$\displaystyle\mathinner{\lVert
u^{\prime\prime}(t)\rVert}\lesssim\mathcal{I}^{\alpha}\mathinner{\lVert
f^{\prime\prime}(\cdot)\rVert}(t)+t^{\alpha-1}\mathinner{\lVert
f^{\prime}(0)\rVert}+t^{\alpha-2}\mathinner{\lVert f(0)\rVert},\;\;t\in(0,T],$
(2.26)
$\displaystyle\|\partial_{t}^{\alpha}u(t)\|\lesssim\epsilon^{-1}t^{\alpha\epsilon/2}\mathinner{\lVert
f\rVert}_{L^{\infty}(0,t;D(A^{\epsilon/2}))}+\|f(t)\|,\;\;t\in[0,T].$ (2.27)

Proof of (2.24). The representation (2.9) for $u_{0}=0$ and Lemma 2.5
establish

$\displaystyle\mathinner{\lVert
u(t)\rVert}_{D(A^{q/2+1-\epsilon/2})}\leq\int_{0}^{t}\mathinner{\\!\left\lVert
F(t-s)f(s)\right\rVert}_{D(A^{q/2+1-\epsilon/2})}\,ds$
$\displaystyle\qquad\qquad\lesssim\int_{0}^{t}(t-s)^{-1+\alpha\epsilon/2}\mathinner{\lVert
f(s)\rVert}_{D(A^{q/2})}\,ds\lesssim\epsilon^{-1}t^{\alpha\epsilon/2}\mathinner{\lVert
f\rVert}_{L^{\infty}(0,t;D(A^{q/2}))}.$

This concludes the proof of (2.24).

Proof of (2.25). Apply Lemma 2.5 with $p=q+2-\epsilon$, $-1\leq q\leq 1$ in
(2.14) to yield, for positive time $t$,

$\displaystyle\begin{aligned} \mathinner{\lVert
u^{\prime}(t)\rVert}_{D(A^{q/2+1-\epsilon/2})}&\leq\int_{0}^{t}\mathinner{\lVert
F(t-s)f^{\prime}(s)\rVert}_{D(A^{q/2+1-\epsilon/2})}\,ds+\mathinner{\lVert
F(t)f(0)\rVert}_{D(A^{q/2+1-\epsilon/2})}\\\
&\lesssim\int_{0}^{t}(t-s)^{-1+\alpha\epsilon/2}\mathinner{\lVert
f^{\prime}(s)\rVert}_{D(A^{q/2})}\,ds+t^{-1+\alpha\epsilon/2}\mathinner{\lVert
f(0)\rVert}_{D(A^{q/2})}\\\
&\lesssim\epsilon^{-1}t^{\alpha\epsilon/2}\mathinner{\lVert
f^{\prime}\rVert}_{L^{\infty}(0,t;D(A^{q/2}))}+t^{-1+\alpha\epsilon/2}\mathinner{\lVert
f(0)\rVert}_{D(A^{q/2})}.\end{aligned}$

This concludes the proof.

Proof of (2.26). Differentiate (2.14) with respect to time and use the first
identity in (2.7) with $k=2$ to obtain

$\displaystyle\begin{aligned}
u^{\prime\prime}(t)&=\int_{0}^{t}F(t-s)f^{\prime\prime}(s)\,ds+F(t)f^{\prime}(0)+F^{\prime}(t)f(0)\\\
&=\int_{0}^{t}F(t-s)f^{\prime\prime}(s)\,ds+F(t)f^{\prime}(0)+t^{\alpha-2}\sum_{j=1}^{\infty}(f(0),\phi_{j})E_{\alpha,\alpha-1}(-\mu_{j}t^{\alpha})\phi_{j}(x).\end{aligned}$
(2.28)

An application of Lemma 2.5 with $p=q=0$ shows

$\displaystyle\begin{aligned}
\mathinner{\lVert\int_{0}^{t}F(t-s)f^{\prime\prime}(s)\,ds\rVert}&\lesssim\int_{0}^{t}(t-s)^{\alpha-1}\mathinner{\lVert
f^{\prime\prime}(s)\rVert}\,ds\lesssim\mathcal{I}^{\alpha}\mathinner{\lVert
f^{\prime\prime}(\cdot)\rVert}(t),\\\ \mathinner{\lVert
F(t)f^{\prime}(0)\rVert}&\lesssim t^{\alpha-1}\mathinner{\lVert
f^{\prime}(0)\rVert}.\end{aligned}$

To bound the last term in (2.28), proceed as the proof of the case $i=1,p=0$
in (2.17) to obtain

$\displaystyle\mathinner{\lVert F^{\prime}(t)f(0)\rVert}\lesssim
t^{\alpha-2}\mathinner{\lVert f(0)\rVert}.$

A combination of the last three estimates in (2.28) completes the proof of
(2.26).

Proof of (2.27). The proof is immediate from (1.1) and (2.24) with
$q=\epsilon$. $\Box$

######  Remark 2.9.

The combination of (2.21) with $\beta\in[1/2,1)$ and (2.25) with
$q\in[2\beta-2+\epsilon,1]$ results in the following estimate for the solution
$u$ of (1.1): For $t\in(0,T]$,

$\displaystyle\|u^{\prime}(t)\|_{D(A^{\beta})}\lesssim
t^{-1+\alpha(1-\beta)}\|u_{0}\|_{D(A)}+\epsilon^{-1}t^{\alpha\epsilon/2}\mathinner{\lVert
f^{\prime}\rVert}_{L^{\infty}(0,t;D(A^{q/2}))}+t^{-1+\alpha\epsilon/2}\mathinner{\lVert
f(0)\rVert}_{D(A^{q/2})}.$

##  3 Semidiscrete scheme

In this section, we describe the lowest-order finite element discretization
schemes for the spatial variable in Subsection 3.1. The Ritz projection
operator and its approximation properties are stated in Subsection 3.2.

###  3.1 Lowest-order finite element discretizations

Let $\mathcal{T}$ denote a shape regular triangulation of the polygonal
Lipschitz domain into compact triangles. Associate its piecewise constant
mesh-size $h_{\mathcal{T}}\in P_{0}(\mathcal{T})$ with
$h_{K}\mathrel{\mathop{\mathchar
58\relax}}=h_{\mathcal{T}}|_{K}\mathrel{\mathop{\mathchar 58\relax}}={\rm
diam}(K)\approx|K|^{1/2}$ in any triangle $K\in\mathcal{T}$ of area $|K|$ and
its maximal mesh-size $h\mathrel{\mathop{\mathchar 58\relax}}={\rm
max}\;h_{\mathcal{T}}$. Let ${\mathcal{V}}$ (resp. ${\mathcal{V}}(\Omega)$ or
${\mathcal{V}}(\partial\Omega)$) denote the set of all (resp. interior or
boundary) vertices in $\mathcal{T}$. Let ${\mathcal{E}}$ (resp.
${\mathcal{E}}(\Omega)$ or ${\mathcal{E}}(\partial\Omega)$) denote the set of
all (resp. interior or boundary) edges. The length of an edge $e$ is denoted
by $h_{e}$. Let the Hilbert space
$H^{m}(\mathcal{T})\equiv\prod_{K\in\mathcal{T}}H^{m}(K)$. The edge-patch
$\omega(e)\mathrel{\mathop{\mathchar 58\relax}}=\text{\rm int}(K_{+}\cup
K_{-})$ of the interior edge $e=\partial K_{+}\cap\partial
K_{-}\in\operatorname{\mathcal{E}}(\Omega)$ is the interior of union
$K_{+}\cup K_{-}$ of the neighboring triangles $K_{+}$ and $K_{-}$; the jump
and average of $\varphi$ are defined by
$\big{[}\\!\\!\big{[}\varphi\big{]}\\!\\!\big{]}\mathrel{\mathop{\mathchar
58\relax}}=\varphi|_{K_{+}}-\varphi|_{K_{-}}$ and
$\big{\\{}\\!\\!\\!\big{\\{}\varphi\big{\\}}\\!\\!\\!\big{\\}}\mathrel{\mathop{\mathchar
58\relax}}=\frac{1}{2}\left(\varphi|_{K_{+}}+\varphi|_{K_{-}}\right)$ across
the interior edge $e$ of the adjacent triangles $K_{+}$ and
$K_{-}\in\mathcal{T}$ in an order such that the unit normal vector
$\nu_{K_{+}}|_{e}=\nu_{e}=-\nu_{K_{-}}|_{e}$ along the edge $e$ has a fixed
orientation and points outside $K_{+}$ and inside $K_{-}$; $\nu_{K}$ is the
outward unit normal of $K$ along $\partial K$. Further for
$e\in\mathcal{E}(\partial\Omega)$, define $\omega(e)\mathrel{\mathop{\mathchar
58\relax}}=\text{int}(K)$, and the jump and average by
$\big{[}\\!\\!\big{[}\varphi\big{]}\\!\\!\big{]}\mathrel{\mathop{\mathchar
58\relax}}=\varphi|_{e}$ and
$\big{\\{}\\!\\!\\!\big{\\{}\varphi\big{\\}}\\!\\!\\!\big{\\}}\mathrel{\mathop{\mathchar
58\relax}}=\varphi|_{e}$.

For functions in $H^{2}_{0}(\Omega)$, the notation
$|\\!|\\!|\cdot|\\!|\\!|\mathrel{\mathop{\mathchar
58\relax}}=|\cdot|_{H^{2}(\Omega)}$ stands for the energy norm. The notation
$|\\!|\\!|\cdot|\\!|\\!|_{\text{pw}}\mathrel{\mathop{\mathchar
58\relax}}=|\cdot|_{H^{2}(\mathcal{T})}\mathrel{\mathop{\mathchar
58\relax}}=\|D^{2}_{\text{pw}}\cdot\|$ refers to the piecewise energy norm
with the piecewise Hessian $D_{\text{pw}}^{2}$. Define the piecewise
polynomials space $P_{r}(\mathcal{T})$ of degree $r\in\mathbb{N}$ by
$P_{r}(\mathcal{T})=\\{v\in L^{2}(\Omega)\mathrel{\mathop{\mathchar
58\relax}}v|_{K}\in P_{r}(K)\text{ for all }K\in\mathcal{T}\\}$.

The nonconforming Morley finite element space $\text{M}(\mathcal{T})$
[ciarlet78] is defined as

$\displaystyle\text{M}(\mathcal{T})\mathrel{\mathop{\mathchar 58\relax}}=\\{$
$\displaystyle\chi_{h}\in P_{2}(\mathcal{T})\mathrel{\mathop{\mathchar
58\relax}}\chi_{h}\text{ is continuous at the vertices and its normal
derivative }{\partial\chi_{h}}/{\partial\nu}\text{ is }$ $\displaystyle\text{
continuous at the midpoints of interior edges},\chi_{h}\text{ vanishes at the
vertices on }\partial\Omega$ $\displaystyle\text{ and its normal derivative
}{\partial\chi_{h}}/{\partial\nu}\text{ vanishes at the midpoints of boundary
edges}\\}.$

On a finite-dimensional space $V_{h}\subset H^{2}(\mathcal{T})$, define a
mesh-dependent broken norm [cgn2015] by

$\displaystyle\begin{aligned}
\mathinner{\lVert\chi_{h}\rVert}_{h}^{2}&=\sum_{K\in\mathcal{T}}|\chi_{h}|_{H^{2}(K)}^{2}+\sum_{e\in\mathcal{E}}\sum_{z\in\mathcal{V}(e)}h_{e}^{-2}\big{|}\big{[}\\!\\!\big{[}\chi_{h}\big{]}\\!\\!\big{]}(z)\big{|}^{2}+\sum_{e\in\mathcal{E}}\Big{|}\fint_{e}\Big{[}\\!\\!\Big{[}\frac{\partial\chi_{h}}{\partial\nu}\Big{]}\\!\\!\Big{]}\,ds\Big{|}^{2},\end{aligned}$
(3.1)

where $\fint_{e}$ denotes the integral mean over the edge $e$. In particular,
for $\chi_{h}\in\mathrm{M}(\mathcal{T})$,
$\|\chi_{h}\|_{h}=|\\!|\\!|\chi_{h}|\\!|\\!|_{\rm pw}$ as the jump terms in
(3.1) vanish. Further, the discrete bilinear form
$a_{h}\mathrel{\mathop{\mathchar
58\relax}}(V_{h}+\text{M}(\mathcal{T}))\times(V_{h}+\text{M}(\mathcal{T}))\rightarrow\mathbb{R}$
in all the examples in this paper has the form

$a_{h}(\cdot,\cdot)=a_{\mathrm{pw}}(\cdot,\cdot)+b_{h}(\cdot,\cdot)+c_{h}(\cdot,\cdot),$

and satisfies (H) below.

(H)

$a_{h}(\cdot,\cdot)$ is symmetric, positive-definite, and continuous on
$V_{h}$ with respect to the discrete norm $\mathinner{\lVert\cdot\rVert}_{h}$,
i.e., $\exists$ constants $\beta_{1},\beta_{2}>0$ independent of $h$ such
that, for all $w_{h},\chi_{h}\in V_{h}$,
$a_{h}(w_{h},\chi_{h})=a_{h}(\chi_{h},w_{h})$ and

$\displaystyle
a_{h}(\chi_{h},\chi_{h})\geq\beta_{1}\mathinner{\lVert\chi_{h}\rVert}_{h}^{2}\text{
and }a_{h}(w_{h},\chi_{h})\leq\beta_{2}\mathinner{\lVert
w_{h}\rVert}_{h}\mathinner{\lVert\chi_{h}\rVert}_{h}.$ (3.2)

The semidiscrete problem that corresponds to the weak formulation in
Definition 2.1 seeks $u_{h}\in W^{1,1}(0,T;V_{h})$ such that

$\displaystyle\partial_{t}^{\alpha}(u_{h}(t),\chi_{h})+a_{h}(u_{h}(t),\chi_{h})$
$\displaystyle=(f(t),\chi_{h})\text{ for all }\chi_{h}\in V_{h},0<t\leq T,$
(3.3) $\displaystyle u_{h}(0)$ $\displaystyle=P_{h}u_{0}\in V_{h}.$

Here $P_{h}\mathrel{\mathop{\mathchar 58\relax}}L^{2}(\Omega)\rightarrow
V_{h}$ denotes the $L^{2}$-projection defined by
$(P_{h}v,\chi_{h})=(v,\chi_{h})$ for all $\chi_{h}\in V_{h}$.

Now we present Examples 3.1-3.3 below for which the discrete bilinear form
$a_{h}(\cdot,\cdot)$ satisfies (H) [cn2022, Section 5].

######  Example 3.1 (Morley).

In this case, $V_{h}\mathrel{\mathop{\mathchar 58\relax}}=\rm M(\mathcal{T})$,
and for all $w_{h},\chi_{h}\in\rm M(\mathcal{T})$,

$\displaystyle\begin{aligned} a_{h}(w_{h},\chi_{h})\mathrel{\mathop{\mathchar
58\relax}}=a_{\rm pw}(w_{h},\chi_{h})\mathrel{\mathop{\mathchar
58\relax}}=\int_{\Omega}D_{\mathrm{pw}}^{2}w_{h}\mathrel{\mathop{\mathchar
58\relax}}D_{\mathrm{pw}}^{2}\chi_{h}\,dx.\end{aligned}$

The discrete Morley norm on $V_{h}$ is defined by
$|\\!|\\!|\cdot|\\!|\\!|_{h}=|\\!|\\!|\cdot|\\!|\\!|_{\mathrm{pw}}=a_{\mathrm{pw}}(\cdot,\cdot)^{1/2}$.
Notice that this norm is equivalent on $V_{h}$ to that introduced in (3.1).

######  Example 3.2 (dG).

Choose $V_{h}\mathrel{\mathop{\mathchar 58\relax}}=P_{2}(\mathcal{T})$, and
for all $w_{h},\chi_{h}\in V_{h}$, let

$\displaystyle\begin{aligned} a_{h}(w_{h},\chi_{h})&\mathrel{\mathop{\mathchar
58\relax}}=a_{\rm pw}(w_{h},\chi_{h})+b_{h}(w_{h},\chi_{h})+c_{\rm
dG}(w_{h},\chi_{h}),\\\ b_{h}(w_{h},\chi_{h})&\mathrel{\mathop{\mathchar
58\relax}}=-\mathcal{J}(w_{h},\chi_{h})-\mathcal{J}(\chi_{h},w_{h}),\mathcal{J}(w_{h},\chi_{h})\mathrel{\mathop{\mathchar
58\relax}}=\sum_{e\in\mathcal{E}}\int_{e}\big{[}\\!\\!\big{[}\nabla
w_{h}\big{]}\\!\\!\big{]}\cdot\big{\\{}\\!\\!\\!\big{\\{}D_{\mathrm{pw}}^{2}\chi_{h}\big{\\}}\\!\\!\\!\big{\\}}\nu\,ds,\\\
c_{\rm dG}(w_{h},\chi_{h})&\mathrel{\mathop{\mathchar
58\relax}}=\sum_{e\in\mathcal{E}}\Big{(}\frac{\sigma_{\rm
dG}^{1}}{h_{e}^{3}}\int_{e}\big{[}\\!\\!\big{[}w_{h}\big{]}\\!\\!\big{]}\big{[}\\!\\!\big{[}\chi_{h}\big{]}\\!\\!\big{]}\,ds+\frac{\sigma_{\rm
dG}^{2}}{h_{e}}\int_{e}\Big{[}\\!\\!\Big{[}\frac{\partial
w_{h}}{\partial\nu}\Big{]}\\!\\!\Big{]}\Big{[}\\!\\!\Big{[}\frac{\partial\chi_{h}}{\partial\nu}\Big{]}\\!\\!\Big{]}\,ds\Big{)},\end{aligned}$

where $\sigma_{\rm dG}^{1}$, $\sigma_{\rm dG}^{2}>0$ are the penalty
parameters and $a_{\rm pw}(\cdot,\cdot)$ is as defined in Example 3.1. The dG
norm $\|\cdot\|_{\rm dG}$ on $P_{2}(\mathcal{T})$ is defined by
$\|w_{h}\|_{\rm dG}=\Big{(}|\\!|\\!|w_{h}|\\!|\\!|_{\mathrm{pw}}^{2}+c_{\rm
dG}(w_{h},w_{h})\Big{)}^{1/2}$, $w_{h}\in V_{h}$. As in the previous example,
this norm is equivalent on $V_{h}$ to that introduced in (3.1).

######  Example 3.3 ($C^{0}$IP).

Choose $V_{h}\mathrel{\mathop{\mathchar 58\relax}}=P_{2}(\mathcal{T})\cap
H_{0}^{1}(\Omega)$ and for all $w_{h},\chi_{h}\in V_{h}$, define

$\displaystyle\begin{aligned} a_{h}(w_{h},\chi_{h})&\mathrel{\mathop{\mathchar
58\relax}}=a_{\rm pw}(w_{h},\chi_{h})+b_{h}(w_{h},\chi_{h})+c_{\rm
IP}(w_{h},\chi_{h}),\\\ c_{\rm IP}(w_{h},\chi_{h})&\mathrel{\mathop{\mathchar
58\relax}}=\sum_{e\in\mathcal{E}}\frac{\sigma_{\rm
IP}}{h_{e}}\int_{e}\Big{[}\\!\\!\Big{[}\frac{\partial
w_{h}}{\partial\nu}\Big{]}\\!\\!\Big{]}\Big{[}\\!\\!\Big{[}\frac{\partial\chi_{h}}{\partial\nu}\Big{]}\\!\\!\Big{]}\,ds,\end{aligned}$

where $\sigma_{\rm IP}$ is a positive parameter, $a_{\rm pw}(\cdot,\cdot)$ and
$b_{h}(\cdot,\cdot)$ are as defined in Examples 3.1 and 3.2. The discrete norm
$\|\cdot\|_{\rm IP}$ on the space $V_{h}$ reads $\|w_{h}\|_{\rm
IP}=\Big{(}|\\!|\\!|w_{h}|\\!|\\!|_{\mathrm{pw}}^{2}+c_{\rm
IP}(w_{h},w_{h})\Big{)}^{1/2}$ for $w_{h}\in V_{h}$. In this example as well,
$\|\cdot\|_{\rm IP}$ is equivalent on $V_{h}$ to that introduced in (3.1).

The equivalence of the common norm $\|\cdot\|_{h}$ with the norms defined in
the Examples 3.1-3.3 (see [cn2022]) is helpful in the proofs of the
approximation properties of the Ritz projection in the next section.

###  3.2 Ritz projection and its approximation properties

We introduce the Ritz projection which is the elliptic projection for the
biharmonic problem and state its approximation properties. The Ritz projection
$\mathcal{R}_{h}\mathrel{\mathop{\mathchar 58\relax}}V\rightarrow V_{h}$ is
defined by

$\displaystyle a_{h}(\mathcal{R}_{h}v,\chi_{h})=a(v,Q\chi_{h})\text{ for all
}\chi_{h}\in V_{h}\text{ and }v\in V,$ (3.4)

where $Q\mathrel{\mathop{\mathchar 58\relax}}=JI_{\text{M}}$ is a smoother
defined from $H^{2}(\mathcal{T})$ to $V$, with
$I_{\text{M}}\mathrel{\mathop{\mathchar
58\relax}}H^{2}(\mathcal{T})\rightarrow\text{M}(\mathcal{T})$ and
$J\mathrel{\mathop{\mathchar 58\relax}}\text{M}(\mathcal{T})\rightarrow V$
denoting the extended Morley interpolation operator and the companion
operator, respectively (see Appendix for the details of the definition and
properties of the interpolation and companion operators). The Lax-Milgram
lemma shows that the projection $\mathcal{R}_{h}$ is a well-defined operator
on $V$.

In the semidiscrete error analysis discussed in Section 4, the error
$u(t)-u_{h}(t)$ is split by introducing the Ritz projection $\mathcal{R}_{h}$
as

$u(t)-u_{h}(t)=(u(t)-\mathcal{R}_{h}u(t))+(\mathcal{R}_{h}u(t)-u_{h}(t)).$

Given the approximation properties for the Ritz projection from Lemma 3.4, the
main task will be to establish bounds for $u_{h}(t)-\mathcal{R}_{h}u(t)$ in
the next section.

######  Lemma 3.4 (Approximation properties).

For any $v\in V$ and the Ritz projection
$\mathcal{R}_{h}\mathrel{\mathop{\mathchar 58\relax}}V\rightarrow V_{h}$
defined in (3.4), the approximation properties in the energy and $L^{2}$ norms
stated below hold.

$\displaystyle\begin{aligned} &\|{v-\mathcal{R}_{h}v}\|+h\mathinner{\lVert
v-\mathcal{R}_{h}v\rVert}_{h}\lesssim h^{1+\gamma}\mathinner{\lVert
v\rVert}_{H^{2+\gamma}(\Omega)}\text{ for all
}\gamma\in[0,\gamma_{0}].\end{aligned}$

The proofs are available in [cn2022]; we provide an outline in the appendix
for continuity in reading. In the error analysis, we also need the following
approximation property of the operator $Q$ for piecewise quadratic functions.

######  Lemma 3.5.

[cn2022, Theorem 4.5 (d)] For any $\chi_{h}\in P_{2}(\mathcal{T})$, the
operator $Q=JI_{\mathrm{M}}$ satisfies

$\displaystyle\|{\chi_{h}-Q\chi_{h}}\|_{H^{s}(\mathcal{T})}\leq
C_{1}h^{2-s}\min_{v\in V}\mathinner{\lVert v-\chi_{h}\rVert}_{h}\text{ for any
}0\leq s\leq 2\text{ and a constant }C_{1}>0.$

##  4 Semidiscrete error estimates

In this section, we derive error bounds for the semidiscrete scheme for both
smooth and nonsmooth initial data. To start with, in Subsection 4.1, we state
some important properties of fractional integrals that are relevant in the
context. This is followed by the main results of this section and their proofs
in Subsections 4.2 and 4.3.

###  4.1 Properties of Riemann-Liouville fractional integrals

For all $\alpha\in(0,\infty)$, all $\beta\in(0,\infty)$ satisfying
$\alpha+\beta\geq 1$, all $v\in L^{1}(0,T)$, the operators ${\cal I}^{\alpha}$
and ${\cal I}^{\beta}$ (cf. (1.3)) satisfy

$\displaystyle{\cal I}^{\alpha}{\cal I}^{\beta}v(t)={\cal
I}^{\alpha+\beta}v(t),\mbox{ for almost all }t\in(0,T).$ (4.1)

If $v\in C([0,T])$, then (4.1) is satisfied at all points $t\in[0,T]$ and
$\alpha,\beta>0$ (cf. [samkokilbasmarichev93, p. 34]). Recall that
$\kappa_{\alpha}(t)\mathrel{\mathop{\mathchar
58\relax}}=t^{\alpha-1}/\Gamma(\alpha)$. In the rest of this subsection, we
assume that $0<\alpha<1$, and so $0<1-\alpha<1$. The three identities below
hold for $v_{1}(t)=tv(t)$ and $v_{2}(t)=t^{2}v(t)$ for $t\in[0,T]$ (see
[karaamustaphapani18, Lemma 2] and [mustapha18, Lemma 2.1] for a proof).

$\displaystyle t\mathcal{I}^{\alpha}v(t)$
$\displaystyle=\mathcal{I}^{\alpha}v_{1}(t)+\alpha\mathcal{I}^{\alpha+1}v(t),$
(4.2) $\displaystyle t\mathcal{I}^{\alpha}v^{\prime}(t)$
$\displaystyle=\mathcal{I}^{\alpha}(v_{1})^{\prime}(t)+(\alpha-1)\mathcal{I}^{\alpha}v(t)-t\kappa_{\alpha}(t)v(0),$
(4.3) $\displaystyle t^{2}\mathcal{I}^{\alpha}v^{\prime}(t)$
$\displaystyle=\mathcal{I}^{\alpha}(v_{2})^{\prime}(t)+2(\alpha-1)\mathcal{I}^{\alpha}v_{1}(t)+\alpha(\alpha-1)\mathcal{I}^{\alpha+1}v(t)-t^{2}\kappa_{\alpha}(t)v(0).$
(4.4)

For $\phi,v\in L^{2}(0,T;L^{2}(\Omega))$, since
$\cos(\alpha\pi/2)-(1-\alpha)\geq 0$, the following continuity property with
any positive $\vartheta$ holds for $\mathcal{I}^{1-\alpha}$
[mustaphaschotzau14, Lemma 3.1 (iii)]:

$\displaystyle\int_{0}^{t}(\mathcal{I}^{1-\alpha}\phi,v)\,ds$
$\displaystyle\leq\vartheta\int_{0}^{t}(\mathcal{I}^{1-\alpha}v,v)\,ds+\frac{1}{4\vartheta(1-\alpha)^{2}}\int_{0}^{t}(\mathcal{I}^{1-\alpha}\phi,\phi)\,ds.$
(4.5)

If $v\mathrel{\mathop{\mathchar 58\relax}}[0,T]\rightarrow L^{2}(\Omega)$ is a
piecewise continuous function in time, then $\mathcal{I}^{\alpha}$ satisfies
[mustaphaschotzau14, Lemma 3.1 (ii)]

$\displaystyle\int_{0}^{T}(\mathcal{I}^{\alpha}v,v)\,dt$
$\displaystyle\geq\cos(\alpha\pi/2)\int_{0}^{T}\|\mathcal{I}^{\alpha/2}v\|^{2}\,dt\geq
0.$ (4.6)

In addition, if $v^{\prime}\mathrel{\mathop{\mathchar
58\relax}}[0,T]\rightarrow L^{2}(\Omega)$ is a piecewise continuous function
in time, then for $t\in(0,T]$, it follows from [lemcleanmustapha16, Lemma 2.1]
that

$\displaystyle\mathinner{\lVert v(t)-v(0)\rVert}^{2}\lesssim
t^{\alpha}\int_{0}^{t}\|\mathcal{I}^{(1-\alpha)/2}v^{\prime}\|^{2}\,ds\lesssim
t^{\alpha}\int_{0}^{t}(\mathcal{I}^{1-\alpha}v^{\prime},v^{\prime})\,ds,$
(4.7)

where the last inequality follows from (4.6).

The proofs of this section use the inequality below frequently. For $a,b\geq
0$, it holds that

$(a+b)^{2}\lesssim(a^{2}+b^{2})\lesssim(a+b)^{2}.\;$ (4.8)

###  4.2 Main results

For the semidiscrete error analysis, we split the error $u(t)-u_{h}(t)$ by
introducing the Ritz projection $\mathcal{R}_{h}$ from (3.4) as

$\displaystyle u(t)-u_{h}(t)=\mathrel{\mathop{\mathchar
58\relax}}\rho(t)+\theta(t),$ (4.9)

with

$\rho(t)\mathrel{\mathop{\mathchar 58\relax}}=u(t)-\mathcal{R}_{h}u(t)\text{
and }\theta(t)\mathrel{\mathop{\mathchar
58\relax}}=\mathcal{R}_{h}u(t)-u_{h}(t).$

Let

$\displaystyle\Lambda_{0}(\epsilon,t)\mathrel{\mathop{\mathchar
58\relax}}=\mathinner{\lVert
u_{0}\rVert}_{D(A)}+\epsilon^{-1}t^{\alpha\epsilon/2}\mathinner{\lVert
f\rVert}_{L^{\infty}(0,T;D(A^{\epsilon/2}))},\;\epsilon\in(0,1),\;t\in[0,T].$
(4.10)

Recall that $D(A)\subset V\cap H^{2+\gamma_{0}}(\Omega)$. For the homogeneous
problem, using (2.17) with $i=0$, $p=1$, and for the nonhomogeneous problem,
applying (2.20) and (2.24) with $q=\epsilon$ we obtain

$\displaystyle\|u(t)\|_{H^{2+\gamma_{0}}(\Omega)}$
$\displaystyle\lesssim\begin{cases}t^{-\alpha}\mathinner{\lVert
u_{0}\rVert}\text{ for }\;u_{0}\in L^{2}(\Omega)\text{ and
}f=0,\;t\in(0,T],\\\ \Lambda_{0}(\epsilon,t)\text{ for }u_{0}\in D(A)\text{
and }f\neq 0,\;t\in[0,T].\end{cases}$

This and Lemma 3.4 establish the estimates for $\rho(t)$ as given below.

$\displaystyle\mathinner{\lVert\rho(t)\rVert}+h\mathinner{\lVert\rho(t)\rVert}_{h}$
$\displaystyle\lesssim\begin{cases}h^{1+\gamma_{0}}t^{-\alpha}\mathinner{\lVert
u_{0}\rVert}\text{ for }u_{0}\in L^{2}(\Omega)\text{ and }f=0,\;t\in(0,T],\\\
h^{1+\gamma_{0}}\Lambda_{0}(\epsilon,t)\text{ for }u_{0}\in D(A)\text{ and
}f\neq 0,\;t\in[0,T].\end{cases}$ (4.11)

Hence the main task in the remaining part of this section is to bound
$\theta(t)$ in the $L^{2}(\Omega)$ and energy norms. For $t\in[0,T]$ and an
$\epsilon\in(0,1)$ (determined by the smoothness of $f$), set

$\displaystyle\begin{aligned}
\Lambda_{1}(\epsilon,t)&\mathrel{\mathop{\mathchar
58\relax}}=t^{\frac{3}{2}}\mathinner{\lVert
u_{0}\rVert}_{D(A)}+t^{\frac{3}{2}-\alpha(1-\frac{\epsilon}{2})}\mathinner{\lVert
f(0)\rVert}_{D(A^{\frac{\epsilon}{2}})}+t\|f\|_{L^{2}(0,t;L^{2}(\Omega))}\\\
&\qquad+\epsilon^{-1}t^{\frac{5}{2}-\alpha(1-\frac{\epsilon}{2})}\mathinner{\lVert
f^{\prime}\rVert}_{L^{\infty}(0,T;D(A^{\frac{\epsilon}{2}}))}\\\
\Lambda_{2}(\epsilon,t)&\mathrel{\mathop{\mathchar
58\relax}}=t^{\frac{3}{2}}\mathinner{\lVert
u_{0}\rVert}_{D(A)}+B(\frac{\alpha\epsilon}{2},1-\alpha)t^{\frac{3}{2}-\alpha(1-\frac{\epsilon}{2})}\mathinner{\lVert
f(0)\rVert}_{D(A^{\frac{\epsilon}{2}})}+t\|f\|_{L^{2}(0,t;L^{2}(\Omega))}\\\
&\qquad+\epsilon^{-1}t^{\frac{5}{2}-\alpha(1-\frac{\epsilon}{2})}\mathinner{\lVert
f^{\prime}\rVert}_{L^{\infty}(0,T;D(A^{\frac{\epsilon}{2}}))}\\\
\mathcal{B}_{0}(\epsilon,t)&\mathrel{\mathop{\mathchar
58\relax}}=\mathinner{\lVert
u_{0}\rVert}_{D(A)}+\epsilon^{-1}t^{\frac{\alpha\epsilon}{2}}\|f\|_{W^{1,\infty}(0,T;D(A^{\epsilon/2}))}\\\
\mathcal{B}_{1}(\epsilon,t)&\mathrel{\mathop{\mathchar
58\relax}}=\Lambda_{1}(\epsilon,t)+\Lambda_{2}(\epsilon,t)+t^{2}\|f(t)\|+t\|f\|_{L^{2}(0,t;L^{2}(\Omega))}+t^{2}\|f^{\prime}\|_{L^{2}(0,t;L^{2}(\Omega))}\\\
&\qquad+t^{\frac{3}{2}}(\mathinner{\lVert f(0)\rVert}+t\mathinner{\lVert
f^{\prime}(0)\rVert})+t^{5/2}\|f^{\prime\prime}\|_{L^{1}(0,T;L^{2}(\Omega))},\end{aligned}$
(4.12)

where $B(\cdot,\cdot)$ denotes the standard beta function. Note that in the
above expression and in the sequel, whenever the norm of function $f$ is
dependent only on the space variable, we denote the dependence of $f$ on $t$
as $f(t)$; if the norm is space-time dependent, the arguments in $f$ are
omitted for notational brevity.

######  Theorem 4.1 (Estimates for $\theta(t)$).

Let $u(t)$ and $u_{h}(t)$ solve (1.1) and (3.3), respectively. Let
$\mathcal{R}_{h}u(t)$ denote the Ritz projection of $u(t)$ defined in (3.4).
Then for $\theta(t)=\mathcal{R}_{h}u(t)-u_{h}(t)$, the estimates in
$(i)$-$(ii)$ below hold.

  * (i)

(Nonsmooth initial data) For $u_{0}\in L^{2}(\Omega)$ and $f=0$,

$\displaystyle\mathinner{\lVert\theta(t)\rVert}+t^{\alpha/2}\|\theta(t)\|_{h}$
$\displaystyle\lesssim\big{(}h^{1+\gamma_{0}}t^{-\alpha}+h^{2}t^{-(1+\alpha)/2}\big{)}\mathinner{\lVert
u_{0}\rVert},\;\;t\in(0,T].$

  * (ii)

(Smooth initial data) For $u_{0}\in D(A)$, $f\in
W^{1,\infty}([0,T];D(A^{\epsilon/2}))\cap W^{2,1}(0,T;L^{2}(\Omega))$, and for
all $\epsilon\in(0,1)$,

$\displaystyle\mathinner{\lVert\theta(t)\rVert}+t^{\alpha/2}\|\theta(t)\|_{h}$
$\displaystyle\lesssim
h^{1+\gamma_{0}}{\mathcal{B}_{0}(\epsilon,t)}+h^{2}t^{\alpha/2-2}\mathcal{B}_{1}(\epsilon,t),\;\;t\in(0,T],$

with $\mathcal{B}_{0}(\epsilon,t)$ and $\mathcal{B}_{1}(\epsilon,t)$ defined
in (4.12).

A combination of the estimates for $\rho(t)$ from (4.11) and $\theta(t)$ from
Theorem 4.1, shows the error estimates for the semidiscrete scheme for both
smooth and nonsmooth initial data in Theorem 4.2.

######  Theorem 4.2 (Error estimates).

Let $u(t)$ and $u_{h}(t)$ be solutions of the continuous and semidiscrete
problems in (1.1) and (3.3), respectively.

  * (i)

(Nonsmooth initial data) For $u_{0}\in L^{2}(\Omega)$, source function $f=0$,
$u_{h}(0)=P_{h}u_{0}$, and $t\in(0,T]$, it holds that

$\displaystyle\mathinner{\lVert u(t)-u_{h}(t)\rVert}$
$\displaystyle\lesssim\big{(}h^{1+\gamma_{0}}t^{-\alpha}+h^{2}t^{-(1+\alpha)/2}\big{)}\mathinner{\lVert
u_{0}\rVert},$ $\displaystyle\mathinner{\lVert u(t)-u_{h}(t)\rVert}_{h}$
$\displaystyle\lesssim h^{\gamma_{0}}t^{-\alpha}\mathinner{\lVert
u_{0}\rVert}+\big{(}h^{1+\gamma_{0}}t^{-3\alpha/2}+h^{2}t^{-(\alpha+1/2)}\big{)}\mathinner{\lVert
u_{0}\rVert}.$

  * (ii)

(Smooth initial data) For $u_{0}\in D(A)$, $f\in
W^{1,\infty}([0,T];D(A^{\epsilon/2}))\cap W^{2,1}(0,T;L^{2}(\Omega))$,
$u_{h}(0)=P_{h}u_{0}$, $t\in(0,T]$, and for all $\epsilon\in(0,1)$, it holds
that

$\displaystyle\mathinner{\lVert u(t)-u_{h}(t)\rVert}$ $\displaystyle\lesssim
h^{1+\gamma_{0}}\Lambda_{0}(\epsilon,t)+h^{1+\gamma_{0}}{\mathcal{B}_{0}(\epsilon,t)}+h^{2}t^{\alpha/2-2}\mathcal{B}_{1}(\epsilon,t),$
$\displaystyle\mathinner{\lVert u(t)-u_{h}(t)\rVert}_{h}$
$\displaystyle\lesssim
h^{\gamma_{0}}\Lambda_{0}(\epsilon,t)+h^{1+\gamma_{0}}t^{-\alpha/2}{\mathcal{B}_{0}(\epsilon,t)}+h^{2}t^{-2}\mathcal{B}_{1}(\epsilon,t),$

where $\Lambda_{0}(\epsilon,t)$ and $\mathcal{B}_{0}(\epsilon,t)$,
$\mathcal{B}_{1}(\epsilon,t)$ are defined, respectively, in (4.10) and (4.12).

###### Proof.

The error decomposition (4.9), the triangle inequality, the first estimate in
(4.11), and Theorem 4.1 $(i)$ lead to the proof of $(i)$. The proof of $(ii)$
follows from (4.9), the triangle inequality, the second estimate in (4.11),
and Theorem 4.1 $(ii)$. ∎

###  4.3 Proof of Theorem 4.1

First of all, a key inequality is proved in Lemma 4.3 and bounds for the terms
appearing in this lemma are established in Lemmas 4.5–4.8. A combination of
these results establishes the proof of Theorem 4.1.

The following notations hold throughout this subsection.

$\displaystyle\begin{aligned} \theta_{1}(t)&\mathrel{\mathop{\mathchar
58\relax}}=t\theta(t),\;\theta_{2}(t)\mathrel{\mathop{\mathchar
58\relax}}=t^{2}\theta(t),\;\;\text{and}\\\
\widehat{v}(t)&\mathrel{\mathop{\mathchar
58\relax}}=\mathcal{I}v(t)=\int_{0}^{t}v(s)\,ds,\;\widehat{\widehat{v}}(t)\mathrel{\mathop{\mathchar
58\relax}}=\mathcal{I}^{2}v(t)=\int_{0}^{t}\int_{0}^{s}v(\tau)\,d\tau\,ds.\end{aligned}$

Recall the smoother $Q$ from Subsection 3.2. For each $\chi_{h}\in V_{h}$,
test (2.8) with $Q\chi_{h}\in V$ and subtract (3.3) from (2.8) to obtain

$\displaystyle(\partial_{t}^{\alpha}u(t),Q\chi_{h})-(\partial_{t}^{\alpha}u_{h}(t),\chi_{h})+(a(u(t),Q\chi_{h})-a_{h}(u_{h}(t),\chi_{h}))=(f(t),(Q-I)\chi_{h}).$

Add and subtract $(\partial_{t}^{\alpha}u(t),\chi_{h})$, utilize (3.4) and
$\theta(t)\mathrel{\mathop{\mathchar 58\relax}}=R_{h}u(t)-u_{h}(t)$ on the
left-hand side of the above expression to obtain

$\displaystyle(\partial_{t}^{\alpha}(u(t)-u_{h}(t)),\chi_{h})+(\partial_{t}^{\alpha}u(t),(Q-I)\chi_{h})+a_{h}(\theta(t),Q\chi_{h})=(f(t),(Q-I)\chi_{h}).$

Recall $u(t)-u_{h}(t)=\rho(t)+\theta(t)$ and the notation from (1.2), that
yields
$\partial_{t}^{\alpha}\theta(t)=\mathcal{I}^{1-\alpha}\theta^{\prime}(t)$,
$\partial_{t}^{\alpha}\rho(t)=\mathcal{I}^{1-\alpha}\rho^{\prime}(t)$, and
$\partial_{t}^{\alpha}u(t)=\mathcal{I}^{1-\alpha}u^{\prime}(t)$. For ease of
notation in the subsequent estimates, we denote
$\varphi(t)\mathrel{\mathop{\mathchar
58\relax}}=\mathcal{I}^{1-\alpha}u^{\prime}(t)$. Then for all $\chi_{h}\in
V_{h}$ and $t\in(0,T]$, the above displayed identity leads to the error
equation in $\theta(t)$ as

$\displaystyle(\mathcal{I}^{1-\alpha}\theta^{\prime}(t),\chi_{h})+a_{h}(\theta(t),\chi_{h})=-(\mathcal{I}^{1-\alpha}\rho^{\prime}(t),\chi_{h})-(\varphi(t),(Q-I)\chi_{h})+(f(t),(Q-I)\chi_{h}).$
(4.13)

Recall the notations

$\displaystyle\theta_{i}(t)=t^{i}\theta(t),\;\rho_{i}(t)=t^{i}\rho(t)\text{
for
}i=1,2;\;\varphi_{2}(t)=t^{2}\varphi(t),\;f_{2}(t)=t^{2}f(t),\widehat{\theta}(t)=\int_{0}^{t}\theta(s)\,ds,$
$\displaystyle\widehat{\rho}(t)=\int_{0}^{t}\rho(s)\,ds,\;\rho_{2}^{\prime}(t)=(t^{2}\rho(t))^{\prime},\;\varphi_{2}^{\prime}(t)=(t^{2}\varphi(t))^{\prime},\;f_{2}^{\prime}(t)=(t^{2}f(t))^{\prime},\text{
and }$
$\displaystyle\lambda(t)=-\rho_{2}^{\prime}(t)+2\alpha\rho_{1}(t)+\alpha(1-\alpha)\widehat{\rho}(t).$

######  Lemma 4.3 (Key inequality).

If $\theta$ satisfies (4.13), then for time $t\in(0,T]$, it holds that

$\displaystyle\mathinner{\lVert\theta(t)\rVert}^{2}+t^{\alpha}\mathinner{\lVert\theta(t)\rVert}_{h}^{2}\lesssim
t^{\alpha-4}\int_{0}^{t}\Big{(}(\mathcal{I}^{1-\alpha}\widehat{\theta},\widehat{\theta})\,ds+(\mathcal{I}^{1-\alpha}\theta_{1},\theta_{1})+(\mathcal{I}^{1-\alpha}\lambda,\lambda)\Big{)}\,ds$
$\displaystyle+h^{4}t^{\alpha-4}\Big{(}\|\varphi_{2}(t)\|^{2}+\int_{0}^{t}(\|\varphi_{2}(s)\|^{2}+\|\varphi_{2}^{\prime}(s)\|^{2})\,ds+\|f_{2}(t)\|^{2}+\int_{0}^{t}(\|f_{2}(s)\|^{2}+\|f_{2}^{\prime}(s)\|^{2})\,ds\Big{)}.$

######  Remark 4.4.

Note the dependency of the variables
$\widehat{\theta},\widehat{\theta_{1}},\lambda$, etc. on the time variable is
suppressed in the above and in the sequel for notational ease when there is no
chance of confusion.

###### Proof.

To prove the assertion, we first show that $\theta_{2}\in W^{1,1}(0,T;V_{h})$
and $\varphi_{2}\in W^{1,1}(0,T;L^{2}(\Omega))$. For nonsmooth initial data,
from (2.17), it follows that

$\|u^{\prime}(t)\|_{V}=\|u^{\prime}(t)\|_{D(A^{1/2})}\lesssim
t^{-(1+\alpha/2)}\|u_{0}\|.$

From which, we deduce that

$\|{\mathcal{R}}_{h}u^{\prime}(t)\|_{h}\lesssim t^{-(1+\alpha/2)}\|u_{0}\|,$

and that $t\mapsto t^{2}{\mathcal{R}}_{h}u(t)$ belongs to
$W^{1,1}(0,T;V_{h})$. We already know that $u_{h}$ belongs to
$W^{1,1}(0,T;V_{h})$. Thus $\theta_{2}$ also belongs to $W^{1,1}(0,T;V_{h})$.
With (2.19), we have

$\|\varphi(t)\|\lesssim t^{-\alpha}\|u_{0}\|.$

Thus $\varphi_{2}$ belongs to $W^{1,1}(0,T;L^{2}(\Omega))$. A similar argument
holds for smooth initial data.

Multiply both sides of (4.13) by $t^{2}$, apply (4.4) (twice) to rewrite both
$t^{2}\mathcal{I}^{1-\alpha}\theta^{\prime}$ and
$t^{2}\mathcal{I}^{1-\alpha}\rho^{\prime}$, utilize
$(\rho(0)+\theta(0),\chi_{h})=0$ from the second identity in (3.3), and the
definition of $\lambda$ to obtain Substitute $\chi_{h}=\theta_{2}^{\prime}(t)$
in the last displayed equation, integrate over $(0,t)$, and apply (4.5) for
the first three terms on the right-hand side with the choice of $\vartheta$ as
$\frac{1}{2(2\alpha+\alpha(1-\alpha)+1)}$ to obtain

$\displaystyle\begin{aligned}
&\int_{0}^{t}(\mathcal{I}^{1-\alpha}\theta_{2}^{\prime},\theta_{2}^{\prime})\,ds+\int_{0}^{t}a_{h}(\theta_{2},\theta_{2}^{\prime})\,ds\leq\frac{1}{2}\int_{0}^{t}(\mathcal{I}^{1-\alpha}\theta_{2}^{\prime},\theta_{2}^{\prime})\,ds+C\int_{0}^{t}\Big{(}(\mathcal{I}^{1-\alpha}\widehat{\theta},\widehat{\theta})\\\
&+(\mathcal{I}^{1-\alpha}\theta_{1},\theta_{1})+(\mathcal{I}^{1-\alpha}\lambda,\lambda)\Big{)}\,ds-\int_{0}^{t}(\varphi_{2},(Q-I)\theta_{2}^{\prime})\,ds+\int_{0}^{t}(f_{2},(Q-I)\theta_{2}^{\prime})\,ds.\end{aligned}$

The constant $C$ in the above inequality depends on $\alpha$. The symmetry of
$a_{h}(\cdot,\cdot)$ from (H) shows
$2a_{h}(\theta_{2}(t),\theta_{2}^{\prime}(t))=\frac{d}{dt}a_{h}(\theta_{2}(t),\theta_{2}(t))$.
This with (3.2), (4.7), and $\theta_{2}(0)=0$ on the left-hand side of the
above inequality establishes

$\displaystyle
t^{-\alpha}\mathinner{\lVert\theta_{2}(t)\rVert}^{2}+\beta_{1}\mathinner{\lVert\theta_{2}(t)\rVert}_{h}^{2}$
$\displaystyle\lesssim\int_{0}^{t}\Big{(}(\mathcal{I}^{1-\alpha}\widehat{\theta},\widehat{\theta})+(\mathcal{I}^{1-\alpha}\theta_{1},\theta_{1})+(\mathcal{I}^{1-\alpha}\lambda,\lambda)\Big{)}\,ds$
$\displaystyle\qquad+|\int_{0}^{t}(\varphi_{2},(Q-I)\theta_{2}^{\prime})\,ds|+|\int_{0}^{t}(f_{2},(Q-I)\theta_{2}^{\prime})\,ds|.$
(4.14)

Now, the task is to bound the last two terms on the right-hand side of (4.3).
Since $\theta_{2}$ and $\varphi_{2}$ belong to $W^{1,1}(0,T;L^{2}(\Omega))$,
with an integration by parts, we have

$\displaystyle T_{1}\mathrel{\mathop{\mathchar
58\relax}}=\int_{0}^{t}(\varphi_{2}(s),(Q-I)\theta_{2}^{\prime}(s))\,ds=(\varphi_{2}(t),(Q-I)\theta_{2}(t))-\int_{0}^{t}(\varphi_{2}^{\prime}(s),(Q-I)\theta_{2}(s))\,ds.$

The Cauchy-Schwarz inequality and Lemma 3.5 show

$\displaystyle\begin{aligned} |T_{1}|&\leq
C_{1}h^{2}\big{(}\mathinner{\lVert\varphi_{2}(t)\rVert}\|\theta_{2}(t)\|_{h}+\int_{0}^{t}\mathinner{\lVert\varphi_{2}^{\prime}(s)\rVert}\|\theta_{2}(s)\|_{h}\,ds\big{)}\\\
&\leq\frac{\beta_{1}}{4}\|\theta_{2}(t)\|_{h}^{2}+C_{1}^{2}\beta_{1}^{-1}h^{4}\mathinner{\lVert\varphi_{2}(t)\rVert}^{2}+\frac{1}{2}\int_{0}^{t}\|\theta_{2}(s)\|_{h}^{2}\,ds+\frac{C_{1}^{2}h^{4}}{2}\int_{0}^{t}\mathinner{\lVert\varphi_{2}^{\prime}(s)\rVert}^{2}\,ds,\end{aligned}$

with an application of Young’s inequality in the last step above. A similar
approach is applied to bound the term $T_{2}\mathrel{\mathop{\mathchar
58\relax}}=\int_{0}^{t}(f_{2},(Q-I)\theta_{2}^{\prime})\,ds$ and leads to

$\displaystyle\begin{aligned} |T_{2}|&\leq C_{1}h^{2}\big{(}\mathinner{\lVert
f_{2}(t)\rVert}\|\theta_{2}(t)\|_{h}+\int_{0}^{t}\mathinner{\lVert
f_{2}^{\prime}(s)\rVert}\|\theta_{2}(s)\|_{h}\,ds\big{)}\\\
&\leq\frac{\beta_{1}}{4}\|\theta_{2}(t)\|_{h}^{2}+C_{1}^{2}\beta_{1}^{-1}h^{4}\mathinner{\lVert
f_{2}(t)\rVert}^{2}+\frac{1}{{2}}\int_{0}^{t}\|\theta_{2}(s)\|_{h}^{2}\,ds+\frac{C_{1}^{2}h^{4}}{2}\int_{0}^{t}\mathinner{\lVert
f_{2}^{\prime}(s)\rVert}^{2}\,ds.\end{aligned}$

Substitute these bounds of $|T_{1}|$ and $|T_{2}|$ in (4.3) to obtain

$\displaystyle
t^{-\alpha}\mathinner{\lVert\theta_{2}(t)\rVert}^{2}+\beta_{1}\mathinner{\lVert\theta_{2}(t)\rVert}_{h}^{2}\lesssim\int_{0}^{t}\Big{(}(\mathcal{I}^{1-\alpha}\widehat{\theta},\widehat{\theta})+(\mathcal{I}^{1-\alpha}\theta_{1},\theta_{1})+(\mathcal{I}^{1-\alpha}\lambda,\lambda)\Big{)}\,ds$
$\displaystyle+h^{4}\Big{(}\|\varphi_{2}(t)\|^{2}+\int_{0}^{t}\|\varphi_{2}^{\prime}(s)\|^{2}\,ds+\|f_{2}(t)\|^{2}+\int_{0}^{t}\|f_{2}^{\prime}(s)\|^{2}\,ds\Big{)}+\int_{0}^{t}\|\theta_{2}(s)\|_{h}^{2}\,ds.$

Apply Gronwall’s lemma in Lemma 2.1 with

$\displaystyle\phi(t)=t^{-\alpha}\mathinner{\lVert\theta_{2}(t)\rVert}^{2}+\beta_{1}\mathinner{\lVert\theta_{2}(t)\rVert}_{h}^{2},\;\chi(t)=C,\text{
and }$
$\displaystyle\psi(t)=C\int_{0}^{t}\big{(}(\mathcal{I}^{1-\alpha}\widehat{\theta},\widehat{\theta})+(\mathcal{I}^{1-\alpha}\theta_{1},\theta_{1})+(\mathcal{I}^{1-\alpha}\lambda,\lambda)\big{)}\,ds$
$\displaystyle\qquad\qquad+Ch^{4}\big{(}\|\varphi_{2}(t)\|^{2}+\int_{0}^{t}\|\varphi_{2}^{\prime}(s)\|^{2}\,ds+\|f_{2}(t)\|^{2}+\int_{0}^{t}\|f_{2}^{\prime}(s)\|^{2}\,ds\big{)},$

utilize (4.6),
$\displaystyle\int_{0}^{t}\int_{0}^{s}g(\tau)\,d\tau\,ds\lesssim\int_{0}^{t}g(s)\,ds$
for $\displaystyle\int_{0}^{t}g(s)\,ds>0$ ($g(t)\geq 0$) and $t\in(0,T]$ for
the terms corresponding to a double integral in time, and recall
$\theta_{2}(t)=t^{2}\theta(t)$ to conclude the proof. ∎

Lemmas 4.5-4.8 bound each term on the right-hand side of the estimate in Lemma
4.3.

The notations $\displaystyle\widehat{\upsilon}(t)=\int_{0}^{t}\upsilon(s)\,ds$
for the choices $\upsilon=\theta,\rho,\varphi,f$ and
$\displaystyle\widehat{\widehat{\upsilon}}(t)=\int_{0}^{t}\int_{0}^{s}\upsilon(\tau)\,d\tau\,ds$
for the choices $\upsilon=\theta,\varphi,f$ are used in the next lemma.

######  Lemma 4.5 (Estimate for
$\int_{0}^{t}(\mathcal{I}^{1-\alpha}\widehat{\theta},\widehat{\theta})\,ds$).

For $t\in(0,T]$ and $\theta$ satisfying (4.13), the bounds stated below hold.

$\displaystyle\int_{0}^{t}(\mathcal{I}^{1-\alpha}\widehat{\theta},\widehat{\theta})\,ds+\|\widehat{\widehat{\theta}}(t)\|_{h}^{2}\lesssim
t^{4-\alpha}\Big{(}h^{1+\gamma_{0}}t^{-\alpha}+h^{2}t^{-(\alpha+1)/2}\Big{)}^{2}\mathinner{\lVert
u_{0}\rVert}^{2},$

if $u_{0}\in L^{2}(\Omega)$ and $f=0$, and

$\displaystyle\int_{0}^{t}(\mathcal{I}^{1-\alpha}\widehat{\theta},\widehat{\theta})\,ds+\|\widehat{\widehat{\theta}}(t)\|_{h}^{2}\lesssim
t^{4-\alpha}\Big{(}h^{1+\gamma_{0}}{\Lambda_{0}(\epsilon,t)}+h^{2}t^{\alpha/2-2}{\Lambda_{1}(\epsilon,t)}\Big{)}^{2},$

for all $\epsilon\in(0,1)$, if $u_{0}\in D(A)$ and $f$ is such that
$\Lambda_{0}(\epsilon,t)$ and $\Lambda_{1}(\epsilon,t)$, defined in (4.10) and
(4.12), are finite.

###### Proof.

The proof is split into four steps. The first step derives an intermediate
bound, and Steps 2-4 estimate the terms derived in Step 1.

Step 1 (An intermediate bound). The definition of $\mathcal{I}^{\beta}(\cdot)$
in (1.3) and (4.1) establish

$\mathcal{I}^{2-\alpha}v^{\prime}(t)=\mathcal{I}^{1-\alpha}\mathcal{I}(v^{\prime}(t))=\mathcal{I}^{1-\alpha}(v(t)-v_{0})=\mathcal{I}^{1-\alpha}v(t)-\kappa_{2-\alpha}(t)v(0).$

Integrate (4.13) over $(0,t)$, utilize the above identity twice, and
$(\rho(0)+\theta(0),\chi_{h})=0$ for all $\chi_{h}\in V_{h}$ from the second
identity in (3.3) to obtain

$\displaystyle(\mathcal{I}^{1-\alpha}\theta,\chi_{h})+a_{h}(\widehat{\theta},\chi_{h})$
$\displaystyle=-(\mathcal{I}^{1-\alpha}\rho,\chi_{h})-(\widehat{\varphi},(Q-I)\chi_{h})+(\widehat{f},(Q-I)\chi_{h})\text{
for all }\chi_{h}\in V_{h}.$

Integrate the above identity in time from $0$ to $t$ and choose
$\chi_{h}=\widehat{\theta}(t)$ to establish

$\displaystyle\begin{aligned}
(\mathcal{I}^{1-\alpha}\widehat{\theta},\widehat{\theta})+a_{h}(\widehat{\widehat{\theta}},\widehat{\theta})&=-(\mathcal{I}^{1-\alpha}\widehat{\rho},\widehat{\theta})-(\widehat{\widehat{\varphi}},(Q-I)\widehat{\theta})+(\widehat{\widehat{f}},(Q-I)\widehat{\theta}).\end{aligned}$

Integrate the above identity once again in time from $0$ to $t$, apply (4.5)
with $\vartheta=1/2$ to
$\int_{0}^{t}(\mathcal{I}^{1-\alpha}(-\widehat{\rho}),\widehat{\theta})\,ds$,
utilize
$2a_{h}(\widehat{\widehat{\theta}},\widehat{\theta})=\frac{d}{dt}a_{h}(\widehat{\widehat{\theta}},\widehat{\widehat{\theta}})$,
and (3.2) to obtain

$\displaystyle\begin{aligned}
\int_{0}^{t}(\mathcal{I}^{1-\alpha}\widehat{\theta},\widehat{\theta})\,ds+\beta_{1}\|\widehat{\widehat{\theta}}\|_{h}^{2}&\lesssim\int_{0}^{t}(\mathcal{I}^{1-\alpha}\widehat{\rho},\widehat{\rho})\,ds\\\
&+\Big{|}\int_{0}^{t}(\widehat{\widehat{\varphi}},(Q-I)\widehat{\theta})\,ds\Big{|}+\Big{|}\int_{0}^{t}(\widehat{\widehat{f}},(Q-I)\widehat{\theta})\,ds\Big{|}.\end{aligned}$

The bounds for the terms
$\big{|}\int_{0}^{t}(\widehat{\widehat{\varphi}},(Q-I)\widehat{\theta})\,ds\big{|}$
and
$\big{|}\int_{0}^{t}(\widehat{\widehat{f}},(Q-I)\widehat{\theta})\,ds\big{|}$
can be established in an analogous way following the steps of the derivation
of the bounds for $|T_{1}|$ and $|T_{2}|$ in Lemma 4.3. In this case, we
obtain

$\displaystyle\int_{0}^{t}(\mathcal{I}^{1-\alpha}\widehat{\theta},\widehat{\theta})\,ds+\|\widehat{\widehat{\theta}}(t)\|_{h}^{2}$
$\displaystyle\lesssim\int_{0}^{t}(\mathcal{I}^{1-\alpha}\widehat{\rho},\widehat{\rho})\,ds+h^{4}\Big{(}\|\widehat{\widehat{\varphi}}(t)\|^{2}+\int_{0}^{t}\|\widehat{\varphi}(s)\|^{2}\,ds$
$\displaystyle\quad+\|\widehat{\widehat{f}}(t)\|^{2}+\int_{0}^{t}\|\widehat{f}(s)\|^{2}\,ds\Big{)}+\int_{0}^{t}\|\widehat{\widehat{\theta}}(s)\|_{h}^{2}\,ds.$
(4.15)

Apply Gronwall’s lemma in (4.3) to obtain

$\displaystyle\int_{0}^{t}(\mathcal{I}^{1-\alpha}\widehat{\theta},\widehat{\theta})\,ds+\|\widehat{\widehat{\theta}}(t)\|_{h}^{2}$
$\displaystyle\lesssim\int_{0}^{t}(\mathcal{I}^{1-\alpha}\widehat{\rho},\widehat{\rho})\,ds+h^{4}\Big{(}\|\widehat{\widehat{\varphi}}(t)\|^{2}+\int_{0}^{t}(\|\widehat{\varphi}(s)\|^{2}+{\|\widehat{\widehat{\varphi}}(s)\|^{2})}\,ds$
$\displaystyle\quad+\|\widehat{\widehat{f}}(t)\|^{2}+\int_{0}^{t}(\|\widehat{f}(s)\|^{2}+{\|\widehat{\widehat{f}}(s)\|^{2})}\,ds\Big{)}.$
(4.16)

Step 2 (Estimate for the term
$\int_{0}^{t}(\mathcal{I}^{1-\alpha}\widehat{\rho},\widehat{\rho})\,ds$ on the
right-hand side of (4.3)). An application of (4.11) results in

$\displaystyle\mathinner{\lVert\widehat{\rho}(s)\rVert}$
$\displaystyle\leq\int_{0}^{s}\|\rho(\tau)\|\,d\tau$
$\displaystyle\lesssim\begin{cases}\displaystyle
h^{1+\gamma_{0}}\mathinner{\lVert
u_{0}\rVert}\int_{0}^{s}\tau^{-\alpha}\,d\tau\lesssim
h^{1+\gamma_{0}}s^{1-\alpha}\mathinner{\lVert u_{0}\rVert}\text{ for
}f=0\text{ and }u_{0}\in L^{2}(\Omega),\\\
h^{1+\gamma_{0}}s\>\Lambda_{0}(\epsilon,s)\text{ for }f\neq 0\text{ and
}u_{0}\in D(A).\end{cases}$

This bound, the definition of $\mathcal{I}^{1-\alpha}$, and $\tau\leq s$ in
the computation of the integral in the second step below lead to

$\displaystyle\mathinner{\lVert\mathcal{I}^{1-\alpha}\widehat{\rho}(s)\rVert}\lesssim\int_{0}^{s}(s-\tau)^{-\alpha}\mathinner{\lVert\widehat{\rho}(\tau)\rVert}\,d\tau\lesssim\begin{cases}h^{1+\gamma_{0}}s^{2-2\alpha}\mathinner{\lVert
u_{0}\rVert}\text{ for }f=0\text{ and }u_{0}\in L^{2}(\Omega),\\\
h^{1+\gamma_{0}}s^{2-\alpha}\Lambda_{0}(\epsilon,s)\text{ for }f\neq 0\text{
and }u_{0}\in D(A).\end{cases}$

The Cauchy-Schwarz inequality and a combination of the last two displayed
bounds show

$\displaystyle\int_{0}^{t}(\mathcal{I}^{1-\alpha}\widehat{\rho},\widehat{\rho})\,ds\leq\int_{0}^{t}\mathinner{\lVert\mathcal{I}^{1-\alpha}\widehat{\rho}\rVert}\mathinner{\lVert\widehat{\rho}\rVert}\,ds\lesssim\begin{cases}h^{2(1+\gamma_{0})}t^{4-3\alpha}\mathinner{\lVert
u_{0}\rVert}^{2}\text{ for }f=0\text{ and }u_{0}\in L^{2}(\Omega),\\\
h^{2(1+\gamma_{0})}t^{4-\alpha}\Lambda_{0}(\epsilon,t)^{2}\text{ for }f\neq
0\text{ and }u_{0}\in D(A).\end{cases}$

Step 3 (Estimate for the term
$(\|\widehat{\widehat{\varphi}}(t)\|^{2}+\int_{0}^{t}(\|\widehat{\varphi}(s)\|^{2}+\|\widehat{\widehat{\varphi}}(s)\|^{2})\,ds)$
on the right-hand side of (4.3)). The approach for the estimates for the
nonsmooth and smooth initial data differs here as direct bounds for the
required estimates offer challenges due to the singularity factor $t^{-1}$
(see (2.17) with $i=1,p=0$) when $u_{0}\in L^{2}(\Omega)$. However, this issue
is resolved below with the help of (4.1).

Nonsmooth initial data ($f=0$ and $u_{0}\in L^{2}(\Omega)$). Recall
$\varphi=\mathcal{I}^{1-\alpha}u^{\prime}$ and apply (4.1) to observe

$\displaystyle\begin{aligned}
\widehat{\varphi}&=\mathcal{I}(\mathcal{I}^{1-\alpha}u^{\prime}(t))=\mathcal{I}^{1-\alpha}\mathcal{I}(u^{\prime}(t))=\mathcal{I}^{1-\alpha}(u(t)-u_{0}),\;\;\widehat{\widehat{\varphi}}(t)=\mathcal{I}^{2-\alpha}(u(t)-u_{0}).\end{aligned}$

The definition of ${\cal I}^{2-\alpha}$ (cf. (1.3)), $(t-s)\leq t$, and (2.17)
with $i=0,p=0$ lead to

$\displaystyle\|\widehat{\widehat{\varphi}}(t)\|^{2}=\|\int_{0}^{t}\frac{(t-s)^{1-\alpha}}{\Gamma(2-\alpha)}(u(s)-u_{0})\,ds\|^{2}\lesssim\Big{(}t^{1-\alpha}\int_{0}^{t}(\|u(s)\|+\|u_{0}\|)\,ds\Big{)}^{2}\lesssim
t^{4-2\alpha}\mathinner{\lVert u_{0}\rVert}^{2}.$

This yields

$\displaystyle\int_{0}^{t}\|\widehat{\widehat{\varphi}}(s)\|^{2}\,ds\lesssim\int_{0}^{t}s^{4-2\alpha}\mathinner{\lVert
u_{0}\rVert}^{2}\,ds=\frac{t^{5-2\alpha}}{5-2\alpha}\mathinner{\lVert
u_{0}\rVert}^{2}\lesssim t^{4-2\alpha}\mathinner{\lVert u_{0}\rVert}^{2}$
(4.17)

with $t\leq T$ (used once) in the last step. The definition of ${\cal
I}^{1-\alpha}$ (cf. (1.3)) and (2.17) with $i=0,p=0$ establish

$\displaystyle\|\widehat{\varphi}(s)\|^{2}=\|\int_{0}^{s}\frac{(s-\tau)^{-\alpha}}{\Gamma(1-\alpha)}(u(\tau)-u_{0})\,d\tau\|^{2}\lesssim\Big{(}\int_{0}^{s}(s-\tau)^{-\alpha}(\|u(\tau)\|+\|u_{0}\|)\,d\tau\Big{)}^{2}\lesssim
s^{2-2\alpha}\mathinner{\lVert u_{0}\rVert}^{2},$

and this shows
$\displaystyle\int_{0}^{t}\|\widehat{\varphi}(s)\|^{2}\,ds\lesssim
t^{3-2\alpha}\mathinner{\lVert u_{0}\rVert}^{2}$. Altogether, we obtain

$\displaystyle\displaystyle\|\widehat{\widehat{\varphi}}(t)\|^{2}+\int_{0}^{t}(\|\widehat{\varphi}(s)\|^{2}+\|\widehat{\widehat{\varphi}}(s)\|^{2})\,ds\lesssim
t^{3-2\alpha}(t+1)\mathinner{\lVert u_{0}\rVert}^{2}.$

Smooth initial data ($f\neq 0$ and $u_{0}\in D(A)$). Using
$\widehat{\widehat{\varphi}}(t)={\cal I}^{3-\alpha}u^{\prime}(t)$, $(t-s)\leq
t$, and (2.21) with $p=0$ and (2.25) with $q=\epsilon$, we have

$\displaystyle\|\widehat{\widehat{\varphi}}(t)\|^{2}$
$\displaystyle\lesssim\Big{(}\int_{0}^{t}(t-s)^{2-\alpha}\mathinner{\lVert
u^{\prime}(s)\rVert}\,ds\Big{)}^{2}\lesssim\Big{(}t^{2-\alpha}\int_{0}^{t}\big{(}s^{\alpha-1}\|u_{0}\|_{D(A)}+s^{-1+\frac{\alpha\epsilon}{2}}\|f(0)\|_{D(A^{\epsilon/2})}$
$\displaystyle\quad+\epsilon^{-1}s^{\frac{\alpha\epsilon}{2}}\|f^{\prime}\|_{L^{\infty}(0,T;D(A^{\epsilon/2}))}\big{)}\,ds\Big{)}^{2}$
$\displaystyle\lesssim\Big{(}t^{2}\mathinner{\lVert
u_{0}\rVert}_{D(A)}+t^{2-\alpha(1-\epsilon/2)}\mathinner{\lVert
f(0)\rVert}_{D(A^{\epsilon/2})}+\epsilon^{-1}t^{3-\alpha(1-\epsilon/2)}\mathinner{\lVert
f^{\prime}\rVert}_{L^{\infty}(0,T;D(A^{\epsilon/2}))}\Big{)}^{2}.$

An integration from $0$ to $t$, (4.8), elementary algebraic manipulations, and
$t\leq T$ (used once) as in (4.17), leads to a similar bound for
$\displaystyle\int_{0}^{t}\|\widehat{\widehat{\varphi}}(s)\|^{2}\,ds$.
Analogous arguments with elementary manipulations lead to

$\displaystyle\int_{0}^{t}\|\widehat{\varphi}(s)\|^{2}\,ds$
$\displaystyle\lesssim\Big{(}t^{\frac{3}{2}}\mathinner{\lVert
u_{0}\rVert}_{D(A)}+t^{\frac{3}{2}-\alpha(1-\frac{\epsilon}{2})}\mathinner{\lVert
f(0)\rVert}_{D(A^{\frac{\epsilon}{2}})}+\epsilon^{-1}t^{\frac{5}{2}-\alpha(1-\frac{\epsilon}{2})}\mathinner{\lVert
f^{\prime}\rVert}_{L^{\infty}(0,T;D(A^{\frac{\epsilon}{2}}))}\Big{)}^{2}.$

A combination of the last two displayed inequalities establishes the required
estimate in this step.

Step 4 (Estimate for the term
$\big{(}\|\widehat{\widehat{f}}(t)\|^{2}+\int_{0}^{t}(\|\widehat{f}(s)\|^{2}+\|\widehat{\widehat{f}}(s)\|^{2})\,ds\big{)}$
on the right-hand side of (4.3)). For nonsmooth initial data, since $f$ is
chosen as zero, the term is zero. For smooth initial data, the definitions of
$\widehat{f}$ and $\widehat{\widehat{f}}$ show

$\displaystyle\|\widehat{\widehat{f}}(t)\|^{2}+\int_{0}^{t}(\|\widehat{f}(s)\|^{2}+\|\widehat{\widehat{f}}(s)\|^{2})\,ds\lesssim\Big{(}t^{3/2}\|f\|_{L^{2}(0,t;L^{2}(\Omega))}+t\|f\|_{L^{2}(0,t;L^{2}(\Omega))}\Big{)}^{2}.$

A combination of Steps 2-4 in (4.3) and algebraic manipulations conclude the
proof. ∎

The notations $\upsilon_{1}(t)=t\upsilon(t)$ for
$\upsilon=\theta,\rho,\varphi,f$ and
$\displaystyle\widehat{\upsilon}(t)=\int_{0}^{t}\upsilon(s)\,ds$ for
$\upsilon=\theta,\theta_{1},\rho,\varphi_{1},f_{1}$ are used in the next
lemma.

######  Lemma 4.6 (Estimate for
$\int_{0}^{t}(\mathcal{I}^{1-\alpha}\theta_{1},\theta_{1})\,ds$).

For $\theta$ that satisfies (4.13) and $t\in(0,T]$, the bounds stated below
hold.

$\displaystyle\int_{0}^{t}(\mathcal{I}^{1-\alpha}\theta_{1},\theta_{1})\,ds+\|\widehat{\theta_{1}}(t)\|_{h}^{2}\lesssim
t^{4-\alpha}\Big{(}h^{1+\gamma_{0}}t^{-\alpha}+h^{2}t^{-(\alpha+1)/2}\Big{)}^{2}\mathinner{\lVert
u_{0}\rVert}^{2},$

if $u_{0}\in L^{2}(\Omega)$ and $f=0$, and

$\displaystyle\int_{0}^{t}(\mathcal{I}^{1-\alpha}\theta_{1},\theta_{1})\,ds+\|\widehat{\theta_{1}}(t)\|_{h}^{2}\lesssim
t^{4-\alpha}\Big{(}h^{1+\gamma_{0}}\Lambda_{0}(\epsilon,t)+h^{2}t^{\alpha/2-2}{\Lambda_{2}(\epsilon,t)}\Big{)}^{2},$

for all $\epsilon\in(0,1)$, if $u_{0}\in D(A)$ and $f$ is such that
$\Lambda_{0}(\epsilon,t)$ and $\Lambda_{2}(\epsilon,t)$, defined in (4.10) and
(4.12), are finite.

###### Proof.

The structure of the proof is similar to that of Lemma 4.5. However, the
majorizations are different and for the readability of the paper, we give the
necessary details. The proof is presented in four steps. To prove the desired
bounds, we proceed following Lemma 4.3 as follows. Step 1 (An intermediate
bound). Multiply (4.13) by $t$, use (4.3), and
$(\rho(0)+\theta(0),\chi_{h})=0$ for all $\chi_{h}\in V_{h}$ to arrive at

$\displaystyle\begin{aligned}
(\mathcal{I}^{1-\alpha}\theta_{1}^{\prime},\chi_{h})+a_{h}(\theta_{1},\chi_{h})&=\alpha(\mathcal{I}^{1-\alpha}\theta,\chi_{h})-(\mathcal{I}^{1-\alpha}\rho_{1}^{\prime},\chi_{h})+\alpha(\mathcal{I}^{1-\alpha}\rho,\chi_{h})\\\
&\quad-(\varphi_{1}(t),(Q-I)\chi_{h})+(f_{1}(t),(Q-I)\chi_{h})\text{ for all
}\chi_{h}\in V_{h}.\end{aligned}$

Integrate the above equality over $(0,t)$, and next choose
$\chi_{h}=\theta_{1}(t)$ to obtain

$\displaystyle\begin{aligned}
(\mathcal{I}^{1-\alpha}\theta_{1},\theta_{1})+a_{h}(\widehat{\theta_{1}},\theta_{1})&=\alpha(\mathcal{I}^{1-\alpha}\widehat{\theta},\theta_{1})-(\mathcal{I}^{1-\alpha}\rho_{1},\theta_{1})+\alpha(\mathcal{I}^{1-\alpha}\widehat{\rho},\theta_{1})\\\
&-(\widehat{\varphi_{1}}(t),(Q-I)\theta_{1})+(\widehat{f_{1}}(t),(Q-I)\theta_{1}).\end{aligned}$

An integration over $(0,t)$ once again, an application of (4.5) to the first
three terms on the right-hand side, and the choice
$\vartheta=\frac{1}{2(2\alpha+1)}$ reveal

$\displaystyle\begin{aligned}
&\int_{0}^{t}(\mathcal{I}^{1-\alpha}\theta_{1},\theta_{1})\,ds+\int_{0}^{t}a_{h}(\widehat{\theta_{1}},\theta_{1})\,ds\leq\frac{1}{2}\int_{0}^{t}(\mathcal{I}^{1-\alpha}\theta_{1},\theta_{1})\,ds+C\int_{0}^{t}\Big{[}(\mathcal{I}^{1-\alpha}\widehat{\theta},\widehat{\theta})+(\mathcal{I}^{1-\alpha}\rho_{1},\rho_{1})\\\
&+(\mathcal{I}^{1-\alpha}\widehat{\rho},\widehat{\rho})\Big{]}\,ds+\Big{|}\int_{0}^{t}(\widehat{\varphi_{1}}(s),(Q-I)\theta_{1})\,ds\Big{|}+\Big{|}\int_{0}^{t}(\widehat{f_{1}}(s),(Q-I)\theta_{1})\,ds\Big{|}.\end{aligned}$

Apply
$2a_{h}(\widehat{\theta_{1}}(t),\theta_{1}(t))=\frac{d}{dt}a_{h}(\widehat{\theta_{1}}(t),\widehat{\theta_{1}}(t))$,
(3.2), approach of Step 1 of Lemma 4.5 to bound the last two terms of the
above displayed estimate, and Gronwall’s lemma to obtain

$\displaystyle\int_{0}^{t}(\mathcal{I}^{1-\alpha}\theta_{1},\theta_{1})\,ds+\|\widehat{\theta_{1}}(t)\|_{h}^{2}\lesssim\int_{0}^{t}(\mathcal{I}^{1-\alpha}\widehat{\theta},\widehat{\theta})\,ds+\int_{0}^{t}(\mathcal{I}^{1-\alpha}\widehat{\rho},\widehat{\rho})\,ds+\int_{0}^{t}(\mathcal{I}^{1-\alpha}\rho_{1},\rho_{1})\,ds$
$\displaystyle\quad+h^{4}\Big{(}\|\widehat{\varphi_{1}}(t)\|^{2}+\int_{0}^{t}(\|\varphi_{1}(s)\|^{2}+\|\widehat{\varphi_{1}}(s)\|^{2})\,ds+\|\widehat{f_{1}}(t)\|^{2}+\int_{0}^{t}(\|f_{1}(s)\|^{2}+\|\widehat{f_{1}}(s)\|^{2})\,ds\Big{)}.$
(4.18)

Observe that the bounds of the first two terms on the right-hand side of (4.3)
are available from the statement and Step 2 of Lemma 4.5. Hence in the steps
below, we bound the remaining terms.

Step 2 (Estimate for the term
$\int_{0}^{t}(\mathcal{I}^{1-\alpha}\rho_{1},\rho_{1})\,ds$ on the right-hand
side of (4.3)). In view of (4.11), we obtain

$\displaystyle\mathinner{\lVert\rho_{1}(s)\rVert}=s\mathinner{\lVert\rho(s)\rVert}\lesssim\begin{cases}\displaystyle
h^{1+\gamma_{0}}s^{1-\alpha}\mathinner{\lVert u_{0}\rVert}\text{ for
}f=0\text{ and }u_{0}\in L^{2}(\Omega),\\\
h^{1+\gamma_{0}}s\>\Lambda_{0}(\epsilon,s)\text{ for }f\neq 0\text{ and
}u_{0}\in D(A).\end{cases}$

Employ this to establish

$\displaystyle\mathinner{\lVert\mathcal{I}^{1-\alpha}\rho_{1}(s)\rVert}\lesssim\int_{0}^{s}(s-\tau)^{-\alpha}\mathinner{\lVert\rho_{1}(\tau)\rVert}\,d\tau\lesssim\begin{cases}h^{1+\gamma_{0}}s^{2-2\alpha}\mathinner{\lVert
u_{0}\rVert}\text{ for }f=0\text{ and }u_{0}\in L^{2}(\Omega),\\\
h^{1+\gamma_{0}}s^{2-\alpha}\Lambda_{0}(\epsilon,s)\text{ for }f\neq 0\text{
and }u_{0}\in D(A).\end{cases}$

The Cauchy-Schwarz inequality plus the last two displayed bounds reveal

$\displaystyle\int_{0}^{t}(\mathcal{I}^{1-\alpha}\rho_{1},\rho_{1})\,ds\leq\int_{0}^{t}\mathinner{\lVert\mathcal{I}^{1-\alpha}\rho_{1}\rVert}\mathinner{\lVert\rho_{1}\rVert}\,ds\lesssim\begin{cases}h^{2(1+\gamma_{0})}t^{4-3\alpha}\mathinner{\lVert
u_{0}\rVert}^{2}\text{ for }f=0\text{ and }u_{0}\in L^{2}(\Omega),\\\
h^{2(1+\gamma_{0})}t^{4-\alpha}\Lambda_{0}(\epsilon,t)^{2}\text{ for }f\neq
0\text{ and }u_{0}\in D(A).\end{cases}$

Step 3 (Estimate for the term
$\big{(}\|\widehat{\varphi_{1}}(t)\|^{2}+\int_{0}^{t}(\|\varphi_{1}(s)\|^{2}+\|\widehat{\varphi_{1}}(s)\|^{2})\,ds\big{)}$
on the right-hand side of (4.3)). Since the nonsmooth data stability estimate
reflects a singularity factor $t^{-1}$ (cf. (2.17) with $i=1,p=0$), we proceed
as follows for this case.

Nonsmooth initial data ($f=0$ and $u_{0}\in L^{2}(\Omega)$). Recall that
$\varphi=\mathcal{I}^{1-\alpha}u^{\prime}$ and apply (4.3) to observe

$\displaystyle\varphi_{1}(s)=s\varphi(s)=s\mathcal{I}^{1-\alpha}u^{\prime}(s)$
$\displaystyle=\mathcal{I}^{1-\alpha}(su(s))^{\prime}-\alpha\mathcal{I}^{1-\alpha}u(s)-s\kappa_{1-\alpha}(s)u_{0}$
$\displaystyle=\mathcal{I}^{1-\alpha}(u(s)+su^{\prime}(s))-\alpha\mathcal{I}^{1-\alpha}u(s)-s\kappa_{1-\alpha}(s)u_{0}.$

The definition of $\mathcal{I}^{1-\alpha}$ in (1.3) and repeated applications
of (2.17) with $i=0,p=0$ and $i=1,p=0$ yield $\|\varphi_{1}(s)\|\lesssim
s^{1-\alpha}\|u_{0}\|$ and hence
$\displaystyle\|\widehat{\varphi_{1}}(t)\|^{2}\lesssim
t^{4-2\alpha}\mathinner{\lVert u_{0}\rVert}^{2}.$ Utilize this bound to arrive
at

$\displaystyle\int_{0}^{t}\|\widehat{\varphi_{1}}(s)\|^{2}\,ds\lesssim\int_{0}^{t}s^{4-2\alpha}\mathinner{\lVert
u_{0}\rVert}^{2}\,ds\lesssim t^{4-2\alpha}\mathinner{\lVert u_{0}\rVert}^{2},$

with $t\leq T$ (used once) in the last step. Analogous arguments show
$\displaystyle\int_{0}^{t}\|\varphi_{1}(s)\|^{2}\,ds\lesssim
t^{3-2\alpha}\mathinner{\lVert u_{0}\rVert}^{2}$. A combination of these
bounds shows

$\displaystyle\displaystyle\|\widehat{\varphi_{1}}(t)\|^{2}+\int_{0}^{t}(\|\varphi_{1}(s)\|^{2}+\|\widehat{\varphi_{1}}(s)\|^{2})\,ds\lesssim
t^{3-2\alpha}(t+1)\mathinner{\lVert u_{0}\rVert}^{2}.$

Smooth initial data ($f\neq 0$ and $u_{0}\in D(A)$). Apply the definition of
$\mathcal{I}^{1-\alpha}$ in (1.3), $s\leq t$, and (2.21) with $p=0$ and (2.25)
with $q=\epsilon$ to obtain

$\displaystyle\|\widehat{\varphi_{1}}(t)\|^{2}$
$\displaystyle=\|\int_{0}^{t}s\mathcal{I}^{1-\alpha}u^{\prime}(s)\,ds\|^{2}\lesssim\Big{(}t\int_{0}^{t}\int_{0}^{s}(s-\tau)^{-\alpha}\|u^{\prime}(\tau)\|\,d\tau\,ds\Big{)}^{2}$
$\displaystyle\lesssim\Big{(}t\int_{0}^{t}\int_{0}^{s}(s-\tau)^{-\alpha}\big{(}\tau^{\alpha-1}\|u_{0}\|_{D(A)}+\tau^{-1+\frac{\alpha\epsilon}{2}}\|f(0)\|_{D(A^{\epsilon/2})}$
$\displaystyle\quad+\epsilon^{-1}\tau^{\frac{\alpha\epsilon}{2}}\|f^{\prime}\|_{L^{\infty}(0,T;D(A^{\epsilon/2}))}\big{)}\,d\tau\,ds\Big{)}^{2}.$
(4.19)

For $\alpha\in(0,1)$, recall the following identities involving the beta
function $B(\cdot,\cdot)$:

$\displaystyle\begin{aligned}
&\int_{0}^{t}(t-s)^{-\alpha}s^{\alpha-1}\,ds=\int_{0}^{1}(1-s)^{\alpha-1}s^{-\alpha}\,ds=B(\alpha,1-\alpha),\\\
&\int_{0}^{t}(t-s)^{-\alpha}s^{\frac{\alpha\epsilon}{2}-1}\,ds=t^{-\alpha(1-\epsilon/2)}\int_{0}^{1}(1-s)^{\frac{\alpha\epsilon}{2}-1}s^{-\alpha}\,ds=B(\alpha\epsilon/2,1-\alpha)t^{-\alpha(1-\epsilon/2)}.\end{aligned}$
(4.20)

Substitute (4.20) in (4.3) to obtain

$\displaystyle\|\widehat{\varphi_{1}}(t)\|^{2}$
$\displaystyle\lesssim\Big{(}t^{2}\mathinner{\lVert
u_{0}\rVert}_{D(A)}+B(\alpha\epsilon/2,1-\alpha)t^{2-\alpha(1-\epsilon/2)}\mathinner{\lVert
f(0)\rVert}_{D(A^{\epsilon/2})}$
$\displaystyle\quad+\epsilon^{-1}t^{3-\alpha(1-\epsilon/2)}\mathinner{\lVert
f^{\prime}\rVert}_{L^{\infty}(0,T;D(A^{\epsilon/2}))}\Big{)}^{2}.$

Utilize $t\leq T$ to conclude that an analogous bound holds for the term
$\int_{0}^{t}\|\widehat{\varphi_{1}}(s)\|^{2}\,ds$. Further, similar
calculations as for $\|\widehat{\varphi_{1}}(t)\|^{2}$ above and algebraic
manipulations with (4.8) reveal

$\displaystyle\int_{0}^{t}\|\varphi_{1}(s)\|^{2}\,ds=\int_{0}^{t}\|s\mathcal{I}^{1-\alpha}u^{\prime}(s)\|^{2}\,ds\lesssim
t^{2}\int_{0}^{t}\Big{(}\int_{0}^{s}(s-\tau)^{-\alpha}\|u^{\prime}(\tau)\|\,d\tau\Big{)}^{2}\,ds\lesssim\Big{(}t^{3/2}\mathinner{\lVert
u_{0}\rVert}_{D(A)}$
$\displaystyle+B(\alpha\epsilon/2,1-\alpha)t^{3/2-\alpha(1-\epsilon/2)}\mathinner{\lVert
f(0)\rVert}_{D(A^{\epsilon/2})}+\epsilon^{-1}t^{5/2-\alpha(1-\epsilon/2)}\mathinner{\lVert
f^{\prime}\rVert}_{L^{\infty}(0,T;D(A^{\epsilon/2}))}\Big{)}^{2}.$

Step 4 (Estimate for the term
$\big{(}\|\widehat{f_{1}}(t)\|^{2}+\int_{0}^{t}(\|f_{1}(s)\|^{2}+\|\widehat{f_{1}}(s)\|^{2})\,ds\big{)}$
on the right-hand side of (4.3)). Since $f=0$ for nonsmooth initial data, the
estimate is trivial. For the case of smooth initial data, the definitions of
$f_{1},\widehat{f_{1}}$ lead to

$\displaystyle\|\widehat{f_{1}}(t)\|^{2}+\int_{0}^{t}(\|f_{1}(s)\|^{2}+\|\widehat{f_{1}}(s)\|^{2})\,ds\lesssim\Big{(}t^{3/2}\|f\|_{L^{2}(0,t;L^{2}(\Omega))}+t\|f\|_{L^{2}(0,t;L^{2}(\Omega))}\Big{)}^{2}.$

A combination of the estimates from Steps 2-4 with (4.3) concludes the proof.
∎

######  Lemma 4.7 (Estimate for
$\int_{0}^{t}(\mathcal{I}^{1-\alpha}\lambda,\lambda)\,ds$).

For
$\lambda=-\rho_{2}^{\prime}+2\alpha\rho_{1}+\alpha(1-\alpha)\widehat{\rho}$
with $\rho_{2}^{\prime}(t)=(t^{2}\rho(t))^{\prime}$, $\rho_{1}(t)=t\rho(t)$,
$\displaystyle\widehat{\rho}(t)=\int_{0}^{t}\rho(s)\,ds$, and $t\in(0,T]$, the
following bounds hold.

$\displaystyle\int_{0}^{t}(\mathcal{I}^{1-\alpha}\lambda,\lambda)\,ds\lesssim
t^{4-\alpha}\Big{(}h^{1+\gamma_{0}}t^{-\alpha}\mathinner{\lVert
u_{0}\rVert}\Big{)}^{2},$

for $u_{0}\in L^{2}(\Omega)$ and $f=0$, and

$\displaystyle\int_{0}^{t}(\mathcal{I}^{1-\alpha}\lambda,\lambda)\,ds\lesssim
t^{4-\alpha}\Big{(}h^{1+\gamma_{0}}\mathcal{B}_{0}(\epsilon,t)\Big{)}^{2},$

for all $\epsilon\in(0,1)$, if $u_{0}\in D(A)$ and $f$ such that
$\mathcal{B}_{0}(\epsilon,t)$ defined in (4.12) is finite.

###### Proof.

Step 1 (Nonsmooth data). For $f=0$ and $u_{0}\in L^{2}(\Omega)$, Lemma 3.4 and
(2.17) (with $i=0,p=1$ and $i=p=1$) result in

$\displaystyle\begin{aligned}
\mathinner{\lVert\lambda(s)\rVert}&\leq\mathinner{\lVert\widehat{\rho}(s)\rVert}+3\mathinner{\lVert\rho_{1}(s)\rVert}+s^{2}\mathinner{\lVert\rho^{\prime}(s)\rVert}\lesssim
h^{1+\gamma_{0}}s^{1-\alpha}\|u_{0}\|+h^{1+\gamma_{0}}s^{2}\mathinner{\lVert
u^{\prime}(s)\rVert}_{H^{2+\gamma_{0}}(\Omega)}\\\ &\lesssim
h^{1+\gamma_{0}}s^{1-\alpha}\|u_{0}\|.\end{aligned}$

This bound leads to

$\displaystyle\begin{aligned}
\mathinner{\lVert\mathcal{I}^{1-\alpha}\lambda(s)\rVert}&\lesssim\int_{0}^{s}(s-\tau)^{-\alpha}\|\lambda(\tau)\|\,d\tau\lesssim
h^{1+\gamma_{0}}s^{2-2\alpha}\|u_{0}\|.\end{aligned}$

The Cauchy-Schwarz inequality and the two bounds displayed above establish

$\displaystyle\begin{aligned}
\int_{0}^{t}(\mathcal{I}^{1-\alpha}\lambda,\lambda)\,ds\leq\int_{0}^{t}\mathinner{\lVert\mathcal{I}^{1-\alpha}\lambda\rVert}\mathinner{\lVert\lambda\rVert}\,ds&\lesssim
h^{2(1+\gamma_{0})}t^{4-3\alpha}\|u_{0}\|^{2}.\end{aligned}$

Step 2 (Smooth data). For $f\neq 0$ and $u_{0}\in D(A)$, apply Lemma 3.4,
(2.20) and (2.24) with $q=\epsilon$, (2.21) with $p=1$ and (2.25) with
$q=\epsilon$ to obtain

$\displaystyle\begin{aligned}
&\mathinner{\lVert\lambda(s)\rVert}\leq\mathinner{\lVert\widehat{\rho}(s)\rVert}+3\mathinner{\lVert\rho_{1}(s)\rVert}+s^{2}\mathinner{\lVert\rho^{\prime}(s)\rVert}\lesssim
h^{1+\gamma_{0}}s\Lambda_{0}(\epsilon,s)+h^{1+\gamma_{0}}s^{2}\mathinner{\lVert
u^{\prime}(s)\rVert}_{H^{2+\gamma_{0}}(\Omega)}\\\ &\lesssim
h^{1+\gamma_{0}}s\Lambda_{0}(\epsilon,s)+h^{1+\gamma_{0}}\Big{(}s\mathinner{\lVert
u_{0}\rVert}_{D(A)}+s^{1+\frac{\alpha\epsilon}{2}}\mathinner{\lVert
f(0)\rVert}_{D(A^{\epsilon/2})}+\epsilon^{-1}s^{2+\frac{\alpha\epsilon}{2}}\mathinner{\lVert
f^{\prime}\rVert}_{L^{\infty}(0,T;D(A^{\epsilon/2})}\Big{)}.\end{aligned}$

An application of this bound in the definition of $\mathcal{I}^{1-\alpha}$
shows

$\displaystyle\begin{aligned}
&\mathinner{\lVert\mathcal{I}^{1-\alpha}\lambda(s)\rVert}{\lesssim\int_{0}^{s}(s-\tau)^{-\alpha}\|\lambda(\tau)\|\,d\tau}\lesssim
h^{1+\gamma_{0}}s^{2-\alpha}\Lambda_{0}(\epsilon,s)\\\
&\quad+h^{1+\gamma_{0}}\Big{(}s^{2-\alpha}\mathinner{\lVert
u_{0}\rVert}_{D(A)}+s^{2-\alpha(1-\epsilon/2}\mathinner{\lVert
f(0)\rVert}_{D(A^{\epsilon/2})}+\epsilon^{-1}s^{3-\alpha(1-\epsilon/2}\mathinner{\lVert
f^{\prime}\rVert}_{L^{\infty}(0,T;D(A^{\epsilon/2})}\Big{)}.\end{aligned}$

The Cauchy-Schwarz inequality and the two bounds displayed above establish

$\displaystyle\begin{aligned}
\int_{0}^{t}(\mathcal{I}^{1-\alpha}\lambda,\lambda)\,ds\leq\int_{0}^{t}\mathinner{\lVert\mathcal{I}^{1-\alpha}\lambda\rVert}\mathinner{\lVert\lambda\rVert}\,ds&\lesssim
h^{2(1+\gamma_{0})}t^{4-\alpha}\Big{(}\Lambda_{0}(\epsilon,t)\\\
&+t^{\alpha\epsilon/2}\mathinner{\lVert
f(0)\rVert}_{D(A^{\epsilon/2})}+\epsilon^{-1}t^{1+\alpha\epsilon/2}\mathinner{\lVert
f^{\prime}\rVert}_{L^{\infty}(0,T;D(A^{\epsilon/2})}\Big{)}^{2}.\end{aligned}$

This concludes the proof. ∎

######  Lemma 4.8 (Estimate for the last term in Lemma 4.3).

For $t\in(0,T]$ and for all $\epsilon\in(0,1)$, it holds that

$\displaystyle\begin{aligned}
h^{4}\Big{(}\|\varphi_{2}(t)\|^{2}+&\int_{0}^{t}(\|\varphi_{2}(s)\|^{2}+\|\varphi_{2}^{\prime}(s)\|^{2})\,ds+\|f_{2}(t)\|^{2}+\int_{0}^{t}(\|f_{2}(s)\|^{2}+\|f_{2}^{\prime}(s)\|^{2})\,ds\Big{)}\\\
&\lesssim
t^{4-\alpha}\begin{cases}\Big{(}h^{2}t^{-(\alpha+1)/2}\mathinner{\lVert
u_{0}\rVert}\Big{)}^{2},\;\;\text{ if }\;u_{0}\in L^{2}(\Omega),f=0,\\\
\Big{(}h^{2}t^{\alpha/2-2}\mathcal{B}_{1}(\epsilon,t)\Big{)}^{2},\;\;\text{ if
}\;u_{0}\in D(A),f\neq 0,\end{cases}\end{aligned}$

where $\varphi_{2}(t)=t^{2}\varphi(t)$,
$\varphi_{2}^{\prime}(s)=(s^{2}\varphi(s))^{\prime}$, $f_{2}(t)=t^{2}f(t)$,
$f_{2}^{\prime}(s)=(s^{2}f(s))^{\prime}$, and $f$ is such that
$\mathcal{B}_{1}(\epsilon,t)$ defined in (4.12) is finite.

###### Proof.

Step 1 (Nonsmooth data). For $f=0$ and $u_{0}\in L^{2}(\Omega)$, we recall
(4.4) to arrive at

$\displaystyle\varphi_{2}(t)$
$\displaystyle=t^{2}\mathcal{I}^{1-\alpha}u^{\prime}(t)=\mathcal{I}^{1-\alpha}(u_{2})^{\prime}(t)-2\alpha\mathcal{I}^{1-\alpha}(u_{1})(t)-\alpha(1-\alpha)\mathcal{I}^{2-\alpha}(u(t))-t^{2}\kappa_{1-\alpha}(t)u_{0}$
$\displaystyle=\mathcal{I}^{1-\alpha}(2tu(t)+t^{2}u^{\prime}(t))-2\alpha\mathcal{I}^{1-\alpha}(tu(t))-\alpha(1-\alpha)\mathcal{I}^{2-\alpha}(u(t))-t^{2}\kappa_{1-\alpha}(t)u_{0}.$

Applications of (1.3) and (2.17) with $i=p=0$ and $i=1,p=0$ in the above
equality yield $\|\varphi_{2}(t)\|\lesssim t^{2-\alpha}\|u_{0}\|$ and this
leads to

$\displaystyle\|\varphi_{2}(t)\|^{2}+\int_{0}^{t}\|\varphi_{2}(s)\|^{2}\,ds\lesssim
t^{4-2\alpha}\|u_{0}\|^{2}.$

To bound the term $\varphi_{2}^{\prime}(s)$, we see that

$\displaystyle\varphi_{2}^{\prime}(s)=2s\varphi(s)+s^{2}\varphi^{\prime}(s)\quad\text{
with }\;\varphi(s)=\mathcal{I}^{1-\alpha}u^{\prime}(s).$ (4.21)

Differentiating $\varphi(s)$ with respect to $s$ we get

$\displaystyle\begin{aligned}
\varphi^{\prime}(s)&=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{s}\tau^{-\alpha}u^{\prime\prime}(s-\tau)\,ds+\frac{s^{-\alpha}u^{\prime}(0)}{\Gamma(1-\alpha)}=\mathcal{I}^{1-\alpha}u^{\prime\prime}(s)+\kappa_{1-\alpha}(s)u^{\prime}(0).\end{aligned}$

Thus, in view of (4.4) we arrive at

$\displaystyle s^{2}\varphi^{\prime}(s)$
$\displaystyle=\mathcal{I}^{1-\alpha}(s^{2}u^{\prime})^{\prime}-2\alpha\mathcal{I}^{1-\alpha}(su^{\prime})-\alpha(1-\alpha)\mathcal{I}^{2-\alpha}(u^{\prime})-s^{2}\kappa_{1-\alpha}(s)u^{\prime}(0)+s^{2}\kappa_{1-\alpha}(s)u^{\prime}(0)$
$\displaystyle=\mathcal{I}^{1-\alpha}(s^{2}u^{\prime})^{\prime}-2\alpha\mathcal{I}^{1-\alpha}(su^{\prime})-\alpha(1-\alpha)\mathcal{I}^{2-\alpha}(u^{\prime}).$

Hence
$\|s^{2}\varphi^{\prime}(s)\|\lesssim\|\mathcal{I}^{1-\alpha}(s^{2}u^{\prime})^{\prime}\|+\|\mathcal{I}^{1-\alpha}(su^{\prime})\|+\|\mathcal{I}^{2-\alpha}(u^{\prime})\|$.
The definition of $\mathcal{I}^{1-\alpha}$ and (2.17) first with $i=1$, $p=0$,
next with $i=0$, $p=0$, and (2.18) lead to

$\displaystyle\begin{aligned}
\|s^{2}\varphi^{\prime}(s)\|&\lesssim\int_{0}^{s}(s-\tau)^{-\alpha}\big{(}\mathinner{\lVert\tau
u^{\prime}(\tau)\rVert}+\mathinner{\lVert\tau^{2}u^{\prime\prime}(\tau)\rVert}+\|u(\tau)\|+\|u_{0}\|\big{)}\,d\tau\lesssim
s^{1-\alpha}\|u_{0}\|.\end{aligned}$

This and Step 3 of Lemma 4.6 in (4.21) show

$\displaystyle\mathinner{\lVert\varphi_{2}^{\prime}(s)\rVert}\lesssim\mathinner{\lVert
s^{2}\varphi^{\prime}(s)\rVert}+\mathinner{\lVert s\varphi(s)\rVert}\lesssim
s^{1-\alpha}\|u_{0}\|.$

This shows
$\displaystyle\int_{0}^{t}\mathinner{\lVert\varphi_{2}^{\prime}(s)\rVert}^{2}\,ds\lesssim
t^{3-2\alpha}\|u_{0}\|^{2}.$

Step 2 (Smooth data). For $f\neq 0$ and $u_{0}\in D(A)$, follow the steps to
bound $\|\widehat{\varphi_{1}}(t)\|^{2}$ for smooth data in the Step 3 in
Lemma 4.6 and recall (4.20) to arrive at

$\displaystyle\|\varphi_{2}(t)\|^{2}\lesssim\Big{(}t^{2}\mathinner{\lVert
u_{0}\rVert}_{D(A)}$
$\displaystyle+B(\alpha\epsilon/2,1-\alpha)t^{2-\alpha(1-\epsilon/2)}\mathinner{\lVert
f(0)\rVert}_{D(A^{\epsilon/2})}$
$\displaystyle\quad+\epsilon^{-1}t^{3-\alpha(1-\epsilon/2)}\mathinner{\lVert
f^{\prime}\rVert}_{L^{\infty}(0,T;D(A^{\epsilon/2}))}\Big{)}^{2}.$

It is easy to establish a similar bound for
$\int_{0}^{t}\|\varphi_{2}(s)\|^{2}\,ds$. The approach in Step 1 of this lemma
and the stability properties in (2.21) with $p=0$, (2.22), (2.25) with
$q=\epsilon$, and (2.26) show

$\displaystyle\int_{0}^{t}\|\varphi_{2}^{\prime}(s)\|^{2}\,ds\lesssim\Big{(}h^{2}t^{\alpha/2-2}B_{1}(\epsilon,t)\Big{)}^{2}.$

The definition of $f_{2}$ and (4.8) establish

$\displaystyle\|f_{2}(t)\|^{2}+\int_{0}^{t}(\|f_{2}(s)\|^{2}+\|f_{2}^{\prime}(s)\|^{2})\,ds\lesssim\Big{(}t^{2}\|f(t)\|+t\|f\|_{L^{2}(0,t;L^{2}(\Omega))}+t^{2}\|f^{\prime}\|_{L^{2}(0,t;L^{2}(\Omega))}\Big{)}^{2}.$

A combination of Steps 1 and 2 establishes the assertion. ∎

We are now ready to prove Theorem 4.1.

Proof of Theorem 4.1. A substitution of the bounds from Lemmas 4.5-4.8 in the
key inequality from Lemma 4.3 and algebraic manipulations establish

$\displaystyle\mathinner{\lVert\theta(t)\rVert}+t^{\alpha/2}\|\theta(t)\|_{h}\lesssim\begin{cases}h^{1+\gamma_{0}}t^{-\alpha}\mathinner{\lVert
u_{0}\rVert}+h^{2}t^{-(\alpha+1)/2}\mathinner{\lVert u_{0}\rVert},\;\;\text{
for }\;u_{0}\in L^{2}(\Omega),f=0,\\\
h^{1+\gamma_{0}}\mathcal{B}_{0}(\epsilon,t)+h^{2}t^{\alpha/2-2}\mathcal{B}_{1}(\epsilon,t),\;\;\text{
for }\;u_{0}\in D(A),f\neq 0,\end{cases}$

and this concludes the proof. $\Box$

##  5 Numerical illustrations

The numerical experiments in this section validate the theoretical orders of
convergences ($OCs$) for the semidiscrete solution established in Section 4.
Let $EOC$ denote the expected orders of convergence with respect to the space
variable. The numerical experiments are performed using Freefem++ [hecht12]
with the following two sets of problem data on $\Omega=(0,1)^{2}$. To compute
the discrete solution, we first triangulate $\overline{\Omega}$, and then
consider a uniform partition of $[0,t]$ with grid points $t_{n}=nk$,
$n=0,1\ldots,N$, where $k=t/N$ is the time step size and $t$ is the time of
interest. The Caputo fractional derivative $\partial_{t}^{\alpha}u(t)$ at
$t=t_{n}$ is approximated by the L1 scheme (cf. [sunwu06, linxu07]) as shown
below:

$\displaystyle\begin{aligned}
&\partial_{t}^{\alpha}u(t_{n})=\frac{1}{\Gamma(1-\alpha)}\sum_{j=0}^{n-1}\int_{t_{j}}^{t_{j+1}}(t_{n}-s)^{-\alpha}\frac{\partial
u(s)}{\partial s}\,ds\\\
&\approx\frac{1}{\Gamma(1-\alpha)}\sum_{j=0}^{n-1}\frac{u(t_{j+1})-u(t_{j})}{k}\int_{t_{j}}^{t_{j+1}}(t_{n}-s)^{-\alpha}\,ds=\sum_{j=0}^{n-1}l_{j}\frac{u(t_{n-j})-u(t_{n-j-1})}{k^{\alpha}}\\\
&=k^{-\alpha}\Big{(}l_{0}u(t_{n})-l_{n-1}u(t_{0})+\sum_{i=1}^{n-1}(l_{j}-l_{j-1})u(t_{n-j})\Big{)},\end{aligned}$

where the weights $l_{j}$ are given by
$l_{j}=((j+1)^{1-\alpha}-j^{1-\alpha})/\Gamma(2-\alpha),\;\;j=0,\ldots,N-1$.

###  5.1 Examples

  * (i)

Choose the nonsmooth initial data as

$\displaystyle u_{0}(x,y)=\begin{cases}1,&x\in(0,1/2],\;y\in(0,1),\\\
-1,&x\in(1/2,1),\;y\in(0,1),\end{cases}$

and $f(x,y,t)=0$. In this case, the exact solution $u(x,y,t)$ is not known.

  * (ii)

The manufactured exact solution $u$ given by
$u(x,y,t)=(t^{\alpha+1}+1)(x(1-x)y(1-y))^{2}$ leads to the smooth initial data
$u_{0}(x,y)=(x(1-x)y(1-y))^{2}$ and source function
$f(x,y,t)=\Gamma(\alpha+2)t(x(1-x)y(1-y))^{2}+(t^{\alpha+1}+1)(24(x^{2}-2x^{3}+x^{4})+24(y^{2}-2y^{3}+y^{4})+2(2-12x+12x^{2})(2-12y+12y^{2}))$.
Note that the initial data in $D(A)$ and source function belongs to the
appropriate function space mentioned in Theorem 4.2 $(ii)$ for any
$\epsilon\in(0,1/4)$.

For the dG method, the numerical computations are performed with the choice of
penalty parameters as $\sigma_{\rm dG}^{1}=\sigma_{\rm dG}^{2}=2$, whereas for
the $C^{0}$IP, we choose $\sigma_{\rm IP}=8$.

###  5.2 Order of spatial convergence for nonsmooth initial data

Since the exact solution is not known in this case, the $OC$ is calculated by
the following formula

$\displaystyle OC=\log(\|W_{2h}-W_{h}\|_{*}/\|W_{h}-W_{h/2}\|_{*})/\log 2,$

where $\|\cdot\|_{*}$ denotes either the $L^{2}(\Omega)$\- or energy norm and
$W_{h}$ is the discrete solution with mesh size $h$. The spatial numerical
experiments are performed for case (i) with the mesh size
$h=\\{1/16,1/32,1/64,1/128,1/256\\}$, and fractional order
$\alpha=0.25,0.50,0.75$ keeping the time step size $k$ fixed at $k=0.001$.

For the Morley method, the numerical results are illustrated in Table 1. The
solution plots for the dG method are displayed in Figure 1 on $256\times 256$
mesh with $\alpha=0.25,0.5,0.75$ with fixed $k=0.001$. Finally, for case
$C^{0}$IP method, the empirical results are shown in Table 2. The energy norm
and $L^{2}$ norm errors demonstrate linear and quadratic orders of
convergence, respectively. These findings agree with the theoretical
convergence rates established in Theorem 4.2 $(i)$.

Table 1: The $L^{2}$-$OC$ and energy-$OC$ for the Morley method in case (i) for $\alpha=0.25,0.50,0.75$ at $t=0.1$ with $k=0.001$ $h$ | $\alpha$ |  errors in $L^{2}$-norm |  $L^{2}$-$OC$ | errors in energy norm |  energy-$OC$
---|---|---|---|---|---
1/16 | 0.25 | 4.50057e-05 | $-$ | 9.30188e-03 | $-$
1/32 |  | 1.18977e-05 | 1.91942 | 4.69592e-03 | 0.98611
1/64 |  | 3.07517e-06 | 1.95195 | 2.34989e-03 | 0.99881
1/128 |  | 7.63686e-07 | 2.00961 | 1.17378e-03 | 1.00143
1/16 | 0.50 | 5.59721e-05 | $-$ | 1.14995e-02 | $-$
1/32 |  | 1.47909e-05 | 1.92000 | 5.80462e-03 | 0.98630
1/64 |  | 3.82260e-06 | 1.95208 | 2.90464e-03 | 0.99884
1/128 |  | 9.49291e-07 | 2.00963 | 1.45088e-03 | 1.00142
1/16 | 0.75 | 5.12578e-05 | $-$ | 1.02222e-02 | $-$
1/32 |  | 1.35183e-05 | 1.92286 | 5.15662e-03 | 0.98720
1/64 |  | 3.49206e-06 | 1.95276 | 2.58010e-03 | 0.99899
1/128 |  | 8.67171e-07 | 2.00969 | 1.28878e-03 | 1.00142
$EOC$ |  |  | 2.0 |  | 1.0

Figure 1: Numerical solutions of dG method on $256\times 256$ mesh with $\alpha=0.25,0.5,0.75$. Table 2: The $L^{2}$-$OC$ and energy-$OC$ for the $C^{0}$IP method in case (i) for $\alpha=0.25,0.50,0.75$ at $t=0.1$ with $k=0.001$ and $\sigma_{\rm IP}=8$ $h$ | $\alpha$ |  errors in $L^{2}$-norm |  $L^{2}$-$OC$ | errors in energy-norm |  energy-$OC$
---|---|---|---|---|---
1/16 | 0.25 | 5.51766e-06 | $-$ | 3.13001e-03 | $-$
1/32 |  | 1.56160e-06 | 1.82103 | 1.59248e-03 | 0.97489
1/64 |  | 4.02967e-07 | 1.95429 | 7.94704e-04 | 1.00278
1/128 |  | 1.05072e-07 | 1.93928 | 4.02570e-04 | 0.98117
1/16 | 0.50 | 6.79515e-06 | $-$ | 3.85380e-03 | $-$
1/32 |  | 1.92315e-06 | 1.82103 | 1.96072e-03 | 0.97489
1/64 |  | 4.96266e-07 | 1.95429 | 9.78472e-04 | 1.00278
1/128 |  | 1.29399e-07 | 1.93929 | 4.95661e-04 | 0.98117
1/16 | 0.75 | 5.92698e-06 | $-$ | 3.35774e-03 | $-$
1/32 |  | 1.67746e-06 | 1.82102 | 1.70834e-03 | 0.97490
1/64 |  | 4.32868e-07 | 1.95428 | 8.52520e-04 | 1.00279
1/128 |  | 1.12867e-07 | 1.93930 | 4.31861e-04 | 0.98116
$EOC$ |  |  | 2.0 |  | 1.0

###  5.3 Order of spatial convergence for smooth data

The spatial numerical experiments are performed with the mesh size
$h=\\{\frac{1}{12},\frac{1}{24},\frac{1}{48},\frac{1}{96}\\}$, and fractional
order $\alpha=0.25,0.50,0.75$ with a fixed time step size $k=0.001$. For the
Morley method, the results are presented in Table 3. The convergence results
for the dG scheme are demonstrated in Table 4. Further, the exact and
numerical solutions plots are also shown in Figure 2. Table 5 and Figure 3
display the numerical results and the exact and numerical solutions for the
$C^{0}$IP method. From the experiments, we observe that the empirical results
are independent of the fractional order $\alpha$. The energy norm and $L^{2}$
norm errors for the Morley, dG, and $C^{0}$IP methods are plotted in a log-log
scale in Figure 4. This validates the theoretical $OCs$ of Theorem 4.2 $(ii)$.

Table 3: The $L^{2}$-$OC$ and energy-$OC$ for the Morley method in case (ii) for $\alpha=0.25,0.50,0.75$ at $t=0.1$ with $k=0.001$ $h$ | $\alpha$ |  errors in $L^{2}$-norm |  $L^{2}$-$OC$ | errors in energy-norm |  energy-$OC$
---|---|---|---|---|---
1/12 | 0.25 | 2.03729e-04 | $-$ | 2.06933e-02 | $-$
1/24 |  | 5.22317e-05 | 1.96365 | 1.04746e-02 | 0.98226
1/48 |  | 1.31474e-05 | 1.99015 | 5.25437e-03 | 0.99531
1/96 |  | 3.29262e-06 | 1.99747 | 2.62936e-03 | 0.99880
1/12 | 0.50 | 1.98919e-04 | $-$ | 2.02103e-02 | $-$
1/24 |  | 5.09999e-05 | 1.96361 | 1.02305e-02 | 0.98222
1/48 |  | 1.28375e-05 | 1.99013 | 5.13193e-03 | 0.99530
1/96 |  | 3.21502e-06 | 1.99746 | 2.56809e-03 | 0.99880
1/12 | 0.75 | 1.96278e-04 | $-$ | 1.99395e-02 | $-$
1/24 |  | 5.03223e-05 | 1.96363 | 1.00933e-02 | 0.98223
1/48 |  | 1.26670e-05 | 1.99013 | 5.06308e-03 | 0.99530
1/96 |  | 3.17243e-06 | 1.99741 | 2.53364e-03 | 0.99880
$EOC$ |  |  | 2.0 |  | 1.0
Table 4: The $L^{2}$-$OC$ and energy-$OC$ for the dG method in case (ii) for $\alpha=0.25,0.50,0.75$ at $t=0.1$ with $k=0.001$ and $\sigma_{\rm dG}^{1}=\sigma_{\rm dG}^{2}=2$ $h$ | $\alpha$ |  errors in $L^{2}$-norm |  $L^{2}$-$OC$ | errors in energy-norm |  energy-$OC$
---|---|---|---|---|---
1/12 | 0.25 | 1.03592e-03 | $-$ | 6.45313e-02 | $-$
1/24 |  | 2.92835e-04 | 1.82275 | 3.05090e-02 | 1.08076
1/48 |  | 7.85626e-05 | 1.89817 | 1.41902e-02 | 1.10434
1/96 |  | 2.03336e-05 | 1.94998 | 6.97876e-03 | 1.02385
1/12 | 0.50 | 1.01130e-03 | $-$ | 6.30133e-02 | $-$
1/24 |  | 2.85915e-04 | 1.82256 | 2.97960e-02 | 1.08054
1/48 |  | 7.67092e-05 | 1.89812 | 1.38593e-02 | 1.10426
1/96 |  | 1.98541e-05 | 1.94996 | 6.81611e-03 | 1.02383
1/12 | 0.75 | 9.97939e-04 | $-$ | 6.21740e-02 | $-$
1/24 |  | 2.82122e-04 | 1.82263 | 2.93972e-02 | 1.08063
1/48 |  | 7.56902e-05 | 1.89814 | 1.36735e-02 | 1.10430
1/96 |  | 1.95904e-05 | 1.94996 | 6.72469e-03 | 1.02384
$EOC$ |  |  | 2.0 |  | 1.0

Figure 2: Exact (left) and numerical (right) solutions of dG method on $96\times 96$ mesh, $\alpha=0.5$. Table 5: The $L^{2}$-$OC$ and energy-$OC$ for the $C^{0}$IP method in case (ii) for $\alpha=0.25,0.50,0.75$ at $t=0.1$ with $k=0.001$ and $\sigma_{\rm IP}=8$ $h$ | $\alpha$ |  errors in $L^{2}$-norm |  $L^{2}$-$OC$ | errors in energy-norm |  energy-$OC$
---|---|---|---|---|---
1/12 | 0.25 | 1.15038e-04 | $-$ | 1.72387e-02 | $-$
1/24 |  | 3.29403e-05 | 1.80419 | 8.87466e-03 | 0.95788
1/48 |  | 8.67545e-06 | 1.92484 | 4.47807e-03 | 0.98681
1/96 |  | 2.21587e-06 | 1.96907 | 2.24645e-03 | 0.99522
1/12 | 0.50 | 1.12328e-04 | $-$ | 1.68369e-02 | $-$
1/24 |  | 3.21638e-05 | 1.80421 | 8.66785e-03 | 0.95788
1/48 |  | 8.47091e-06 | 1.92485 | 4.37372e-03 | 0.98681
1/96 |  | 2.16363e-06 | 1.96906 | 2.19411e-03 | 0.99522
1/12 | 0.75 | 1.10834e-04 | $-$ | 1.66111e-02 | $-$
1/24 |  | 3.17360e-05 | 1.80421 | 8.55157e-03 | 0.95788
1/48 |  | 8.35815e-06 | 1.92487 | 4.31505e-03 | 0.98681
1/96 |  | 2.13470e-06 | 1.96915 | 2.16467e-03 | 0.99522
$EOC$ |  |  | 2.0 |  | 1.0

Figure 3: Exact (left) and numerical (right) solutions of $C^{0}$IP method on
$96\times 96$ mesh, $\alpha=0.5$.

Figure 4: $L^{2}$ errors (left) and energy errors (right) plots for case (ii)
for the Morley, dG, and $C^{0}$IP schemes at $t=0.1$ with $k=0.001$

##  6 Concluding remarks

In this paper, we have studied an initial-boundary value problem for a time-
fractional biharmonic problem with clamped boundary conditions in a bounded
convex polygonal domain in $\mathbb{R}^{2}$. After defining a weak solution,
we have proved the well-posedness result of the problem and derived several
regularity results of the solutions of both homogeneous and nonhomogeneous
problems which are useful in the error analysis. Using an energy argument a
spatially semidiscrete scheme that covers various popular lowest-order
nonstandard piecewise quadratic finite element schemes (e.g., the Morley,
discontinuous Galerkin, and $C^{0}$ interior penalty) is developed and
analyzed. The convergence analysis is carried out for smooth and nonsmooth
initial data cases, including the initial data $u_{0}\in L^{2}(\Omega)$.
Numerical results are provided to validate the theoretical convergence rates
of the discrete solution. One future direction is the convergence analysis for
time discretization, which will be addressed in another work.

Acknowledgements. The second and third authors acknowledge the support from
the IFCAM project “Analysis, Control and Homogenization of Complex Systems”.
The second author acknowledges the support from SERB project on finite element
methods for phase field crystal equation, code (RD/0121-SERBF30-003).

## References

##  Appendix A Appendix

In this section, we provide the details of the interpolation and companion
operators and outline the proof of the error estimates of the Ritz projection
defined in (3.4). For details of the proof of the estimates, we refer to
[cn2022].

######  Lemma A.1 (Morley interpolation and properties [cn2022]).

The (extended) Morley interpolation operator $I_{\rm
M}\mathrel{\mathop{\mathchar 58\relax}}H^{2}(\mathcal{T})\rightarrow{\rm
M}(\mathcal{T})$ is defined in terms of degrees of freedom as follows:

$(I_{\rm M}v_{\mathrm{pw}})(z)\mathrel{\mathop{\mathchar
58\relax}}=|\mathcal{T}(z)|^{-1}\sum_{K\in\mathcal{T}(z)}(v_{\mathrm{pw}}|_{K})(z)\text{
and }\fint_{e}\frac{\partial(I_{\rm
M}v_{\mathrm{pw}})}{\partial\nu}\,ds\mathrel{\mathop{\mathchar
58\relax}}=\fint_{e}\Big{\\{}\\!\\!\\!\Big{\\{}\frac{\partial
v_{\mathrm{pw}}}{\partial\nu}\Big{\\}}\\!\\!\\!\Big{\\}},\;\;v_{\mathrm{pw}}\in
H^{2}(\mathcal{T}),$

where, for an interior vertex $z$, $\mathcal{T}(z)$ denotes the set of
attached triangles with cardinality $|\mathcal{T}(z)|$. The rest of the
degrees of freedom at the vertices and edges on the boundary have no
contribution due to homogeneous boundary conditions. The following properties
hold for $I_{\rm M}$:

$\displaystyle|\\!|\\!|I_{\mathrm{M}}v_{h}|\\!|\\!|_{\mathrm{pw}}$
$\displaystyle\lesssim\|v_{h}\|_{h}\text{ for all }\;v_{h}\in V_{h}$ (A.1)
$\displaystyle\|v_{h}-I_{\mathrm{M}}v_{h}\|_{h}$
$\displaystyle\lesssim\|v_{h}-v\|_{h}\text{ for all }v_{h}\in V_{h}\text{ and
all }v\in V.$ (A.2)

Note that the standard Morley interpolation operator $I_{\rm
M}\mathrel{\mathop{\mathchar 58\relax}}V\rightarrow{\rm M}(\mathcal{T})$
defined by $(I_{\rm M}v)(z)=v(z)$ and $\fint_{e}\frac{\partial I_{\rm
M}v}{\partial\nu}\,ds=\fint_{e}\frac{\partial v}{\partial\nu}\,ds$ for any
$z\in\mathcal{V}(\Omega)$ and $e\in{\mathcal{E}(\Omega)}$ satisfies ([cn2022])

$\displaystyle a_{\rm pw}(v-I_{\rm M}v,w_{h})=0\text{ for all }v\in V\text{
and all }w_{h}\in P_{2}(\mathcal{T}),$ (A.3) $\displaystyle|\\!|\\!|v-I_{\rm
M}v|\\!|\\!|_{\rm pw}\lesssim h^{\gamma}\|v\|_{H^{2+\gamma}(\Omega)}\text{ for
all }v\in H^{2+\gamma}(\Omega),0\leq\gamma\leq 1.$ (A.4)

######  Lemma A.2 (Companion operator and properties [carsput2020, gal2015,
veeserzan2018, cn2022]).

Let ${\rm HCT}(\mathcal{T})$ denote the Hsieh-Clough-Tocher FEM [ciarlet78,
Chapter 6]. There exists a linear mapping $J\mathrel{\mathop{\mathchar
58\relax}}{\rm M}(\mathcal{T})\to({\rm
HCT}(\mathcal{T})+P_{8}(\mathcal{T}))\cap V$ such that any $w_{\rm M}\in{\rm
M}(\mathcal{T})$ satisfies

$\displaystyle\begin{aligned} &\text{(i) }Jw_{\rm M}(z)=w_{\rm M}(z)\text{ for
any }z\in\mathcal{V}\;\text{ (ii) }\nabla({J}w_{\rm
M})(z)=|\mathcal{T}(z)|^{-1}\sum_{K\in\mathcal{T}(z)}(\nabla w_{\rm
M}|_{K})(z)\text{ for }z\in\mathcal{V}(\Omega)\\\ &\text{ (iii)
}\fint_{e}\partial Jw_{\rm M}/\partial\nu\,ds=\fint_{e}\partial w_{\rm
M}/\partial\nu\,ds\text{ for any }e\in\mathcal{E}\;\text{ (iv)
}|\\!|\\!|w_{\rm M}-Jw_{\rm M}|\\!|\\!|_{\rm pw}\lesssim\min_{v\in
V}|\\!|\\!|w_{\rm M}-v|\\!|\\!|_{\rm pw}\\\ &\text{ (v) }w_{\rm M}-Jw_{\rm
M}\perp P_{2}(\mathcal{T})\text{ in }L^{2}(\Omega)\;\text{ (vi) }I_{\rm
M}J=\rm id\text{ in }{\rm M}({\mathcal{T}})\\\ &\text{ (vii)
}\sup_{w_{\mathrm{M}}\in\mathrm{M}(\mathcal{T})\setminus\\{0\\}}|\\!|\\!|w_{\mathrm{M}}-Jw_{\mathrm{M}}|\\!|\\!|_{\mathrm{pw}}/|\\!|\\!|w_{\mathrm{M}}|\\!|\\!|_{\mathrm{pw}}\lesssim
1.\end{aligned}$

######  Lemma A.3 (Transfer operator and its properties [cn2022, cgn2015]).

The linear transfer operator $I_{h}\mathrel{\mathop{\mathchar 58\relax}}{\rm
M}(\mathcal{T})\rightarrow V_{h}$ is chosen as the identity map for the Morley
and dG methods, and is defined for $C^{0}$IP scheme as follows. For all
$w_{\rm M}\in{\rm M}(\mathcal{T})$, for all $z\in\mathcal{V}$, $(I_{h}w_{\rm
M})(z)=w_{\rm M}(z)$, for $z=\text{mid}(e)$ with $e\in\mathcal{E}(\Omega)$,
$(I_{h}w_{\rm M})(z)=\big{\\{}\\!\\!\\!\big{\\{}w_{\rm
M}\big{\\}}\\!\\!\\!\big{\\}}(z)$, and for $z=\text{mid}(e)$ with
$e\in\mathcal{E}(\partial\Omega)$, $(I_{h}w_{\rm M})(z)=0$. It holds that

$\displaystyle\|w_{\rm M}-I_{h}w_{\rm M}\|_{h}\lesssim|\\!|\\!|w_{\rm
M}-v|\\!|\\!|_{\rm pw}\text{ for all }w_{\rm M}\in{\rm M}(\mathcal{T})\text{
and }v\in V.$ (A.5)

######  Lemma A.4 (Relation between $|\\!|\\!|\cdot|\\!|\\!|_{\rm
pw}\mathrel{\mathop{\mathchar 58\relax}}=a_{\rm pw}(\cdot,\cdot)^{1/2}$ and
$\|\cdot\|_{h}$ [cn2022]).

The norms $\|\cdot\|_{h}$ in $H^{2}(\mathcal{T})$ and
$|\\!|\\!|\cdot|\\!|\\!|_{\rm pw}$ in $V+{\rm M}(\mathcal{T})$ are related by

$\displaystyle({\rm i})\;\;|\\!|\\!|\cdot|\\!|\\!|_{\rm
pw}=\|\cdot\|_{h}\text{ in }V+{\rm M}(\mathcal{T})\;\;\text{ and }\;\;({\rm
ii})\;\;|\\!|\\!|\cdot|\\!|\\!|_{\rm pw}\leq\|\cdot\|_{h}\text{ in
}H^{2}(\mathcal{T}).$

The bilinear form $a_{h}(\cdot,\cdot)\mathrel{\mathop{\mathchar
58\relax}}(V_{h}+\mathrm{M}(\mathcal{T}))\times(V_{h}+\mathrm{M}(\mathcal{T}))\rightarrow{\mathbb{R}}$
in (3.3) has the general form

$\displaystyle a_{h}(v_{h},w_{h})\mathrel{\mathop{\mathchar
58\relax}}=a_{\mathrm{pw}}(v_{h},w_{h})+b_{h}(v_{h},w_{h})+c_{h}(v_{h},w_{h})\text{
for all }v_{h},w_{h}\in V_{h}+\mathrm{M}(\mathcal{T}),$ (A.6)

where
$a_{\mathrm{pw}}(\cdot,\cdot),b_{h}(\cdot,\cdot),c_{h}(\cdot,\cdot)\mathrel{\mathop{\mathchar
58\relax}}(V_{h}+\mathrm{M}(\mathcal{T}))\times(V_{h}+\mathrm{M}(\mathcal{T}))\rightarrow{\mathbb{R}}$
are all bounded bilinear forms on $(V_{h}+\mathrm{M}(\mathcal{T}))$. The
sufficient conditions in (H1)-(H4) from [cn2022, Subsections 5.4 and 6.1] help
to establish the approximation properties of the Ritz projection. It is
verified in [cn2022] that all the lowest-order methods considered in this
article satisfy these hypotheses.

(H1)

For all $v_{\mathrm{M}}\in\mathrm{M}(\mathcal{T})$, $w_{h}\in V_{h}$, and
$v,w\in V$,

$\displaystyle
a_{\mathrm{pw}}(v_{\mathrm{M}},w_{h}-I_{\mathrm{M}}w_{h})+b_{h}(v_{\mathrm{M}},w_{h}-I_{\mathrm{M}}w_{h})\lesssim|\\!|\\!|v_{\mathrm{M}}-v|\\!|\\!|_{\mathrm{pw}}\|w_{h}-w\|_{h}.$

(H2)

For all $v_{\mathrm{M}},w_{\mathrm{M}}\in\mathrm{M}(\mathcal{T})$ and
$v_{h},w_{h}\in V_{h}$,

$\displaystyle\text{ (i) }b_{h}(v_{\mathrm{M}},w_{\mathrm{M}})=0\;\;\text{
(ii)
}\;b_{h}(v_{h}+v_{\mathrm{M}},w_{h}+w_{\mathrm{M}})\lesssim\|v_{h}+v_{\mathrm{M}}\|_{h}\|w_{h}+w_{\mathrm{M}}\|_{h}.$

(H3)

For all $v_{h},w_{h}\in V_{h}$ and $v,w\in V$,
$c_{h}(v_{h},w_{h})\lesssim\|v-v_{h}\|_{h}\|w-w_{h}\|_{h}$.

(H4)

For all $v_{h}\in V_{h}$, $w_{\mathrm{M}}\in\mathrm{M}(\mathcal{T})$, and all
$v,w\in V$,

$a_{\mathrm{pw}}(v_{h}-I_{\mathrm{M}}v_{h},w_{\mathrm{M}})+b_{h}(v_{h},w_{\mathrm{M}})\lesssim\|v-v_{h}\|_{h}|\\!|\\!|w-w_{\mathrm{M}}|\\!|\\!|_{\mathrm{pw}}.$

######  Theorem A.5 (Energy norm bound).

Under the hypotheses (H1)–(H3), the Ritz projection $\mathcal{R}_{h}$
satisfies

$\displaystyle\|u-\mathcal{R}_{h}u\|_{h}\lesssim
h^{\gamma}\|u\|_{H^{2+\gamma}(\Omega)},\;\;\gamma\in[0,\gamma_{0}]\text{ with
}\gamma_{0}\in(1/2,1].$

Though the elliptic regularity index for biharmonic problems with clamped
boundary conditions guarantee $\gamma\in[1/2,1)$, Theorems A.5 and A.6 hold
for any $\gamma\in(0,1]$.

###### Proof.

Let $\zeta_{h}\mathrel{\mathop{\mathchar 58\relax}}=I_{h}I_{\rm
M}u-\mathcal{R}_{h}u$ and split the error $u-\mathcal{R}_{h}u$ as a sum of
three terms given by $u-\mathcal{R}_{h}u=(u-I_{\rm M}u)+(I-I_{h})I_{\rm
M}u+\zeta_{h}.$ Lemma A.4 (i) and (A.5) show

$\displaystyle\begin{aligned}
\|u-\mathcal{R}_{h}u\|_{h}&\lesssim|\\!|\\!|u-I_{\mathrm{M}}u|\\!|\\!|_{\mathrm{pw}}+\|\zeta_{h}\|_{h}.\end{aligned}$
(A.7)

For $Q=JI_{\mathrm{M}}$, the coercivity of $a_{h}(\cdot,\cdot)$ from (3.2) and
(3.4) reveal

$\displaystyle\begin{aligned} &\|\zeta_{h}\|_{h}^{2}\lesssim
a_{h}(\zeta_{h},\zeta_{h})=a_{h}(I_{h}I_{\rm
M}u,\zeta_{h})-a_{h}(\mathcal{R}_{h}u,\zeta_{h})=a_{h}(I_{h}I_{\rm
M}u,\zeta_{h})-a(u,Q\zeta_{h})\\\
&=a_{\mathrm{pw}}(u,(I-J)I_{\mathrm{M}}\zeta_{h})-a_{\mathrm{pw}}(I_{\mathrm{M}}u,I_{\mathrm{M}}\zeta_{h})+a_{h}(I_{h}I_{\rm
M}u,\zeta_{h})\end{aligned}$

with (A.3) in the last step. Lemma A.2 (vi) plus (A.3) show
$a_{\mathrm{pw}}(I_{\mathrm{M}}u,(I-J)I_{\mathrm{M}}\zeta_{h})=0$. This,
(A.6), and (H2) (i) in the above displayed expression show

$\displaystyle\begin{aligned} \|\zeta_{h}\|_{h}^{2}&\lesssim
a_{\mathrm{pw}}(u-I_{\mathrm{M}}u,(I-J)I_{\mathrm{M}}\zeta_{h})+(a_{\mathrm{pw}}+b_{h})((I_{h}-I)I_{\mathrm{M}}u,\zeta_{h})\\\
&\quad+((a_{\mathrm{pw}}+b_{h})(I_{\mathrm{M}}u,\zeta_{h}-I_{\mathrm{M}}\zeta_{h})+c_{h}(I_{h}I_{\mathrm{M}}u,\zeta_{h}))=\mathrel{\mathop{\mathchar
58\relax}}T_{1}+T_{2}+T_{3}.\end{aligned}$ (A.8)

The Cauchy-Schwarz inequality, Lemma A.2 (vii), and (A.1) establish

$\displaystyle
T_{1}=a_{\mathrm{pw}}(u-I_{\mathrm{M}}u,(I-J)I_{\mathrm{M}}\zeta_{h})\leq|\\!|\\!|u-I_{\mathrm{M}}u|\\!|\\!|_{\mathrm{pw}}|\\!|\\!|(I-J)I_{\mathrm{M}}\zeta_{h}|\\!|\\!|_{\mathrm{pw}}\lesssim|\\!|\\!|u-I_{\mathrm{M}}u|\\!|\\!|_{\mathrm{pw}}\|\zeta_{h}\|_{h}.$

The Cauchy-Schwarz inequality once again, (H2) (ii), Lemma A.4 (ii), and (A.5)
imply

$\displaystyle T_{2}$
$\displaystyle=a_{\mathrm{pw}}((I_{h}-I)I_{\mathrm{M}}u,\zeta_{h})+b_{h}((I_{h}-I)I_{\mathrm{M}}u,\zeta_{h})\lesssim\|(I-I_{h})I_{\mathrm{M}}u\|_{h}\|\zeta_{h}\|_{h}\lesssim\|u-I_{\mathrm{M}}u\|_{\mathrm{pw}}\|\zeta_{h}\|_{h}.$

To bound the first term in $T_{3}$, apply (H1) with $w=0$ to arrive at

$\displaystyle(a_{\mathrm{pw}}+b_{h})(I_{\mathrm{M}}u,\zeta_{h}-I_{\mathrm{M}}\zeta_{h})\lesssim|\\!|\\!|u-I_{\mathrm{M}}u|\\!|\\!|_{\mathrm{pw}}\|\zeta_{h}\|_{h}.$

Utilize (H3) with $w=0$, the triangle inequality, (A.5), and Lemma A.4 (i) to
obtain

$\displaystyle\begin{aligned}
c_{h}(I_{h}I_{\mathrm{M}}u,\zeta_{h})&\lesssim\|u-I_{h}I_{\mathrm{M}}u\|_{h}\|\zeta_{h}\|_{h}\leq(\|u-I_{\mathrm{M}}u\|_{h}+\|(I-I_{h})I_{\mathrm{M}}u\|_{h})\|\zeta_{h}\|_{h}\lesssim|\\!|\\!|u-I_{\mathrm{M}}u|\\!|\\!|_{\mathrm{pw}}\|\zeta_{h}\|_{h}.\end{aligned}$

A combination of the bounds for $T_{1}$, $T_{2}$, and $T_{3}$ (from last two
displayed inequalities) in (A.8) and then in (A.7), with (A.4) to bound
$|\\!|\\!|u-I_{\mathrm{M}}u|\\!|\\!|_{\mathrm{pw}}$ concludes the proof. ∎

######  Theorem A.6 (Weaker and piecewise Sobolev norm bound).

Under the hypotheses (H1)-(H4), the Ritz projection $\mathcal{R}_{h}$
satisfies

$\displaystyle\mathinner{\lVert
u-\mathcal{R}_{h}u\rVert}_{H^{s}(\mathcal{T})}\lesssim
h^{2-s}\mathinner{\lVert u-\mathcal{R}_{h}u\rVert}_{h}\text{ for }2-\gamma\leq
s\leq 2,0<\gamma\leq 1.$

In particular,

$\|u-\mathcal{R}_{h}u\|\leq\|u-\mathcal{R}_{h}u\|_{H^{1}(\mathcal{T})}\lesssim
h\mathinner{\lVert u-\mathcal{R}_{h}u\rVert}_{h}.$

The proof is based on duality arguments. Here we sketch only the main ideas
and for details, we refer to [cn2022, Section 6]. Note that the $L^{2}$
estimates in the last displayed expression cannot be improved further for
$P_{2}$ approximations [cn2021cmam].

Outline of the proof. For $u\in V$ and $v=Q\mathcal{R}_{h}u\in V$, the duality
of $H^{-s}(\Omega)$ and $H^{s}_{0}(\Omega)$ reveals

$\displaystyle\|u-v\|_{H^{s}(\Omega)}=\sup_{0\neq\mathcal{F}\in
H^{-s}(\Omega)}\frac{\mathcal{F}(u-v)}{\|\mathcal{F}\|_{H^{-s}(\Omega)}}=\mathcal{F}(u-v).$

A corollary of the Hahn-Banach theorem ensures that the supremum is attained
for some $\mathcal{F}\in H^{-s}(\Omega)\subset V^{*}$ with
$\|\mathcal{F}\|_{H^{-s}(\Omega)}=1$. The functional
$a(z,\cdot)=\mathcal{F}\in V^{*}$ has a unique Riesz representation $z\in V$
in the Hilbert space $(V,a)$; $z\in V$ is the weak solution of
$\Delta^{2}z=\mathcal{F}$. By virtue of the elliptic regularity, $z\in V\cap
H^{4-s}(\Omega)$, $2\leq 4-s\leq 2+\gamma$ with
$\|z\|_{H^{4-s}(\Omega)}\lesssim\|\mathcal{F}\|_{H^{-s}(\Omega)}\lesssim 1.$
Hence

$\displaystyle\|u-v\|_{H^{s}(\Omega)}=a(u-v,z)\text{ and
}|\\!|\\!|z-I_{\mathrm{M}}z|\\!|\\!|_{\mathrm{pw}}\lesssim
h^{2-s}\|z\|_{H^{4-s}(\Omega)}\lesssim h^{2-s}.$ (A.9)

For $z_{h}\mathrel{\mathop{\mathchar 58\relax}}=I_{h}I_{\mathrm{M}}z\in V_{h}$
and $\xi\mathrel{\mathop{\mathchar 58\relax}}=Qz_{h}\in V$, the definition of
the Ritz projection $\mathcal{R}_{h}$ (cf. (3.4)) and elementary algebra show
the key identity [cn2022, Lemma 6.3]:

$\displaystyle a(u-v,z)=a(u-v,z-\xi)+a_{\mathrm{pw}}(\mathcal{R}_{h}u-v,\xi-
z_{h})+a_{\mathrm{pw}}((I-J)I_{\mathrm{M}}\mathcal{R}_{h}u,z_{h}-I_{\mathrm{M}}z_{h})$
$\displaystyle+a_{\mathrm{pw}}((I-I_{\mathrm{M}})\mathcal{R}_{h}u,I_{\mathrm{M}}z_{h}-\xi)+a_{h}(\mathcal{R}_{h}u,z_{h})-a_{\mathrm{pw}}(\mathcal{R}_{h}u,I_{\mathrm{M}}z_{h})+a_{\mathrm{pw}}((I-I_{\mathrm{M}})\mathcal{R}_{h}u,z_{h}).$

All the terms on the right-hand side of the identity above can be bounded in
modulus by a constant times $\mathinner{\lVert
u-\mathcal{R}_{h}u\rVert}_{h}|\\!|\\!|z-I_{\mathrm{M}}z|\\!|\\!|_{\mathrm{pw}}$.
The proof of this part is skipped for brevity; we refer to Lemmas 6.4 and 6.5
in [cn2022] for details. Hence with (A.9), we have

$a(u-v,z)\lesssim h^{2-s}\mathinner{\lVert u-\mathcal{R}_{h}u\rVert}_{h}.$

To conclude the proof, split $u-\mathcal{R}_{h}u$ as
$u-\mathcal{R}_{h}u=(u-v)+(v-\mathcal{R}_{h}u)$. Invoke (A.9) and the above
displayed bound to establish

$\|u-v\|_{H^{s}(\Omega)}\lesssim h^{2-s}\|u-\mathcal{R}_{h}u\|_{h}.$

The definition of the Sobolev-Slobodeckii semi-norm shows
$\|u-v\|_{H^{s}(\mathcal{T})}^{2}=\sum_{K\in\mathcal{T}}|u-v|_{H^{s}(K)}^{2}\leq|u-v|_{H^{s}(\Omega)}^{2}\text{
for }\;1<s<2.$ This and the last displayed estimate provide

$\displaystyle\|u-v\|_{H^{s}(\mathcal{T})}\lesssim
h^{2-s}\|u-\mathcal{R}_{h}u\|_{h}\text{ for all }2-\gamma\leq s\leq 2.$

Since $\mathcal{R}_{h}u\in P_{2}(\mathcal{T})$, Lemma 3.5 provides the
estimate $\|v-\mathcal{R}_{h}u\|_{H^{s}(\mathcal{T})}\lesssim
h^{2-s}\|u-\mathcal{R}_{h}u\|_{h}.$ A triangle inequality and a combination of
the last two estimates conclude the proof. $\Box$