\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2024 (2024), No. 58, pp. 1--14.\newline
ISSN: 1072-6691. URL: https://ejde.math.txstate.edu, https://ejde.math.unt.edu \newline
DOI: 10.58997/ejde.2024.58}
\thanks{\copyright 2024. This work is licensed under a CC BY 4.0 license.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2024/58\hfil
 Equations modeling thin film growth]
{Mild solutions to fourth-order parabolic equations modeling 
thin film growth with time fractional derivative}


\author[Q. Liu, W. Zhu, H. Ye \hfil EJDE-2024/58\hfilneg]
{Qiang Liu, Wanyu Zhu, Hailong Ye}

\address{Qiang Liu \newline
School of Mathematical Sciences, 
Shenzhen University, Shenzhen, 518060, China}
\email{matliu@szu.edu.cn}

\address{Wanyu Zhu \newline
School of Mathematical Sciences, 
Shenzhen University, Shenzhen, 518060, China}
\email{2200201025@email.szu.edu.cn} 

\address{Hailong Ye (corresponding author) \newline
School of Mathematical Sciences, 
Shenzhen University, Shenzhen, 518060, China}
\email{yhl@szu.edu.cn}

\thanks{Submitted March 30, 2024. Published October 4, 2024.}
\subjclass[2020]{35G25, 35K90}
\keywords{Thin-film equation; Caputo fractional derivative;
mild solution;
\hfill\break\indent  Mittag-Leffler functions}

\begin{abstract}
In this article, we study initial-boundary problems for  
fourth-order nonlinear parabolic equations modeling thin film growth with 
Caputo-type time fractional derivative. 
By means of the theory of abstract fractional calculus and 
$L^p-L^q$ estimates, we establish the existence and uniqueness of local
mild solutions in the spaces 
$C([0,T]; L^{\frac{\beta N}{2-\beta}}(\Omega))$ with $1<\beta<2$.
Moreover, the local solutions can be extended globally if the initial 
data is sufficiently small.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{proposition}[theorem]{Proposition}
\allowdisplaybreaks

\section{Introduction and main result}

Thin film growth processes play a crucial role in various scientific and 
technological applications, ranging from semiconductor manufacturing to
material sciences \cite{BG1, BG2, BG3}. 
Understanding the dynamics of thin film growth is essential for optimizing 
the quality, stability, and functionality of thin films in these applications. 
Mathematical modeling provides a powerful tool to capture the intricate
dynamics involved in such processes and to develop predictive models that 
guide experimental design and optimization.

In recent years, there has been a growing interest in utilizing fractional 
calculus to describe and analyze complex phenomena that exhibit non-local 
and memory effects. 
In this article, we focus on boundary value problems of 
fourth-order parabolic equation modeling thin film growth with time 
fractional derivative,
\begin{equation}\label{problem}
 \begin{gathered}
 {}_cD_t^{\alpha} u(t)+\Delta^2u=\nabla\cdot f(\nabla u), 
 \quad x\in\Omega,\; t>0,
 \\
 \partial_{\nu}u|_{\partial\Omega}=\partial_{\nu}\Delta u|_{\partial\Omega}=0,
 \quad t>0,
 \\
 u(x,0)=\varphi(x), \quad x\in\Omega.
 \end{gathered}
\end{equation}
Here $\Omega\in\mathbb{R}^{N}$ with $N\ge2$ is a bounded smooth domain, 
$\partial\Omega$ denotes the boundary of $\Omega$ and $\nu$ is the unit 
outer vector  normal to $\Omega$, ${}^cD_t^{\alpha}$ is the Caputo 
derivative with order $\alpha \in (0,1)$,  which is defined by
$$
{}_cD_t^{\alpha} u(t)=\frac{1}{\Gamma(1-\alpha)} 
\frac{d}{d t} \int_0^t(t-s)^{-\alpha}(u(s)-u(0))\,ds,
$$
where the Gamma function defined by 
$\Gamma(\lambda):=\int_0^{\infty} t^{\lambda-1} e^{-t} dt$. 
Here the term $u$ represents the scaled film height, and 
$\Delta^2u$ denotes the capillarity-driven surface diffusion whereas 
$\nabla\cdot f(\nabla u)$ denotes the 
upward hopping of atoms. Throughout this paper, we will assume that 
$f\in C^1(\mathbb{R}^{N}, \mathbb{R}^{N})$ with $f(0)=Df(0)=0$ and 
for some $\beta>1$, 
and $f$ satisfies the following growth condition
\begin{equation}\label{growth condition}
 |f'(\xi_1)-f'(\xi_2)|
 \le C(|\xi_1|^{\beta-1}+|\xi_2|^{\beta-1})|\xi_1-\xi_2|
\end{equation}
for any $\xi_1, \xi_2\in\mathbb{R}^{N}$. As a simple example 
of \eqref{growth condition},
we can take $f(\xi)=|\xi|^{\beta}\xi$.

When $\alpha=1$, the fractional time Caputo derivative is replaced with the 
classical  integer derivative $u_t$. The problem \eqref{problem} becomes 
the usual thin film  growth model, which is well studied by many researchers. 
 King et al.\ \cite{theo1} 
proved the existence, uniqueness, regularity and large time behavior of 
solutions in Sobolev function space. Using Kato's Method, 
Sandjo et al.\ \cite{Sandjo-Moutari-Gningue2015} 
established existence, uniqueness and regularity of the solution in spaces 
of $C^0([0,T]; L^p(\Omega))$ with $p = \frac{n\beta}{2-\beta}$, $1 < \beta 
< 2$. 
Furthermore, they illustrated the qualitative behavior of the approximate 
solution  through some numerical simulations. 
 Ishige et al.\ \cite{theo3} give sufficient 
conditions on the existence of global solutions, the maximal existence 
time and blow up rate for $0<\beta\le 2$. Y. Feng et al.\ \cite{theo4} 
studied the  existence of local mild solutions for any initial data lies 
in $L^2$ on the two-dimensional torus with and without advection.
We refer the interested reader to the references therein.

Fractional models extend the classical model with non-integer orders
differentiation and integration, allowing the incorporation of memory 
and long-range dependencies into models. 
There are many numerical simulations for the time fractional thin film 
growth equations. T. Tang et al.\ \cite{TTao} proved the time-fractional 
molecular  beam epitaxy models admit an energy dissipation law as integral 
cases and proposed a class of finite difference schemes inherited the 
theoretical energy stability.
 Chen et al.\ \cite{numerical3} developed an efficient and accurate, 
full discrete, linear numerical approximation for the time-fractional thin 
film model with the classical Caputo fractional derivative of order $\alpha$ 
and shown  the models possessed an energy dissipation law. 
Wang et al.\ \cite{numerical2} 
proposed a variable-step $L^1$ scheme for the time-fractional molecular beam 
epitaxy model and also investigated the stability and convergence of the
stabilized convex splitting scheme. However, to the best of our knowledge, 
the well-posedness  for the solutions of the time fractional thin film growth 
equation is not clear, which is the main motivation of the present work. 
We refer \cite{wang} for well posedness of the linear and semilinear time 
fractional order Cauchy problem with 
almost sectorial operators, and \cite{ourwork} for well posedness of the 
time fractional order Cahn-Hilliard equation in $\mathbb{R}^3$.

In this article, we focus on the existence and uniqueness of mild solutions 
on problem \eqref{problem}. Because of the observation that the biharmonic 
operator can be regarded as a sectorial operator on some spaces,
we follow some ideas in \cite{Carvalho-Neto2013}, properly adapted to our 
problem.

We denote $\mathcal{A}=\Delta^2$ defined on $L^p(\Omega)$, $1<p<\infty$ 
with its domain
$$
D(\mathcal{A})=\{u\in W^{4,p}(\Omega):
 \partial_{\nu}u|_{\partial\Omega}=\partial_{\nu}\Delta u|_{\partial\Omega}=0
 \}.
$$
It is clear \cite{Sandjo-Dissertation} that the following homogeneous
boundary 
value problem is {\it normally elliptic} in $\Omega$,
\begin{equation}\label{linear problem}
\begin{gathered}
\partial_tu=-\mathcal{A} u, \quad x\in\Omega,\; t>0, \\
u(x,0) = \varphi(x), \quad x\in\Omega.
\end{gathered}
\end{equation}
Hence, the fourth-order operator $-\Delta^2$ with corresponding Neumann 
boundary $(\frac{\partial}{\partial_{\nu}},
 \frac{\partial\Delta}{\partial_{\nu}})$ is the 
infinitesimal generator of an analytic semigroup $T(t)=e^{-t\Delta^2}$ in 
$L^p(\Omega)$. For a detailed discussion on these results, we refer 
the reader to \cite{Amann1990,Amann1993}.

The time fractional version of \eqref{linear problem} with lower order 
term can be written as
\begin{equation}\label{time fractional linear equation}
\begin{gathered} 
\sideset{_c}{^\alpha_t}{\mathop{D}} u =-\mathcal{A} u+F(u), 
\quad x\in \Omega,\; t>0\\
u(x, 0)=\varphi(x), \quad x\in \Omega.
\end{gathered}
\end{equation}
We formally give a mild solution for 
 \eqref{time fractional linear equation}, 
where the rigorous deduction could be found 
in \cite{Carvalho-Neto2013, Kemppainen-SZ2017, liu-ode,wang}. 
We denote by $\mathcal{L}$ Laplace transform operator. 
By the convolution property of Laplace transform, we have
$$
\mathcal{L}({\sideset{_c}{^\alpha_t}{\mathop{D}}u})
=\lambda^{\alpha}\mathcal{L}(u)-\lambda^{\alpha-1}\varphi,
$$
as in \cite{liu-ode}. So taking the Laplace transform to 
\eqref{time fractional linear equation} gives
$$
\mathcal{L}(u) = \lambda^{\alpha-1}(\lambda^\alpha+\mathcal{A})^{-1} 
\varphi + (\lambda^\alpha + \mathcal{A})^{-1} \mathcal{L}(F(u)).
$$
Application of Laplace inversion \cite{Carvalho-Neto2013} implies
$$
u(t)=E_{\alpha}\left(-t^{\alpha} \mathcal{A}\right) \varphi
+\int_0^{t}(t-s)^{\alpha-1} E_{\alpha, \alpha}
\left(-(t-s)^{\alpha} \mathcal{A}\right) F(u)(s) \mathrm{d} s.
$$
where $E_{\alpha}(-t^{\alpha}\mathcal{A})$ and 
$E_{\alpha,\alpha}(-t^{\alpha}\mathcal{A})$ 
are the Mittag-Leffler operators (see Section 2). 
This formal computation then motivates 
the definition of the mild solution of \eqref{problem} as follows:

\begin{definition} \rm
Let $0<\alpha<1$ and $T>0$ .

(i) A function $u$ such that 
$u\in C([0,T]; L^{\frac{\beta N}{2-\beta}}(\Omega))$ defined by
\begin{equation}\label{integro-differential equation}
u(x,t)=E_{\alpha}(-t^{\alpha}\mathcal{A})\varphi
+\int_0^t(t-s)^{\alpha-1}E_{\alpha, \alpha}(-(t-s)^{\alpha}\mathcal{A})
\nabla\cdot f(\nabla u)\,\mathrm{d}s,
\end{equation}
is called a local mild solution of \eqref{problem}.

(ii) If $T=\infty$, we say that $u$ is a global mild solution of 
\eqref{problem}.
\end{definition}


We are now in a position to state the main result of this article.

\begin{theorem}\label{main result}
Suppose $\Omega\subseteq\mathbb{R}^N$ is a bounded domain with $C^4$ boundary, 
$0<\alpha<1<\beta<2$  and $\varphi\in L^{\frac{\beta N}{2-\beta}}(\Omega)$. 
Then there exists $T>0$ such that  \eqref{problem} admits a unique mild 
solution $u\in C([0,T]; L^{\frac{\beta N}{2-\beta}}(\Omega))$  satisfying
$$
 \max\Big\{\sup_{0\le t\le T}t^{\frac{\alpha}{2\beta}
 -\frac{\alpha\gamma}{4\beta^2}}\|\nabla u(t)\|_{\frac{\beta^2 N}{\gamma}},\,
 \sup_{0\le t\le T}t^{\frac{\alpha}{2}}\|\nabla^2 u(t)
  \|_{\frac{\beta N}{2-\beta}}\Big\} <\infty,
$$
where $0<\gamma<\min\{2\beta, \beta(N+1)-2\}$. Furthermore, 
if $\|\varphi\|_{\frac{\beta N}{2-\beta}}$ is sufficiently small,
the solution $u$ can be extended to be global, that is, $T=\infty$.
\end{theorem}

To explain the meaning of the result, we take $f(\xi)=|\xi|^{\beta}\xi$ in 
$\eqref{problem}_1$ as an example. Notice that a smooth function 
$u(x,t)$ solves the equation in \eqref{problem} for $t>0$ if and only if
$$
u_{\lambda}(x,t)=\lambda^{\frac{2\alpha}{\beta}
-\alpha}u(\lambda^{\alpha}x, \lambda^4t)
$$
does so too with each given constant $\lambda$. 
In addition, for the initial data $\varphi$, under the transformation 
$\varphi\mapsto\varphi_{\lambda}$, the $L^{\frac{\beta N}{2-\beta}}(\Omega)$ 
norm is invariant. Therefore, we expect the global existence and uniqueness 
of solutions when the initial data is sufficiently small in the critical 
space $L^{\frac{\beta N}{2-\beta}}(\Omega)$.

It is worth mentioning that the singularity together with
strong nonlinearity arising from the time fractional derivative 
$\sideset{_c}{^\alpha_t}{\mathop{D}} u$ and the nonlinear term 
$\nabla\cdot f(\nabla u)$  make its mathematical analysis more difficult 
in comparison with the fourth-order parabolic equation 
$\eqref{linear problem}_1$.  For example, the Mittag-Leffler 
operators $E_{\alpha}(-t^{\alpha}\mathcal{A})$ and 
$E_{\alpha,\alpha}(-t^{\alpha}\mathcal{A})$
do not satisfy the semigroup properties. So we need to overcome these 
essential difficulties to get some a priori estimates and extend the 
local solution to a  global one. For more details, one can refer 
to Section 3.

This article is organized as follows. In the next section, we introduce some 
elementary properties of Mittag-Leffler operators, which are essential 
throughout the whole paper and give the main results of this paper.

In Section 3, we establish the existence and uniqueness of mild solutions 
using  appropriate functional spaces and Banach fixed pointed theorem.

\section{A priori estimates}

Throughout this paper, $C$ stands for a generic positive constant
which may vary from line to line.
For the analytic semigroup $T(t)$, we have the $L^p-L^q$ estimates, 
which are state as follows.

\begin{proposition}[\cite{Sandjo-Dissertation, Weisler1979}]
\label{L^p-L^q of T(t)}
 Let $1<q\le p<\infty$ and $j=0,1,2,3$. For $u\in L^{q}(\Omega)$, we have
\[
 \|\nabla^jT(t)u\|_{L^p(\Omega)}\le C\, t^{-\frac{N}{4}(\frac{1}{q}
 -\frac{1}{p})-\frac{j}{4}}\|u\|_{L^q(\Omega)},  \quad t>0.
\]
\end{proposition}

Let us recall some properties of Mittag-Leffler operators. 
For $\alpha\in(0,1)$, we denote the entire function 
$M_{\alpha}: \mathbb{C}\to\mathbb{C}$ the Mainardi  function by
$$
M_{\alpha}(z):=\sum_{n=0}^{\infty}\frac{(-z)^n}{n!\Gamma(1-\alpha(1+n))},
$$
which is a particular case of the Wright type function introduced by
Mainardi in \cite{Mainardi1994}  to characterize the fundamental solutions 
for some standard boundary value problems in physics. 
The following classical result gives some essential relations used in 
this article to obtain the main estimates.

\begin{proposition}[\cite{wang}]\label{Mainardi function}
Let $0<\alpha<1$ and $-1<\gamma<\infty$. If we restrict $M_{\alpha}$ to the 
positive real line, then it holds that
$$
M_{\alpha}\in \mathcal{S}([0,\infty)), \quad M_{\alpha}(t)\ge0 \quad
\text{for\ all $t\ge0$ and }
 \int_0^{\infty}t^{\gamma}M_{\alpha}(t)\,\mathrm{d}t
 =\frac{\Gamma(\gamma+1)}{\Gamma(\alpha \gamma+1)},
$$
where $\mathcal{S}([0,\infty))$ is the Schwartz space on $[0, \infty)$.
\end{proposition}

Now, for each $\alpha\in(0,1)$, we define the Mittag-Leffler families 
\begin{gather*}
E_{\alpha}(-t^{\alpha}\mathcal{A})
=\int_0^{\infty}M_{\alpha}(s)T(st^{\alpha})\,\mathrm{d}s, \\
E_{\alpha,\alpha}(-t^{\alpha}\mathcal{A})
=\int_0^{\infty}\alpha sM_{\alpha}(s)T(st^{\alpha})\,\mathrm{d}s.
\end{gather*}
It is interesting to notice
that the Mainardi functions act as a bridge between the fractional and the
classical abstract theories, for more details see 
\cite{wang,Lunardi1995,Carvalho-Neto2013}.

The next result comprises the main assertions about the theory of abstract 
fractional calculus.

\begin{proposition}[\cite{wang}] \label{uniform continuity}
$E_{\alpha}(-t^{\alpha}\mathcal{A})$ and 
$E_{\alpha,\alpha}(-t^{\alpha}\mathcal{A})$ 
are well defined from $L^p(\Omega)$ to $L^p(\Omega)$, $p\in(1, \infty)$. 
Moreover, for $t\ge0$, $E_{\alpha}(-t^{\alpha}\mathcal{A})$ and 
$E_{\alpha,\alpha}(-t^{\alpha}\mathcal{A})$ 
are uniformly continuous in the uniform operator topology on $L^p(\Omega)$.
\end{proposition}

Our proof relies on the well-known Weierstrass M-test in Banach space 
and  $L^p-L^q$  estimates for Mittag-Leffler operators.

\begin{lemma}[Weierstrass M-test \cite{Jeffreys2000}]
\label{Weierstrass M-test}
 Let $X$ denote a Banach space equipped with the norm $\|\cdot\|_{X}$. 
 Suppose $\{\omega_j\}_{j\ge0}$  is a sequence of continuous functions from 
 $[0,T]$ to $X$ such that
 $$
 \sup_{0\le t\le T}\|\omega_j(t)\|_{X}\le M_j, \quad j=0,1,2,\dots
 $$
 where $0<T\le\infty$ and the sequence $\{M_j\}_{j\ge0}$ satisfies 
 $\sum_{j=0}^{\infty}M_j<\infty$, then the sequence $\{\omega_j\}_{j\ge0}$
  converges uniformly on $[0,T]$, that is,
 $$
 \sum_{j=0}^{\infty}\omega_j\in C([0,T]; X).
 $$
\end{lemma}

Using Proposition \ref{L^p-L^q of T(t)}, we can obtain similar 
$L^p-L^q$ estimates for both  families of Mittag-Leffler operators,
 which also depend on the exponent of differentiation $\alpha$.

\begin{proposition}\label{L^p-L^q 2}
Let $1<q\le p<\infty$ and $\frac{1}{q}-\frac{1}{p}<\frac{4-j}{N}$ with any
nonnegative integer $j<4$. For $u\in L^{q}(\Omega)$, we have
\begin{equation}\label{L^p-L^q 2-1}
\|\nabla^jE_{\alpha}(-t^{\alpha}\mathcal{A})u\|_{L^p(\Omega)}
\le C t^{-\frac{N\alpha}{4}(\frac{1}{q}-\frac{1}{p})
 -\frac{j\alpha}{4}}\|u\|_{L^q(\Omega)},\quad t>0.
\end{equation}
\end{proposition}

\begin{proof}
Noting that if
$0\le\frac{1}{q}-\frac{1}{p}<\frac{4-j}{N}$ with $j<4$, we have
$$
-\frac{N}{4}(\frac{1}{q}-\frac{1}{p})-\frac{j}{4}>-1.
$$
By Proposition \ref{L^p-L^q of T(t)} and \ref{Mainardi function}, we obtain
\begin{align*}
\|\nabla^jE_{\alpha}(-t^{\alpha}\mathcal{A})u\|_{L^p(\Omega)}
& \le \int_0^{\infty}M_{\alpha}(s)\|\nabla^jT(st^{\alpha})u\|_{L^p(\Omega)}
\,\mathrm{d}s \\
& \le C\Big(\int_0^{\infty}M_{\alpha}(s)s^{-\frac{N}{4}(\frac{1}{q}
 -\frac{1}{p})-\frac{j}{4}}\,\mathrm{d}s\Big)
 \left(t^{-\frac{N\alpha }{4}(\frac{1}{q}-\frac{1}{p})
  -\frac{j\alpha}{4}}\|u\|_{L^q(\Omega)}\right)
\\
& \le C t^{-\frac{N\alpha}{4}(\frac{1}{q}-\frac{1}{p})
 -\frac{j\alpha}{4}}\|u\|_{L^q(\Omega)}.
\end{align*}
The proof is complete.
\end{proof}

\begin{proposition}\label{L^p-L^q 3}
Let $1<q\le p<\infty$ and $\frac{1}{q}-\frac{1}{p}<\frac{8-j}{N}$ 
with a nonnegative integer $j<4$. Then we have
\[
\|\nabla^jE_{\alpha,\alpha}(-t^{\alpha}\mathcal{A})u\|_{L^p(\Omega)}
\le C t^{-\frac{N\alpha}{4}(\frac{1}{q}-\frac{1}{p})
 -\frac{j\alpha}{4}}\|u\|_{L^q(\Omega)},\quad t>0.
\]
\end{proposition}

The proof of the above propostion is similar to that of 
Proposition \ref{L^p-L^q 2}, so we omit it.

\section{Proof of Theorem \ref{main result}}

In this section, based on the $L^p-L^q$ estimates in section 2, we will 
prove Theorem \ref{main result} by using the method of successive 
approximation.
To simplify the notation, we denote the norm 
$\|\cdot\|_{L^p(\Omega)}$ by $\|\cdot\|_p$.
Recalling the integral equation \eqref{integro-differential equation}, 
we define  a sequence $\{u_j\}_{j\ge0}$ as follows,
\begin{equation}\label{approximation}
\begin{gathered}
u_0(x,t)=E_{\alpha}(-t^{\alpha}\mathcal{A})\varphi,\\
u_j(x,t)=E_{\alpha}(-t^{\alpha}\mathcal{A})\varphi+\int_0^t(t-s)^{\alpha-1}
E_{\alpha, \alpha}(-(t-s)^{\alpha}\mathcal{A})\nabla\cdot f(\nabla u_{j-1})
\,\mathrm{d}s
\end{gathered}
\end{equation}
for positive integers $j$ and $t>0$. In additional, we define
\begin{equation}\label{definition of R_j}
 R(t)_j=\max\Big\{\sup_{0<s\le t}s^{\frac{\alpha}{2}}\|\nabla^2 u_j(s)
 \|_{\frac{\beta N}{2-\beta}},\,
 \sup_{0<s\le t}s^{\frac{\alpha}{2\beta}-\frac{\alpha\gamma}{4\beta^2}}
 \|\nabla u_j(s)\|_{\frac{\beta^2 N}{\gamma}}\Big\}
\end{equation}
for $t>0$ and $j\ge0$, where $0<\gamma<\min\{2\beta, \beta(N+1)-2\}$.

Next, we show  $u_j(t)$ belongs 
 $L^{\frac{\beta N}{2-\beta}}(\Omega)$ 
and is continuous under some smallness conditions. 
The proofs require the following iteration lemma.

\begin{lemma}[\cite{Miao1986}]\label{miao}
 Let $\lambda, \beta>0$ and $b_j$ be a nonnegative sequence such that
 $$
 b_j\le b_0+\lambda b_{j-1}^{1+\beta}
 $$
 for all positive integer $j$. If
 $ 2\lambda(2b_0)^{\beta}<1$,
 then for each nonnegative integer $j$, we have
 $$
 b_j\le\frac{b_0}{1-\lambda(2b_0)^{\beta}}.
 $$
\end{lemma}


\begin{lemma} \label{lemmp1}
For any $T>0$, there exists a constant $\varepsilon_0>0$ independent of 
$T$ such that if $R(T)_0\le\varepsilon_0$, each $u_j(t)$ is well defined 
as an element of $L^{\frac{\beta N}{2-\beta}}(\Omega)$ for any $t>0$, and
\begin{equation}
 R(T)_j\le2R(T)_0, \quad j\ge0.
\end{equation}
\end{lemma}

\begin{proof}
Because $f'(\xi)$ behaves like $|\xi|^{\beta}$, we have
$$
\|f'(\nabla u)\|_{s}\le C\|\nabla u\|_{s}^{\beta}.
$$
By H\"{o}lder's inequality, it is easy to verify that
\begin{equation}\label{nonlinear term}
\begin{aligned}
&\|\nabla\cdot f(\nabla u_{j}(s))\|_{\frac{\beta N}{2-\beta+\gamma}}\\
& = \|\nabla^2u_j(s)\cdot f'(\nabla u_{j}(s))\|_{\frac{\beta N}
{2-\beta+\gamma}} 
\\
& \le \|\nabla^2 u_j(s)\|_{\frac{\beta N}{2-\beta}}
\|f'(\nabla u_j(s))\|_{\frac{\beta N}{\gamma}} 
\\
& = s^{-\alpha+\frac{\alpha\gamma}{4\beta}}\Big(s^{\frac{\alpha}{2}}\|\nabla^2 u_j(s)\|_{\frac{\beta N}{2-\beta}}\Big)
 \Big(s^{\frac{\alpha}{2}-\frac{\alpha\gamma}{4\beta}}
 \|\nabla u_j(s)\|^{\beta}_{\frac{\beta^2 N}{\gamma}}\Big) 
\\
& \le s^{-\alpha+\frac{\alpha\gamma}{4\beta}}
\Big(\sup_{0<s\le t}s^{\frac{\alpha}{2}}\|\nabla^2 u_j(s)
\|_{\frac{\beta N}{2-\beta}}\Big)
 \Big(\sup_{0<s\le t}s^{\frac{\alpha}{2\beta}
 -\frac{\alpha\gamma}{4\beta^2}}\|\nabla u_j(s)
 \|_{\frac{\beta^2 N}{\gamma}}\Big)^{\beta}
\\
& \le s^{-\alpha+\frac{\alpha\gamma}{4\beta}}R(t)_j^{1+\beta}.
\end{aligned}
\end{equation}

Applying the differential operator $\nabla$ to \eqref{approximation}, and 
applying Propositions \ref{L^p-L^q 2} and \ref{L^p-L^q 3}, 
 and \eqref{nonlinear term}, we obtain that
\begin{align}
&\|\nabla u_{j+1}\|_{\frac{\beta^2 N}{\gamma}} \nonumber \\
& \le \|\nabla E_{\alpha}(-t^{\alpha}\mathcal{A})\varphi\|_{\frac{\beta^2 N}{\gamma}}
 +\int_0^t(t-s)^{\alpha-1}\|\nabla E_{\alpha, \alpha}(-(t-s)^{\alpha}\mathcal{A})\nabla\cdot
 f(\nabla u_{j})\|_{\frac{\beta^2 N}{\gamma}}\,\mathrm{d}s 
\nonumber \\
& \le \|\nabla u_0\|_{\frac{\beta^2 N}{\gamma}}
 +C\int_0^t(t-s)^{\alpha-1+\frac{\alpha\gamma}{4\beta^2}-\frac{\alpha}{2\beta}-\frac{\alpha\gamma}{4\beta}}
 \|\nabla\cdot f(\nabla u_{j})\|_{\frac{\beta N}{2-\beta+\gamma}}\,\mathrm{d}s 
\label{Du} \\
& \le \|\nabla u_0\|_{\frac{\beta^2 N}{\gamma}}
 +CR(t)_j^{1+\beta}\int_0^t(t-s)^{\alpha-1+\frac{\alpha\gamma}{4\beta^2}-\frac{\alpha}{2\beta}-\frac{\alpha\gamma}{4\beta}}
 s^{-\alpha+\frac{\alpha\gamma}{4\beta}}\,\mathrm{d}s 
\nonumber \\
& \le \|\nabla u_0\|_{\frac{\beta^2 N}{\gamma}}
 +Ct^{\frac{\alpha\gamma}{4\beta^2}-\frac{\alpha}{2\beta}}R(t)_j^{1+\beta},
\nonumber 
\end{align}
and similarly
\begin{align}
&\|\nabla^2 u_{j+1}\|_{\frac{\beta N}{2-\beta}}\nonumber \\
& \le \|\nabla^2 E_{\alpha}(-t^{\alpha}\mathcal{A})
 \varphi\|_{\frac{\beta N}{2-\beta}}
 +\int_0^t(t-s)^{\alpha-1}\|\nabla^2E_{\alpha, \alpha}
  (-(t-s)^{\alpha}\mathcal{A})\nabla\cdot
 f(\nabla u_{j})\|_{\frac{\beta N}{2-\beta}}\,\mathrm{d}s 
\nonumber \\
& \le \|\nabla^2 u_0\|_{\frac{\beta N}{2-\beta}}
 +C\int_0^t(t-s)^{\frac{\alpha}{2}-1-\frac{\alpha\gamma}{4\beta}}
 \|\nabla\cdot f(\nabla u_{j})
 \|_{\frac{\beta N}{2-\beta+\gamma}}\,\mathrm{d}s 
\nonumber \\
& \le \|\nabla^2 u_0\|_{\frac{\beta N}{2-\beta}}
 +CR(t)_j^{1+\beta}\int_0^t(t-s)^{\frac{\alpha}{2}-1-\frac{\alpha\gamma}
 {4\beta}} s^{-\alpha+\frac{\alpha\gamma}{4\beta}}\,\mathrm{d}s 
\label{D^2u} \\
& \le \|\nabla^2 u_0\|_{\frac{\beta N}{2-\beta}}
+Ct^{-\frac{\alpha}{2}}R(t)_j^{1+\beta}. \nonumber 
\end{align}
Combining \eqref{Du} and \eqref{D^2u}, for any fixed $T>0$, we have
$$
R(T)_{j+1}\le R(T)_0+CR(T)_j^{1+\beta},
$$
where $C>0$ is independent of $T$.
According to Lemma \ref{miao}, there exists a constant $\varepsilon_0>0$ 
such that 
$R(T)_0\le\varepsilon_0$, and $2C(2R(T)_0)^\beta < 1$, which yields
$$
R(T)_j\le2R(T)_0, \quad j\ge0.
$$
At last, we obtain
\begin{align*}
\|u_{j+1}\|_{\frac{\beta N}{2-\beta}}
& \le \|E_{\alpha}(-t^{\alpha}\mathcal{A})
 \varphi\|_{\frac{\beta N}{2-\beta}}
 +\int_0^t(t-s)^{\alpha-1}\|E_{\alpha, \alpha}(-(t-s)^{\alpha}\mathcal{A})
 \nabla\cdot  f(\nabla u_{j})\|_{\frac{\beta N}{2-\beta}}\,\mathrm{d}s 
\\
& \le R(T)_0
 +C\int_0^T(t-s)^{\alpha-1-\frac{\alpha\gamma}{4\beta}}\|\nabla\cdot 
  f(\nabla u_{j}) \|_{\frac{\beta N}{2-\beta+\gamma}}\,\mathrm{d}s 
\\
& \le R(T)_0
 +CR(T)_j^{1+\beta}\int_0^T(t-s)^{\alpha-1-\frac{\alpha\gamma}{4\beta}}
 s^{-\alpha+\frac{\alpha\gamma}{4\beta}}\,\mathrm{d}s 
\\
& \le R(T)_0+CR(T)_j^{1+\beta} 
\\
& \le CR(T)_0,
\end{align*}
which means $u_j(t)\in L^p(\Omega)$ for any $0<t\le T$ and $j\ge0$.
\end{proof}

Next, we shall show the strong continuity of $u_j(t)$ for $t\in[0, T]$.

\begin{lemma}\label{lemmp2} 
Under the assumption $R(T)_0\le\varepsilon_0$ in Lemma \ref{lemmp1}, 
for any $j\ge 0$, we have
$$
 u_{j}\in C([0,T]; L^{\frac{\beta N}{2-\beta}}(\Omega)), 
 \quad T>0,~~ j=0,1,2,\cdots.
$$
\end{lemma}

\begin{proof}
 Fixing $t_0\in(0,T)$, we have the  estimate
\begin{equation} \label{continuity of u_j(t)>0}
\begin{aligned}
&\|u_{j+1}(t)-u_{j+1}(t_0)\|_{\frac{\beta N}{2-\beta}} \\
& \le \|(E_{\alpha}(-t^{\alpha}\mathcal{A})-E_{\alpha}
(-t_0^{\alpha}\mathcal{A}))\varphi\|_{\frac{\beta N}{2-\beta}} 
\\
& \quad  +\int_{t_0}^t(t-s)^{\alpha-1}\left\|E_{\alpha,
  \alpha}(-(t-s)^{\alpha}\mathcal{A})\nabla\cdot
 f(\nabla u_{j})\right\|_{\frac{\beta N}{2-\beta}}\,\mathrm{d}s 
\\
& \quad +\int_0^{t_0}\big\|(t-s)^{\alpha-1}E_{\alpha,
  \alpha}(-(t-s)^{\alpha}\mathcal{A})\nabla\cdot
 f(\nabla u_{j}) 
\\
&\quad -(t_0-s)^{\alpha-1} E_{\alpha, \alpha}(-(t_0-s)^{\alpha}
\mathcal{A})\nabla\cdot
 f(\nabla u_{j})\big\|_{\frac{\beta N}{2-\beta}}\,\mathrm{d}s 
\\
& = I_1+I_2+I_3,
\end{aligned}
\end{equation}
where $t_0<t\le T$.

Using the strong continuity of $E_{\alpha}(-t^{\alpha}\mathcal{A})$ on 
$L^p(\Omega)$ for $t\in[0,\infty)$,
we deduce easily that the first term $I_1$ goes to zero as $t\to t_0^+$.
By Proposition \ref{L^p-L^q 3}, the estimate of the second term yields
\begin{align*}
I_2
& = \int_{t_0}^t(t-s)^{\alpha-1}\left\|E_{\alpha, \alpha}(-(t-s)^{\alpha}\mathcal{A})\nabla\cdot
 f(\nabla u_{j})\right\|_{\frac{\beta N}{2-\beta}}\,\mathrm{d}s
\\
& \le C\int_{t_0}^t(t-s)^{\alpha-1-\frac{\alpha\gamma}{4\beta}}\|\nabla\cdot f(\nabla u_{j})
 \|_{\frac{\beta N}{2-\beta+\gamma}}\,\mathrm{d}s
\\
& \le CR(T)_0^{1+\beta}\int_{t_0}^t(t-s)^{\alpha-1-\frac{\alpha\gamma}{4\beta}}
 s^{-\alpha+\frac{\alpha\gamma}{4\beta}}\,\mathrm{d}s
\\
& \le CR(T)_0^{1+\beta}\Big(\frac{t-t_0}{t_0}\Big)^{1-\alpha+\frac{\alpha\gamma}{4\beta}}.
\end{align*}
Obviously, this term vanishes as $t\to t_0^+$.
For the estimate of the third term, we first denote
\begin{align*}
g(t,s)
&=\big\|(t-s)^{\alpha-1} E_{\alpha, \alpha}(-(t-s)^{\alpha}\mathcal{A})
 \nabla\cdot f(\nabla u_{j})
\\
&\quad-(t_0-s)^{\alpha-1} E_{\alpha, \alpha}(-(t_0-s)^{\alpha}
 \mathcal{A})\nabla\cdot f(\nabla u_{j})\big\|_{\frac{\beta N}{2-\beta}},
\end{align*}
where $0<s<t_0<t$.
Combining this with \eqref{nonlinear term}, it is easy to see that
\begin{equation}
\begin{aligned}
g(t,s)
& \le |(t-s)^{\alpha-1}-(t_0-s)^{\alpha-1}|\cdot\big\| E_{\alpha, \alpha}
 (-(t-s)^{\alpha}\mathcal{A})\nabla\cdot f(\nabla u_{j}(s))
 \big\|_{\frac{\beta N}{2-\beta}} 
 \\
& \quad +(t_0-s)^{\alpha-1}\big\| E_{\alpha,\alpha}(-(t-s)^{\alpha}
 \mathcal{A})  \nabla\cdot f(\nabla u_{j}(s)) 
 \\
& \quad - E_{\alpha, \alpha}(-(t_0-s)^{\alpha}\mathcal{A})\nabla\cdot 
f(\nabla u_{j}(s))\big\|_{\frac{\beta N}{2-\beta}} 
 \\
& \le C|(t-s)^{\alpha-1}-(t_0-s)^{\alpha-1}|(t-s)^{-\frac{\alpha\gamma}
 {4\beta}}\cdot\|\nabla\cdot f(\nabla u_{j}(s))\|_{\frac{\beta N}
 {2-\beta+\gamma}}
 \\
& \quad +(t_0-s)^{\alpha-1}\big\| E_{\alpha,\alpha}
 (-(t-s)^{\alpha}\mathcal{A})\nabla\cdot f(\nabla u_{j}(s))
 \\
& \quad - E_{\alpha, \alpha}(-(t_0-s)^{\alpha}\mathcal{A})\nabla\cdot 
 f(\nabla u_{j}(s))\big\|_{\frac{\beta N}{2-\beta}}
 \\
& \le C|(t-s)^{\alpha-1}-(t_0-s)^{\alpha-1}|(t-s)^{-\frac{\alpha\gamma}
 {4\beta}}s^{-\alpha+\frac{\alpha\gamma}{4\beta}}R(T)_0 
 \\
& \quad +(t_0-s)^{\alpha-1}\big\| E_{\alpha,\alpha}
 (-(t-s)^{\alpha}\mathcal{A})\nabla\cdot f(\nabla u_{j}(s)) 
 \\
& \quad - E_{\alpha, \alpha}(-(t_0-s)^{\alpha}\mathcal{A})\nabla\cdot 
 f(\nabla u_{j}(s))\big\|_{\frac{\beta N}{2-\beta}}. 
\end{aligned}\label{20230901313a}
\end{equation}
Then combining with the strong continuity of 
$E_{\alpha,\alpha}(-t^{\alpha}\mathcal{A})$ on
$L^{\frac{\beta N}{2-\beta}}(\Omega)$ for $t\in(0,\infty)$ , we have
$$
\lim_{t\to t_0^+}g(t,s)=0
$$
for each fixed $s\in(0,t_0)$. In addition, the last term of 
\eqref{20230901313a} is integrable.
Applying dominated convergence theorem yields the third term $I_3$ 
also tends to zero 
as $t\to t_0^+$. Indeed, for $0<s<t_0<t$, we have
\begin{align*}
I_3
& = \int_0^{t_0} g(t,s)\,\mathrm{d}s
\\
& \le C\int_0^{t_0}(t_0-s)^{\alpha-1-\frac{\alpha\gamma}{4\beta}}
 \|\nabla\cdot f(\nabla u_{j}) \|_{\frac{\beta N}
 {2-\beta+\gamma}}\,\mathrm{d}s,
\\
& \le CR(T)_j^{1+\beta}\int_0^{t_0}(t_0-s)^{\alpha-1
 -\frac{\alpha\gamma}{4\beta}}
 s^{-\alpha+\frac{\alpha\gamma}{4\beta}}\,\mathrm{d}s
\\
& \le CR(T)_0^{1+\beta}.
\end{align*}

Similarly, we can also prove the same limit as 
$t\to t_0^-$ with $t_0\in(0, T]$. Thus,
$$
\lim_{t\to t_0}\|u_{j+1}(t)-u_{j+1}(t_0)\|_{\frac{\beta N}{2-\beta}}=0,
\quad t_0\in(0,T].
$$
As for the continuity up to $t=0$ of $u_{j+1}$, Observe that
\begin{align*}
&\|u_{j+1}(t)-u_{j+1}(0)\|_{\frac{\beta N}{2-\beta}}\\
& \le \|E_{\alpha}(-t^{\alpha}\mathcal{A})\varphi
 -\varphi\|_{\frac{\beta N}{2-\beta}}
 +\int_0^t(t-s)^{\alpha-1}\left\| E_{\alpha, \alpha}
  (-(t-s)^{\alpha}\mathcal{A})\nabla\cdot
 f(\nabla u_{j}(s))\right\|_{\frac{\beta N}{2-\beta}}\,\mathrm{d}s 
\\
& \le \|E_{\alpha}(-t^{\alpha}\mathcal{A})\varphi
 -\varphi\|_{\frac{\beta N}{2-\beta}}
 +C\int_0^{t}(t-s)^{\alpha-1-\frac{\alpha\gamma}{4\beta}}\|\nabla\cdot 
  f(\nabla u_{j})
 \|_{\frac{\beta N}{2-\beta+\gamma}}\,\mathrm{d}s r
\\
& \le \|E_{\alpha}(-t^{\alpha}\mathcal{A})\varphi-\varphi
 \|_{\frac{\beta N}{2-\beta}}
 +CR(t)_j^{1+\beta}\int_0^{t}(t-s)^{\alpha-1-\frac{\alpha\gamma}{4\beta}}
 s^{-\alpha+\frac{\alpha\gamma}{4\beta}}\,\mathrm{d}s
\\
& \le \|E_{\alpha}(-t^{\alpha}\mathcal{A})\varphi
 -\varphi\|_{\frac{\beta N}{2-\beta}}
 +CR(t)_0^{1+\beta}.
\end{align*}
Obviously, if
\begin{equation} \label{R small near t=0}
\lim_{t\to0^+}R(t)_0=0,
\end{equation}
then combining with the strong continuity of 
$E_{\alpha}(-t^{\alpha}\mathcal{A})$,
we  obtain that
$$
\lim_{t\to0^+}\|u_{j+1}(t)-u_{j+1}(0)\|_{\frac{\beta N}{2-\beta}} =0.
$$
To prove \eqref{R small near t=0}, we notice that since 
$\varphi\in L^{\frac{\beta N}{2-\beta}}(\Omega)$, 
for any $\varepsilon>0$, there exists 
$\widetilde{\varphi}\in C_0^{\infty}(\Omega)$ such that
$$
\|\varphi-\widetilde{\varphi}\|_{\frac{\beta N}{2-\beta}}<\varepsilon.
$$
Then we have
\begin{align*}
&t^{\frac{\alpha}{2\beta}-\frac{\alpha\gamma}{4\beta^2}}\|\nabla u_0(t)
 \|_{\frac{\beta^2 N}{\gamma}}\\
& = t^{\frac{\alpha}{2\beta}-\frac{\alpha\gamma}{4\beta^2}}
 \|\nabla E_{\alpha}(-t^{\alpha}\mathcal{A})\varphi\|_{\frac{\beta^2 N}
  {\gamma}}
\\
& \le t^{\frac{\alpha}{2\beta}-\frac{\alpha\gamma}{4\beta^2}}\|
 \nabla E_{\alpha}(-t^{\alpha}\mathcal{A})(\varphi-\widetilde{\varphi})
  \|_{\frac{\beta^2 N}{\gamma}}
 +t^{\frac{\alpha}{2\beta}-\frac{\alpha\gamma}{4\beta^2}}\|
  \nabla E_{\alpha}(-t^{\alpha}\mathcal{A})\widetilde{\varphi}
  \|_{\frac{\beta^2 N}{\gamma}}
\\
& \le \|\varphi-\widetilde{\varphi}\|_{\frac{\beta N}{2-\beta}}
 +t^{\frac{\alpha}{4}}\|\nabla \widetilde{\varphi}
 \|_{\frac{\beta N}{2-\beta}}
 \le \varepsilon
\end{align*}
for small $t>0$, which implies
\begin{equation} \label{small near t=0-1}
\lim_{t\to0^+}t^{\frac{\alpha}{2\beta}
-\frac{\alpha\gamma}{4\beta^2}}\|\nabla u_0(t)
\|_{\frac{\beta^2 N}{\gamma}}=0.
\end{equation} 
Similarly, we can  obtain
\begin{equation} \label{small near t=0-2}
\lim_{t\to0^+}t^{\frac{\alpha}{2}}\|\nabla^2 u_0(t)
 \|_{\frac{\beta N}{2-\beta}}=0.
\end{equation}
Combining \eqref{small near t=0-1} with \eqref{small near t=0-2}, 
we derive \eqref{R small near t=0}.
Summing up, we see that
$$
u_{j}\in C([0,T]; L^{\frac{\beta N}{2-\beta}}(\Omega))
$$
for any $T>0$ and $j\ge0$. 
\end{proof}

\begin{proof}[Proof of Theorem \ref{main result}]
To prove the main result,
we only need to show the uniform convergence of the sequence 
$\{u_j\}_{j\ge0}$ under  the assumption that $R(T)_0\le\varepsilon_0$.
By Lemma \ref{Weierstrass M-test}, we define the sequence
\begin{equation}\label{approximation-2}
\begin{gathered}
\omega_0(x,t)=u_0(x,t),\\
\omega_j(x,t)=u_j(x,t)-u_{j-1}(x,t), \quad j\ge1.
\end{gathered}
\end{equation}
Obviously, $\omega_j\in C([0,T]; L^{\frac{\beta N}{2-\beta}}(\Omega))$ for 
any $T>0$ and $j\ge0$.
In addition, we have the identity
\begin{align}
&\omega_{j+1}(x,t) \nonumber\\
& = u_{j+1}(x,t)-u_{j}(x,t) \nonumber\\
& = \int_0^t(t-s)^{\alpha-1}E_{\alpha, \alpha}(-(t-s)^{\alpha}\mathcal{A})
 \Big(\nabla\cdot f(\nabla u_{j})-\nabla\cdot f(\nabla u_{j-1})\Big)\,
 \mathrm{d}s 
\nonumber\\
& = \int_0^t(t-s)^{\alpha-1}E_{\alpha, \alpha}(-(t-s)^{\alpha}\mathcal{A})
 \Big(\nabla^2u_{j}\cdot f'(\nabla u_{j})-\nabla^2u_{j-1}\cdot f'(\nabla u_{j-1})\Big)\,\mathrm{d}s 
\nonumber\\
& = \int_0^t(t-s)^{\alpha-1}E_{\alpha, \alpha}(-(t-s)^{\alpha}\mathcal{A})
 \nabla^2\omega_{j}\cdot f'(\nabla u_{j})\,\mathrm{d}s 
 \label{omega}\\
&\quad +\int_0^t(t-s)^{\alpha-1}E_{\alpha, \alpha}
(-(t-s)^{\alpha}\mathcal{A})
 \nabla^2u_{j-1}\cdot \Big(f'(\nabla u_{j})
 -f'(\nabla u_{j-1})\Big)\,\mathrm{d}s. \nonumber
\end{align}
As in Lemma \ref{lemmp1},  to estimate $\omega_j$, we to derive 
a priori estimates for $\nabla \omega_j$ and $\nabla^2\omega_j$.
So we define
\begin{equation} \label{definition of R_j-2}
\widetilde{R(t)}_j=\max\Big\{\sup_{0<s\le t}s^{\frac{\alpha}{2}}
 \|\nabla^2 \omega_j(s)\|_{\frac{\beta N}{2-\beta}},\,
 \sup_{0<s\le t}s^{\frac{\alpha}{2\beta}-\frac{\alpha\gamma}{4\beta^2}}
 \|\nabla \omega_j(s)\|_{\frac{\beta^2 N}{\gamma}}\Big\}
\end{equation}
for $t>0$ and $j\ge0$, where $0<\gamma<\min\{2\beta, \beta(N+1)-2\}$. 
Obviously, $\widetilde{R(t)}_0=R(t)_0$.

Let $R(T)_0\le\varepsilon_0$. By H\"{o}lder's inequality, we obtain that 
for $0<s\le t\le T$,
\begin{equation} \label{nonlinear term-2}
\begin{aligned}
\|\nabla^2\omega_j(s)\cdot f'(\nabla u_{j}(s))
 \|_{\frac{\beta N}{2-\beta+\gamma}} 
& \le \|\nabla^2 \omega_j(s)\|_{\frac{\beta N}{2-\beta}}
 \|f'(\nabla u_j(s))\|_{\frac{\beta N}{\gamma}} 
\\
& = s^{-\alpha+\frac{\alpha\gamma}{4\beta}}
 \Big(s^{\frac{\alpha}{2}}\|\nabla^2 \omega_j(s)
  \|_{\frac{\beta N}{2-\beta}}\Big)
 \Big(s^{\frac{\alpha}{2}-\frac{\alpha\gamma}{4\beta}}
  \|\nabla u_j(s)\|^{\beta}_{\frac{\beta^2 N}{\gamma}}\Big) 
\\
& \le s^{-\alpha+\frac{\alpha\gamma}{4\beta}}
 \widetilde{R(T)}_jR(T)_j^{\beta} 
\\
& \le Cs^{-\alpha+\frac{\alpha\gamma}{4\beta}}
 \widetilde{R(T)}_jR(T)_0^{\beta}.
\end{aligned}
\end{equation}
In addition, the growth condition on $f'$ and H\"{o}lder's inequality yield, 
for $0<s\le t\le T$,
\begin{align*}
&\|f'(\nabla u_{j}(s))-f'(\nabla u_{j-1}(s))\|_{\frac{\beta N}{\gamma}}\\
& \le \Big\|\nabla \omega_j(s)\big(|\nabla u_{j}(s)|^{\beta-1}
 +|\nabla u_{j-1}(s)|^{\beta-1}\big)\Big\|_{\frac{\beta N}{\gamma}}
\\
& \le \|\nabla \omega_j(s)\|_{\frac{\beta^2 N}{\gamma}} 
 \Big\||\nabla u_{j}(s)|^{\beta-1}
 +|\nabla u_{j-1}(s)|^{\beta-1}\Big\|_{\frac{\beta^2 N}{\gamma(\beta-1)}}
\\
& \le \|\nabla \omega_j(s)\|_{\frac{\beta^2 N}{\gamma}} 
 \Big(\|\nabla u_{j}(s)\|_{\frac{\beta^2 N}{\gamma}}^{\beta-1}
 +\|\nabla u_{j-1}(s)\|_{\frac{\beta^2 N}{\gamma}}^{\beta-1}\Big)
\\
& \le s^{-\frac{\alpha}{2}+\frac{\alpha\gamma}{4\beta}}
 \widetilde{R(T)}_j\Big(R(T)_j^{\beta-1}+R(T)_{j-1}^{\beta-1}\Big).
\end{align*}
Then for $R(T)_0\le\varepsilon_0$, we have
\begin{equation} \label{nonlinear term-3}
\begin{aligned}
& \|\nabla^2u_{j-1}(s)\cdot \big(f'(\nabla u_{j}(s))-f'(\nabla u_{j-1}(s))
\big)\|_{\frac{\beta N}{2-\beta+\gamma}} 
\\
& \le \|\nabla^2u_{j-1}(s)\|_{\frac{\beta N}{2-\beta}}
 \|f'(\nabla u_{j}(s))-f'(\nabla u_{j-1}(s))\|_{\frac{\beta N}{\gamma}} 
\\
& \le s^{-\alpha+\frac{\alpha\gamma}{4\beta}}\widetilde{R(T)}_jR(T)_{j-1}
\Big(R(T)_j^{\beta-1}+R(T)_{j-1}^{\beta-1}\Big)
\\
& \le Cs^{-\alpha+\frac{\alpha\gamma}{4\beta}}
 \widetilde{R(T)}_jR(T)_0^{\beta}.
\end{aligned}
\end{equation}
Applying the differential operator $\nabla$ to \eqref{omega}, and 
combining this with \eqref{nonlinear term-2}--\eqref{nonlinear term-3}, 
we have
\begin{equation}\label{D omega}
\begin{aligned}
&\|\nabla \omega_{j+1}\|_{\frac{\beta^2 N}{\gamma}}\\
& \le \int_0^t(t-s)^{\alpha-1}\|\nabla E_{\alpha, 
 \alpha}(-(t-s)^{\alpha}\mathcal{A})
 \nabla^2\omega_{j}\cdot f'(\nabla u_{j})\|_{\frac{\beta^2 N}{\gamma}}
  \,\mathrm{d}s 
\\
& \quad +\int_0^t(t-s)^{\alpha-1}\|\nabla E_{\alpha, 
 \alpha}(-(t-s)^{\alpha}\mathcal{A})
 \nabla^2u_{j-1}\cdot \big(f'(\nabla u_{j})
  -f'(\nabla u_{j-1})\big)\|_{\frac{\beta^2 N}{\gamma}}\mathrm{d}s 
\\
& \le C\int_0^t(t-s)^{\alpha-1+\frac{\alpha\gamma}{4\beta^2}
 -\frac{\alpha}{2\beta}-\frac{\alpha\gamma}{4\beta}}
 \|\nabla^2\omega_{j}\cdot f'(\nabla u_{j})\|_{\frac{\beta N}
  {2-\beta+\gamma}}\,\mathrm{d}s 
\\
& \quad +C\int_0^t(t-s)^{\alpha-1+\frac{\alpha\gamma}{4\beta^2}
 -\frac{\alpha}{2\beta}-\frac{\alpha\gamma}{4\beta}}
 \|\nabla^2u_{j-1}\cdot \big(f'(\nabla u_{j})-f'(\nabla u_{j-1})
  \big)\|_{\frac{\beta N}{2-\beta+\gamma}}\mathrm{d}s 
\\
& \le C\widetilde{R(T)}_jR(T)_0^{\beta}
 \int_0^t(t-s)^{\alpha-1+\frac{\alpha\gamma}{4\beta^2}
  -\frac{\alpha}{2\beta}-\frac{\alpha\gamma}{4\beta}}
 s^{-\alpha+\frac{\alpha\gamma}{4\beta}}\,\mathrm{d}s 
\\
& \le Ct^{\frac{\alpha\gamma}{4\beta^2}
 -\frac{\alpha}{2\beta}}\widetilde{R(T)}_jR(T)_0^{\beta}, \quad j\ge1.
\end{aligned}
\end{equation}
Similarly, we  have
\begin{align}
&\|\nabla^2 \omega_{j+1}\|_{\frac{\beta N}{2-\beta}}
\nonumber \\
& \le \int_0^t(t-s)^{\alpha-1}\|\nabla^2 E_{\alpha, 
 \alpha}(-(t-s)^{\alpha}\mathcal{A})
 \nabla^2\omega_{j}\cdot f'(\nabla u_{j})\|_{\frac{\beta N}{2-\beta}}\,
  \mathrm{d}s 
\nonumber \\
& \quad  +\int_0^t(t-s)^{\alpha-1}\|\nabla^2 E_{\alpha, \alpha}
 (-(t-s)^{\alpha}\mathcal{A})
 \nabla^2u_{j-1}\cdot \big(f'(\nabla u_{j})-f'(\nabla u_{j-1})
  \big)\|_{\frac{\beta N}{2-\beta}}\mathrm{d}s 
\nonumber \\
& \le C\int_0^t(t-s)^{\frac{\alpha}{2}-1-\frac{\alpha\gamma}{4\beta}}
 \|\nabla^2\omega_{j}\cdot f'(\nabla u_{j})\|_{\frac{\beta N}
  {2-\beta+\gamma}}\,\mathrm{d}s 
\nonumber \\
& \quad  +C\int_0^t(t-s)^{\frac{\alpha}{2}-1-\frac{\alpha\gamma}{4\beta}}
 \|\nabla^2u_{j-1}\cdot \big(f'(\nabla u_{j})
  -f'(\nabla u_{j-1})\big)\|_{\frac{\beta N}{2-\beta+\gamma}}\mathrm{d}s 
\nonumber \\
& \le C\widetilde{R(T)}_jR(T)_0^{\beta}
 \int_0^t(t-s)^{\frac{\alpha}{2}-1-\frac{\alpha\gamma}{4\beta}}
 s^{-\alpha+\frac{\alpha\gamma}{4\beta}}\,\mathrm{d}s 
\nonumber \\
& \le Ct^{-\frac{\alpha}{2}}\widetilde{R(T)}_jR(T)_0^{\beta},\quad j\ge1.
 \label{D^2 omega}
\end{align}
Combining \eqref{definition of R_j-2}, \eqref{D omega} with 
\eqref{D^2 omega}, for any fixed $T>0$, we have
\begin{equation} \label{widetilde{R(T)}j>1}
\widetilde{R(T)}_{j+1}\le C\widetilde{R(T)}_jR(T)_0^{\beta}, 
\quad j=1,2,\dots
\end{equation} 
where $C>0$ is independent of $T$. Noting that
$$
\omega_1=\int_0^t(t-s)^{\alpha-1}
E_{\alpha, \alpha}(-(t-s)^{\alpha}\mathcal{A})\nabla\cdot 
f(\nabla u_0)\,\mathrm{d}s,
$$
following the same idea in the proof \eqref{D omega} and \eqref{D^2 omega},
 we  obtain
\begin{equation} \label{widetilde{R(T)}j=1}
\widetilde{R(T)}_{1}\le C\widetilde{R(T)}_0
R(T)_0^{\beta}=CR(T)_0^{1+\beta}.
\end{equation}
Combining \eqref{widetilde{R(T)}j>1} with \eqref{widetilde{R(T)}j=1},
we finally derive
$$
\widetilde{R(T)}_{j}\le R(T)_0\Big(CR(T)_0^{\beta}\Big)^{j}, \quad 
j=0,1,2,\dots
$$
provided that $R(T)_0\le\varepsilon_0$.

With the help of a priori estimates above, now we can estimate 
$\omega_j$,
\begin{equation} \label{D^0-omega}
\begin{aligned}
&\|\omega_{j+1}\|_{\frac{\beta N}{2-\beta}}\\
& \le \int_0^t(t-s)^{\alpha-1}\|E_{\alpha, \alpha}
 (-(t-s)^{\alpha}\mathcal{A})
 \nabla^2\omega_{j}\cdot f'(\nabla u_{j})\|_{\frac{\beta N}{2-\beta}}
 \,\mathrm{d}s
 \\
& \quad +\int_0^t(t-s)^{\alpha-1}\|E_{\alpha, 
 \alpha}(-(t-s)^{\alpha}\mathcal{A})
 \nabla^2u_{j-1}\cdot \Big(f'(\nabla u_{j})-f'(\nabla u_{j-1})
 \Big)\|_{\frac{\beta N}{2-\beta}}\,\mathrm{d}s
 \\
& \le C\int_0^t(t-s)^{\alpha-1-\frac{\alpha\gamma}{4\beta}}
 \|\nabla^2\omega_{j}\cdot f'(\nabla u_{j})\|_{\frac{\beta N}{2-\beta+\gamma}}
  \,\mathrm{d}s 
\\
& \quad +C\int_0^t(t-s)^{\alpha-1-\frac{\alpha\gamma}{4\beta}}
 \|\nabla^2u_{j-1}\cdot \big(f'(\nabla u_{j})-f'(\nabla u_{j-1})
 \big)\|_{\frac{\beta N}{2-\beta+\gamma}}\mathrm{d}s 
\\
& \le C\widetilde{R(T)}_jR(T)_0^{\beta}\int_0^t(t-s)^{\alpha-1
-\frac{\alpha\gamma}{4\beta}}
 s^{-\alpha+\frac{\alpha\gamma}{4\beta}}\,\mathrm{d}s 
\\
& \le R(T)_0\Big(CR(T)_0^{\beta}\Big)^{j+1}, \quad j\ge1.
\end{aligned}
\end{equation}
Similarly, we obtain that
$$
\|\omega_{1}\|_{\frac{\beta N}{2-\beta}}\le CR(T)_0^{1+\beta}.
$$
Therefore, 
$$
M_j:=\sup_{0\le t\le T}\|\omega_{j}(t)\|_{\frac{\beta N}{2-\beta}}
\le R(T)_0\Big(CR(T)_0^{\beta}\Big)^{j}, \quad j\ge 0
$$
and the sequence $\{M_j\}_{j\ge0}$ is summable provided that
\begin{equation}\label{small R(T)}
R(T)_0<\min\{\varepsilon_0, C^{-\frac{1}{\beta}}\}.
\end{equation}

Since
$\lim_{T\to0^+}R(T)_0=0$,
we can choose $T>0$ small enough such that \eqref{small R(T)} holds.
Then Lemma \ref{Weierstrass M-test} implies that $\{u_j\}_{j\ge0}$ is a 
Cauchy sequence in $C([0,T]; L^{\frac{\beta N}{2-\beta}}(\Omega))$, 
and converges to a unique solution 
$u\in C([0,T]; L^{\frac{\beta N}{2-\beta}}(\Omega))$ of the 
integral equation \eqref{integro-differential equation}.

Moreover, recalling that
$$
u_0(x,t)=E_{\alpha}(-t^{\alpha}\mathcal{A})\varphi,
$$
and combining this with Proposition \ref{L^p-L^q 2}, we obtain
$$
R(T)_0\le \max\Big\{\sup_{0\le t\le T}t^{\frac{\alpha}{2\beta}
 -\frac{\alpha\gamma}{4\beta^2}}\|\nabla u_0(t)\|_{\frac{\beta^2 N}{\gamma}},\,
\sup_{0\le t\le T}t^{\frac{\alpha}{2}}\|\nabla^2 u_0(t)
 \|_{\frac{\beta N}{2-\beta}}\Big\}
\le C\|\varphi\|_{\frac{\beta N}{2-\beta}}.
$$
So if $\|\varphi\|_{\frac{\beta N}{2-\beta}}$ is sufficiently small, 
the solution $u$ can be extended to be global.
The proof of Theorem \ref{main result} is complete.
\end{proof}


\subsection*{Acknowledgements}
This work is supported by the NSFC (12271186, 12271178), 
by the Guangdong Basic and Applied 
Basic Research Foundation (2022B1515120009, 2022A1515010348),
by the Science and Technology 
Program of Shenzhen, China (20231121110406001, 20231120205244001),
and by the China Scholarship Council (202308440378).

\begin{thebibliography}{99}

\bibitem{Amann1990}
H. Amann;
\emph{Dynamic theory of quasilinear parabolic equations II. 
Reaction-diffusion systems}, 
Differential Integral Equations, 4 (1990), 13--75.

\bibitem{Amann1993}
H. Amann;
\emph{Nonhomogeneous linear and quasilinear elliptic and parabolic boundary 
value problem}, 
in: H. T. H. J. Schmeisser (Ed.), Function Spaces,
Differential Operators and Nonlinear Analysis, Teubner, Stuttgart, 
Leipzig, (1993), 9--126.


\bibitem{Carvalho-Neto2013}
P. M. Carvalho-Neto;
\emph{Fractional differential equations: a novel study of local and global 
solutions in Banach spaces}, 
PhD thesis, Universidade de Sao Paulo, Sao Carlos, 2013.

\bibitem{numerical3} 
L. Z. Chen, J. Zhang, J. Zhao, W. X. Cao, H. Wang, J. W. Zhang;
\emph{An accurate and efficient algorithm for the
time-fractional molecular beam epitaxy model with slope selection}, 
Comput. Phys. Commun., 245 (2019), 106842.

\bibitem{theo4} 
Y. Feng, B. Y. Hu, X. Q. X.;
\emph{Suppression of epitaxial thin film growth by mixing}, 
J. Differential Equations, 317 (2022), 561-602.

\bibitem{theo3} 
K. Ishige, N. Miyake, S. Okabe;
\emph{Blowup for a fourth-order parabolic equation with gradient nonlinearity},
SIAM J. Math. Anal., 52 (2020), no. 1, 927--953.

\bibitem{Jeffreys2000}
H. Jeffreys, B. Jeffreys;
\emph{Methods of Mathematical Physics} (3rd edition), 
Cambridge University Press, 2000.

\bibitem{Kemppainen-SZ2017}
J. Kemppainen, J. Siljander, R. Zacher;
\emph{Representation of solutions and large-time behavior for fully nonlocal 
diffusion equations},
J. Differential Equations, 263 (2017), 149--201.

\bibitem{theo1} 
B. B. King, M. Winkler;
\emph{A fourth-order parabolic equation modeling epitaxial thin film growth}, 
J. Math. Anal.Appl., 286 (2003), 459--490.

\bibitem{liu-ode} 
L. Li,  J. G. Liu;
\emph{A generalized definition of caputo derivatives and its application 
to fractional ODEs}, 
SIAM J. Math. Anal., 50 (2018), 2867--2900.

\bibitem{Lunardi1995}
A. Lunardi;
\emph{Analytic Semigroups and Optimal Regularity in Parabolic Problems}, 
Springer-Verlag, 1995.


\bibitem{Mainardi1994}
F. Mainardi;
\emph{On the initial value problem for the fractional diffusion-wave equation}, 
Ser. Adv. Math. Appl. Sci., 23 (1994), 246--251.

\bibitem{Miao1986}
C. Miao;
\emph{Weak solution of class of nonlinear heat equation systems and 
application to the Navier-Stokes system},
J. Differential Equations, 61 (1986), 141--151.

\bibitem{BG1} 
Udo W. Pohl;
\emph{Epitaxy of Semiconductors: Introduction to Physical Principles}. 
Springer-Verlag Berlin Heidelberg, 2013.

\bibitem{Sandjo-Dissertation}
A. N. Sandjo;
\emph{Solutions for fourth-order parabolic equation modeling epitaxial thin 
film growth}, Dissertation, 
RWTH Aachen University, 2011.

\bibitem{Sandjo-Moutari-Gningue2015}
A. N. Sandjo, S. Moutari, Y. Gningue;
\emph{Solutions of fourth-order parabolic equation modeling thin film growth}, 
J. Differential Equations, 259 (2015), 7260--7283.


\bibitem{BG2} 
R. L. Schwoebel, E. J. Shipsey Step;
\emph{Motion on crystal surfaces}, 
J. Appl. Phys., 37 (1966), 3682--3686. 

\bibitem{TTao} 
T. Tang, H. Yu, T. Zhou;
\emph{On energy dissipation theory and numerical stability for 
time-fractional phase-field equations},
SIAM J. Sci. Comput., 41(6)(2019), 3757--3778.


\bibitem{wang} 
R. N. Wang, D. H. Chen, T. J. Xiao;
\emph{Abstract fractional Cauchy problems with almost sectorial operators}, 
J. Differential Equations, 252 (2012), 202--235.

\bibitem{numerical2} 
J. D. Wang, Y. Yang, B. Q. Ji;
\emph{Two energy stable variable-step L1 schemes for the time-fractional 
MBE model without slope selection}, 
J. Comput. Appl. Math., 419 (2023), 1--15.

\bibitem{Weisler1979}
F. B. Weisler;
\emph{Semilinear evolution equations in Banach spaces}, 
J. Funct. Anal., 32 (1979), 277--296.

\bibitem{ourwork} 
H. L. Ye, Q. Liu, Z. M. Chen;
\emph{Global existence of solutions of the time fractional Cahn-Hilliard 
equation in $\mathbb{R}^3$},  J. Evol. Equ., 21 (2021), no.2, 2377--2411.

\bibitem{BG3} 
A. Zangwill;
\emph{Some causes and a consequence of epitaxial roughening}, 
J. Cryst. Growth, 163 (1996), 8--21.

\end{thebibliography}

\end{document}
