Discrete maximal regularity of time-stepping schemes for fractional evolution equations

In this work, we establish the maximal \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell ^p$$\end{document}ℓp-regularity for several time stepping schemes for a fractional evolution model, which involves a fractional derivative of order \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \in (0,2)$$\end{document}α∈(0,2), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \ne 1$$\end{document}α≠1, in time. These schemes include convolution quadratures generated by backward Euler method and second-order backward difference formula, the L1 scheme, explicit Euler method and a fractional variant of the Crank–Nicolson method. The main tools for the analysis include operator-valued Fourier multiplier theorem due to Weis (Math Ann 319:735–758, 2001. doi:10.1007/PL00004457) and its discrete analogue due to Blunck (Stud Math 146:157–176, 2001. doi:10.4064/sm146-2-3). These results generalize the corresponding results for parabolic problems.


Introduction
Maximal L p -regularity is an important mathematical tool in studying the existence, uniqueness and regularity of solutions of nonlinear partial differential equations of parabolic type. A generator A of an analytic semigroup on a Banach space X is said to have maximal L p -regularity, if the solution u of the following parabolic differential equation satisfies the following estimate u L p (R + ;X ) + Au L p (R + ;X ) ≤ c p,X f L p (R + ;X ) ∀ f ∈ L p (R + ; X ), (1.2) with 1 < p < ∞. On a Hilbert space X , every generator of a bounded analytic semigroup has maximal L p -regularity [45], and Hilbert spaces are only spaces for which this holds true [24]. Beyond Hilbert spaces, an important and very useful characterization of the maximal L p -regularity was given by Weis [48] on UMD spaces in terms of the R-boundedness of a family of operators using the resolvent R(z; A) := (z − A) −1 ; see Theorem 1 in Sect. 2 for details. An important question from the perspective of numerical analysis is whether such maximal regularity estimates carry over to time-stepping schemes for discretizing the parabolic problem (1.1), which have important applications in numerical analysis of nonlinear parabolic problems [1,2,17,29,33]. This question has been studied in a number of works from different aspects [4,5,15,16,31,32,34]. Ashyralyev et al. [4] showed the following discrete maximal regularity: for all f n ∈ X, n = 1, 2, . . . , for the time-discrete solutions u n , n = 1, 2, . . . , given by the implicit Euler method, where τ is the time step size and the constant c p,X is independent of τ . A variant of the maximal p -regularity for the Crank-Nicolson method was also shown in [4].
Recently, Kovács et al. [28] proved the discrete maximal regularity for the Crank-Nicolson, BDF and A-stable Runge-Kutta methods. Kemmochi and Saito [25,26] proved the maximal p -regularity for the θ -method. In these works, the main tools are the maximal L p -regularity characterization due to Weis [48] and its discrete analogue due to Blunck [10]. Independently, Leykekhman and Vexler [31] proved the maximal L p -regularity of discontinuous Galerkin methods without using Blunck's multiplier technique. The maximal p -regularity of fully discrete numerical solutions have been investigated in [25,26,31] and [35] for parabolic equations with time-independent and time-dependent coefficients, respectively; also see [28, section 6]. The maximal L p -regularity has also been studied for the following fractional evolution equation ∂ α t u(t) = Au(t) + f ∀t > 0, (1.3) together with the following initial condition(s) In the model (1.3), the notation ∂ α t u denotes the Caputo fractional derivative of order α of u with respect to time t, defined by [27, pp. 91]  Throughout, we use only the notation ∂ α t u to denote either derivative. When α = 1, the fractional derivative ∂ α t u(t) coincides with the usual first-order derivative u (t), and accordingly, the fractional model (1.3) recovers the standard parabolic equation (1.1). In this paper we focus on the fractional cases 0 < α < 1 and 1 < α < 2, which are known as the subdiffusion and diffusion-wave equation, respectively. In analogy with Brownian motion for normal diffusion (1.1), the model (1.3) with 0 < α < 1 is the macroscopic counterpart of continuous time random walk.
The fractional model (1.3) has received much attention in recent years, since it can adequately capture the dynamics of anomalous diffusion processes. For example, the subdiffusion equation, i.e., α ∈ (0, 1), has been employed to describe transport in column experiments, thermal diffusion in media with fractal geometry, and flow in highly heterogeneous aquifers. See [40] for an extensive list of applications. The diffusion-wave equation, i.e., α ∈ (1, 2), can be used to model mechanical wave propagation in viscoelastic media.
In a series of interesting works [6][7][8], Bazhalekov and collaborators have established the following maximal L p -regularity for the fractional model (1.3): for any 1 < p < ∞, u ∈ L p (R + ; D(A)) and under suitable conditions on the operator A (see Theorem 3 in Sect. 2 for details). Further, they applied the theory to analyze nonautonomous and semilinear problems [6,8]. See also [44] for closely related maximal regularity results for Volterra integrodifferential equations. The discrete analogue of (1.4) is important for the numerical analysis of nonautonomous and nonlinear fractional evolution problems. The only existing result we are aware of is the very recent work of Lizama [37]. Specifically, Lizama studied the following fractional difference equation with 0 < α < 1: where u 0 = 0 and Δ α is a certain fractional difference operator. The author established the maximal p -regularity for the problem, under the condition that the set {δ(z)(δ(z)− T ) −1 : |z| = 1, z = 1} is R-bounded, with δ(z) = z 1−α (1 − z) α , following the work of Blunck [10]. It can be interpreted as a time-stepping scheme: upon letting T = τ α A and f n = τ α g n , we get τ −α Δ α u n = Au n + g n . Hence, it amounts to a convolution quadrature generated by the kernel z 1−α (1 − z) α . However, this scheme lacks the maximal p -regularity if A = Δ, the Dirichlet Laplacian operator in bounded domains.
In this work, we address the following question: Under which conditions do the time discretizations of (1.3) preserve the maximal p -regularity, uniformly in the time step size τ ? We provide an analysis for several time-stepping schemes, including the convolution quadratures generated by the implicit Euler method and second-order backward difference formula [11,21], the L1 scheme [36,46], the explicit Euler method [50] and a fractional variant of the Crank-Nicolson method. Amongst them, the convolution quadrature is relatively easy to analyze. In contrast, the L1 scheme and explicit Euler method are easy to implement, but challenging to analyze. The explicit Euler method requires a bounded numerical range of the operator A and the step size τ to be small enough. The maximal p -regularity of the Crank-Nicolson method behaves like the implicit Euler scheme when 0 < α < 1 and like the explicit Euler scheme when 1 < α < 2. Our proof strategy follows closely the recent works [25,28] and employs the (discrete) Fourier multiplier technique of Blunck [10].
The rest of the paper is organized as follows. In Sect. 2 we recall basic tools for showing maximal p -regularity, including R-boundedness, UMD spaces, and Fourier multiplier theorems. Then four classes of time-stepping schemes, i.e., convolution quadrature, L1 scheme, explicit Euler method and a variant of the fractional Crank-Nicolson method, are discussed in Sects. 3-6, respectively. In Sect. 7, we discuss the extension to nonzero initial data. Last, in Sect. 8, we illustrate the results with several concrete examples.
We conclude the introduction with some notation. For a Banach space X , we denote by B(X ) the set of all bounded linear operators from X into itself. For a linear operator A on X , we denote by σ (A) and ρ(A) its spectrum and resolvent set, respectively. We denote the unit circle in the complex plane C by D = {z : |z| = 1}, and D = {z : |z| = 1, z = ±1}. Given any θ ∈ (0, π), the notation Σ θ denotes the open sector Σ θ := {z ∈ C : | arg(z)| < θ, z = 0}, where arg(z) denotes the argument of z ∈ C\{0} in the range (−π, π]. Throughout, the notation c and C, with or without a subscript/superscript, denote a generic constant, which may differ at different occurrences, but it is always independent of the time step size τ and the number N of time steps.

Preliminaries
In this section we collect basic results on the maximal L p -regularity and related concepts, especially R-boundedness, UMD spaces, and Fourier multiplier theorems, used in the fundamental work of Weis [48], where he characterized the maximal L pregularity of an operator A in terms of its resolvent operator R(z; A) := (z − A) −1 . We refer readers to the review [30] for details.

R-boundedness
The concept of R-boundedness plays a crucial role in Weis' operator-valued Fourier multiplier theorem and its discrete analogue. A collection of operators M = {M(λ) : In particular, if ⊂ {z ∈ C : |z| ≤ c 0 } for some c 0 > 0, then the set {λI : λ ∈ } is R-bounded with an R-bound 2c 0 . This fact will be used extensively below.
There are a number of basic properties of R-bounded sets, summarized below. They follow from definition and the proofs can be found in [30].
The following useful result is a slight extension of [10, Corollary 3.5].

Operator-valued multiplier theorems on UMD spaces
Now we recall the concept of UMD spaces, which is essential for multiplier theorems. Let S(R; X ) denote the space of rapidly decreasing X -valued functions. A Banach space X is said to be a UMD space if the Hilbert transform defined on the space S(R; X ), can be extended to a bounded operator on L p (R; X ) for all 1 < p < ∞. Equivalently, this can be characterized by unconditional martingale differences, hence the abbreviation UMD. Examples of UMD spaces include Hilbert spaces, finite-dimensional Banach spaces, and L q (Ω, dμ) ((Ω, μ) is a σ -finite measure space, 1 < q < ∞), and its closed subspaces (e.g., Sobolev spaces W m, p (Ω), 1 < p < ∞), and the product space of UMD spaces. Throughout, X always denotes a UMD space. Next we recall the concept of R-sectoriality operator. The definition below is equivalent to [30,Section 1.11] by changing A to −A and changing θ to π − θ . Similarly, A is said to be R-sectorial of angle θ if (i), (ii) and the following condition hold:

Definition 1 An operator
The following theorem is a simple consequence of Dore [13, Theorem 2.1] and Weis [48,Theorem 4.2].

Theorem 1 A densely defined closed operator A on a UMD space X has maximal parabolic L p -regularity (1.2) if and only if A is R-sectorial of angle π/2.
The "if" direction in Theorem 1 is a consequence of the following operator-valued Fourier multiplier theorem [48,Theorem 3.4], where F denotes the Fourier transform on R, i.e., Theorem 2 Let X be a UMD space. Let M : R\{0} → B(X ) be differentiable such that the set with an R-bound c R . Then M f := F −1 (M(·)(F f )(·)) extends to a bounded operator Further, there exists c p,X > 0 independent of M such that the operator norm of M is bounded by c R c p,X .
Using Theorem 2, one can similarly show the following maximal regularity result for the fractional model (1.3) [6][7][8], which naturally extends the "if" part of Theorem 1 to the fractional case.

Theorem 3 Let A be an R-sectorial operator of angle απ/2 on a UMD space X .
Then the solution of (1.3) satisfies the maximal L p -regularity estimate (1.4) for any 1 < p < ∞.
In this work, we discuss the discrete analogue of Theorem 3 for a number of timestepping schemes for solving (1.3), under the same condition on the operator A, using a discrete version of Theorem 2 due to Blunck [10]. We slightly abuse F for the Fourier transform on Z + := {n ∈ Z : n ≥ 0}, which maps a sequence ( f n ) ∞ n=0 to its Fourier series on the interval (0, 2π), i.e., and let F −1 θ denote the inverse Fourier transform with respect to θ , i.e., The following result is an immediate consequence of [10, Theorem 1.3], and will be used extensively; see also [25] for a proof with a more explicit constant. The statement is equivalent to Blunck's original theorem via the transformation ξ = e −iθ , but avoids introducing a different notation M(θ ).

Theorem 4 Let X be a UMD space, and let M
Further, there exists a c p,X > 0 independent of M such that the operator norm of M is bounded by c R c p,X .
To simplify the notations, for a given sequence {M n } ∞ n=0 of operators on a UMD space X , we define the generating function The operator M is then given by (M f ) n = n j=0 M n− j f j , n = 0, 1, . . .. The generating function satisfies the convolution rule

Convolution quadrature
The convolution quadrature of Lubich (see the review [38] and references therein) presents one versatile framework for developing time-stepping schemes for the model (1.3). One salient feature is that it inherits excellent stability property (of that for ODEs). We shall consider convolution quadrature generated by backward Euler (BE) and second-order backward difference formula (BDF2), whose error analysis has been carried out in [11,21].

BE scheme
We first illustrate basic ideas to prove discrete maximal regularity on BE scheme in time t, with a constant time step size τ > 0. The BE scheme for (1.3) is given by: where the BE approximation∂ α τ u n to ∂ α t u(t n ) is given bȳ where δ(ξ ) = 1 − ξ is the characteristic function of the BE method. Now we can state the discrete maximal regularity of the BE scheme (3.1).
Theorem 5 Let X be a UMD space, 0 < α < 1 or 1 < α < 2, and let A be an R-sectorial operator on X of angle απ/2. Then the BE scheme (3.1) has the following maximal p -regularity where the constant c p,X is independent of N , τ and A, and c R denotes the R-bound of the set of operators {z R(z; A) : z ∈ Σ απ/2 }.
Proof By multiplying both sides of (3.1) by ξ n and summing over n, we have It suffices to compute the term ∞ n=1 ξ n∂ α τ u n . By noting u 0 = 0, the definition of the BE approximation (3.2) and discrete convolution rule (2.5), we deduce Consequently, upon letting f 0 = 0, we arrive at The Rsectoriality of angle απ/2 of the operator A ensures that the operator τ −α δ(ξ ) α − A is invertible. Meanwhile, the generating function of the BE approximation∂ α τ u is given by Appealing to the R-sectoriality of A again gives that z R(z; A) is analytic and R-bounded in the sector Σ απ/2 , which imply that M(ξ ) is differentiable and R-bounded for ξ ∈ D . Direct computation yields which together with Lemma 1 (iii)-(iv) implies the R-boundedness of the set (2.2). Then the desired result follows from Theorem 4.

Remark 1
The BE scheme (3.2) is identical with the Grünwald-Letnikov formula, a popular difference analogue of the Riemann-Liouville fractional derivative ∂ α t u [43], which has been customarily employed for discretizing (1.3).

Second-order BDF scheme
Next we consider the convolution quadrature generated by the second-order backward difference formula (BDF2) for discretizing the model (1.3): where the BDF2 approximation∂ α τ u n to ∂ α t u(t n ), t n = nτ , is given bȳ with the characteristic function δ(ξ ) We approximate the fractional derivative ∂ α t u(t n ) by the BDF2 convolution quadrature (3.4), and consider the scheme (3.3) with zero starting values u 0 = u 1 = 0. Note that for the BDF2 scheme (and other higher-order linear multistep methods), the initial steps have to be corrected properly in order to achieve the desired accuracy [11,21]. The next result gives the discrete maximal regularity of the scheme (3.3).
Theorem 6 Let X be a UMD space, 0 < α < 1 or 1 < α < 2, and let A be an R-sectorial operator on X of angle απ/2. Then the BDF2 scheme (3.3) satisfies the following discrete maximal regularity where the constant c p,X is independent of N , τ and A, and c R denotes the R-bound of the set of operators {z R(z; A) : z ∈ Σ απ/2 }.
Proof In a straightforward manner, upon letting f 0 = f 1 = 0, we obtain The R-sectoriality of the operator A implies the R-boundedness of the set {M(ξ ) : Since d(ξ ) is bounded on D , Lemma 1 (iii)-(iv) and Theorem 4 give the desired assertion.

L1 scheme
Now we discuss one time-stepping scheme of finite difference type for simulating subdiffusion-the L1 scheme [36,46]-which is easy to implement and converges robustly for nonsmooth data, hence very popular. However, unlike convolution quadrature, finite difference type methods are generally challenging to analyze. For the subdiffusion case, i.e., α ∈ (0, 1), it approximates the (Caputo) fractional derivative ∂ α t u(t n ) with a time step size τ by where the weights b j are given by For the case α ∈ (1, 2), the L1 scheme reads [46] where δ t u j− 1 2 = u j − u j−1 denotes central difference, and a j = ( j + 1) 2−α − j 2−α , and we have abused the notation∂ α τ u n for approximating ∂ α t u(t n− 1 2 ). Formally, it can be obtained by applying (4.1) to the first derivative ∂ t u, in view of the identity , and then discretizing the ∂ t u with the Crank-Nicolson type method. The scheme requires ∂ t u(0), in addition to the initial condition u(0). Accordingly, we approximate the right hand side of (1.3) by a Crank-Nicolson type scheme. In sum, the L1 scheme reads Remark 2 For α ∈ (0, 1), Lin and Xu [36] proved that the L1 scheme is uniformly O(τ 2−α ) accurate for C 2 solutions; and for α ∈ (1, 2), Sun and Wu [46] showed that it is uniformly O(τ 3−α ) accurate for C 3 solutions. It is worth noting that even for smooth initial data and source term, the solution of fractional-order PDEs may not be C 2 in general. In fact, the L1 scheme is generally only first-order [20,23].
For the analysis, we recall the polylogarithmic function Li p (z), p ∈ R and z ∈ C, defined by The function Li p (z) is well defined on {z : |z| < 1}, and it can be analytically continued to the split complex plane C\[1, ∞) [14]. With z = 1, it recovers the Riemann zeta function ζ( p) = Li p (1). First we state the solution representation.

4)
with the generating functions Proof We first show the representation for α ∈ (0, 1), and the case α ∈ (1, 2) is analogous. Multiplying both sides of (4.3) by ξ n and summing over n yield Using the polylogarithmic function Li p (z), b(ξ ) is given by from which the desired solution representation follows directly.
We shall need the following result, which is of independent interest.

Lemma 5 For the function δ(ξ ) defined by (4.4), there holds
where d(ξ ) is uniformly bounded on D .
Now we can give the discrete maximal regularity for the L1 scheme (4.1).
Theorem 7 Let X be a UMD space, 0 < α < 1 or 1 < α < 2, and let A be an R-sectorial operator on X of angle απ/2. Then the L1 scheme (4.3) satisfies the following discrete maximal regularity

where the constant c p,X is independent of N , τ and A, and c R denotes the R-bound of the set of operators {z R(z; A) : z ∈ Σ απ/2 }.
Proof First we consider the case 0 < α < 1. Upon setting f 0 = 0, Lemmas 3 and 4 yield where δ(ξ ) is defined by (4.4). By Lemma 4, we have where the latter set is R-bounded by assumption. Meanwhile, where, by Lemma 5, d(ξ ) is uniformly bounded on D . By Lemma 1, the set Thus we deduce from Theorem 4 the desired assertion.
Next we consider the case of 1 < α < 2. In this case, we let g 0 = 0 and g n = f n , n ≥ 1, to obtain In view of the relation δ(ξ ) = The rest of the proof follows like before, using Lemma 5. ∈ (0, 1), the piecewise constant discontinuous Galerkin (PCDG) method in [39] leads to a time-stepping scheme identical with the L1 scheme. The PCDG is given by: find u n such that Next we derive the explicit expression for the discrete approximation∂ α τ u n

Remark 3 For the model (1.3) with α
Thus it is identical with the L1 scheme, and Theorem 7 applies.

Explicit Euler method
Now we analyze the explicit Euler method for discretizing (1.3) in time: where the approximation∂ α τ u n denotes the BE approximation (3.2). A variant of the scheme was analyzed in [50]. By multiplying (5.1) by ξ n and summing up the results for n = 1, 2, . . . , we obtain Recall that the numerical range S(A) of an operator A is defined by [42, pp. 12] Theorem 8 Let X be a UMD space, 0 < α < 1 or 1 < α < 2, and let A be an Rsectorial operator of angle απ/2 such that S(A) ⊂ C\Σ ϕ for some ϕ ∈ (απ/2, π]. Then, under the condition (for small > 0) the scheme (5.1) satisfies the following discrete maximal regularity where the constant c p,X depends only on , ϕ and α (independent of τ and A), and c R denotes the R-bound of the set {z R(z; A) : z ∈ Σ απ/2 }.

Remark 4
The constant in condition (5.3) is sharp. The scaling factor τ α is one notable feature of the model (1.3), and for α ∈ (0, 1), the exponent α agrees with that in the stability condition in [50]. Hence, the smaller the fractional order α is, the smaller the step size τ should be taken.

Remark 5
The condition (5.3) covers bounded operators, e.g., finite element approximations of a self-adjoint second-order elliptic operator. For a self-adjoint discrete approximation, the numerical range S(A) is the closed interval spanned by the largest and smallest eigenvalues, but in general, the numerical range S(A) has to be approximated [19,Section 5.6]. analytic in a neighborhood of the ball B(z 0 , ρ), centered at z 0 with radius ρ, then the set of operators {F(z) : λ ∈ B(z 0 , ρ/2)} is R-bounded on X , and its R-bound is at most

Fractional Crank-Nicolson method
By the fractional Crank-Nicolson method, we mean the following scheme: where the approximation∂ α τ u n denotes the BE approximation (3.2). When α = 1, (6.1) coincides with the standard Crank-Nicolson method. For any 0 < α < 2, one can verify that it is second-order in time, provided that the solution is sufficiently smooth [22]. By multiplying (5.1) by ξ n and summing up the results for n = 1, 2, . . . , we obtain First, we prove the maximal p -regularity for (6.1) in the case 0 < α < 1.

Theorem 9
Let X be a UMD space, 0 < α < 1, and let A be an R-sectorial operator on X of angle απ/2. Then the scheme (6.1) satisfies the following discrete maximal regularity where the constant c p,X depends only on α (independent of τ and A), and c R denotes the R-bound of the set of operators {z R(z; A) : z ∈ Σ απ/2 }.
Proof It suffices to prove that the family of operators where the functions ρ(θ) and ψ(θ) are defined respectively by and It is straightforward to compute Thus α 2 θ − ψ(θ) is an increasing function of θ , taking values from 0 to απ as θ changes from 0 to 2π . Thus τ −α δ(e iθ ) ∈ Σ απ/2 , and by Lemma 1, the set Let the function ψ be defined in (6.3), and θ ϕ ∈ (0, π) be the unique root of the equation Then we have the following result for the case 1 < α < 2.
Theorem 10 Let X be a UMD space, 1 < α < 2, and let A be an R-sectorial operator on X of angle απ/2 such that S(A) ⊂ C\Σ ϕ for some ϕ ∈ (απ/2, π). Then, under the condition (for small > 0) the scheme (6.1) satisfies the following discrete maximal regularity where the constant c p,X depends only on , ϕ and α (independent of τ and A), and c R denotes the R-bound of the set {z R(z; A) : z ∈ Σ απ/2 }.

Inhomogeneous initial condition
In this section, we consider maximal p -regularity for the problem with nontrivial initial conditions: We focus on the BE scheme since other schemes can be analyzed similarly. For (7.2), the BE scheme reads [21,22]: with u 0 = v, find u n such that where∂ α τ denotes the BE convolution quadrature (3.2). We shall need the scaled L p -norm and weak L p -norm (cf. [9, section 1.3]) The main result of this section is the following theorem.

Theorem 11
Let X be a Banach space, 0 < α < 1, and let A be a sectorial operator on X of angle απ/2. Then the BE scheme (7.3) has the following maximal p -regularity where the constant c p depends on the bound of the set of operators {z R(z; A) : z ∈ Σ απ/2 }, independent of N , τ and A.
Proof By multiplying both sides of (7.3) by ξ n and summing over n, we have Au n ξ n = 0.

Remark 6
The proof shows that in the absence of a source term f , the maximal pregularity of (7.3) only requires the sectorial property of A, rather than the R-sectorial property. The general case (with nonzero source and nonzero initial data) is a linear combination of (1.3) and (7.2).

Remark 7
We have focused our discussions on the Caputo fractional derivative, since it allows specifying the initial condition as usual, and thus is very popular among practitioners. In the Riemann-Liouville case, generally it requires integral type initial condition(s) [27], for which the physical interpretation seems unclear.
In the proof of Theorem 11, we first prove two end-point cases, p = 1/α and p = ∞. Then we use real interpolation method for the case 1/α < p < ∞. The real interpolation method also holds for 0 < p < 1 ([9, Theorem 5.2.1]). Thus, we have the following theorem in the case 1 < α < 2. The proof is omitted, since it is almost identical with the proof of Theorem 11. Theorem 12 Let X be a Banach space, 1 < α < 2, and let A be a sectorial operator on X of angle απ/2. Then the BE scheme (7.3) has the following maximal p -regularity: where the constant c p depends on the bound of the set of operators {z R(z; A) : z ∈ Σ απ/2 }, independent of N , τ and A.

Examples and application to error estimates
In this section, we present a few examples of fractional evolution equations which possess the maximal L p -regularity, and investigate conditions under which the timestepping schemes in Sects. 3-6 satisfy the maximal p -regularity.
Example 1 (Continuous problem) Consider the following time fractional parabolic equation in a bounded smooth domain Ω ⊂ R d (d ≥ 1): where T > 0 is given and Δ denotes the Laplacian operator. In the appendix, we show that the L q realization Δ q in X = L q (Ω) of Δ is an R-sectorial operator in X with angle θ ∈ (0, π), and that Δ q v coincides with Δv in the domain D(Δ q ) of Δ q : Thus Theorem 3 implies that the solution u q of satisfies u q (·, t) ∈ D(Δ q ) for almost all t ∈ R + and ∂ α t u q L p (0,T ;L q (Ω)) + Δ q u q L p (0,T ;L q (Ω)) ≤ ∂ α t u q L p (R + ;L q (Ω)) + Δ q u q L p (R + ;L q (Ω)) ≤ c p,X f L p (R + ;L q (Ω)) = c p,X f L p (0,T ;L q (Ω)) , ∀ 1 < p, q < ∞.
As an application of the maximal p -regularity, we present error estimates for the numerical solutions by the BE scheme (3.1), with the scaled L p -norm (7.4). Other time-stepping schemes can be analyzed similarly.
When q > d, error estimates in such strong norms as W 2,q (Ω) → W 1,∞ (Ω) can be used to control some strong nonlinear terms in the numerical analysis of nonlinear parabolic problems [1,2,17]. We will explore such an analysis in the future.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
In particular, we have The Gaussian estimate (A.2) yields That is, E q (t/2)v ∈ L 2 (Ω) for v ∈ L q (Ω) and t > 0. Hence, (A.1) implies where the last inclusion is due to the analyticity of the semigroup E 2 (t) [ Since lim t↓0 E q (t)v − v L q (Ω) = 0, the last identity implies That is, Δ q v coincides with the distributional partial derivative Δv in the sense of distributions, i.e., Hence, we have the characterization D(Δ q ) = {Δ −1 v : v ∈ L q (Ω)}.