Regularity theory for time-fractional advection–diffusion–reaction equations

Regularity theory for time-fractional advection-diffusion-reaction equations ✩ arXiv:1902.00850v1 [math.AP] 3 Feb 2019 William McLean School of Mathematics and Statistics, The University of New South Wales, Sydney 2052, Australia Kassem Mustapha, Raed Ali Department of Mathematics and Statistics, KFUPM, Dhahran, 31261, KSA Omar M. Knio Computer, Electrical, Mathermatical Sciences and Engineering Division, KAUST, Thuwal 23955, KSA Abstract We investigate the behavior of the time derivatives of the solution to a linear timefractional, advection-diffusion-reaction equation, allowing space- and time-dependent coefficients as well as initial data that may have low regularity. Our focus is on proving estimates that are needed for the error analysis of numerical methods. The nonlocal nature of the fractional derivative creates substantial difficulties compared with the case of a classical parabolic PDE. In our analysis, we rely on novel energy methods in combination with a fractional Gronwall inequality and certain properties of fractional integrals. Keywords: Fractional PDE, regularity analysis, energy arguments, fractional Gronwall inequality 2010 MSC: 26A33; 35B45, 35B65, 35D30, 35K57, 35Q84, 35R11. 1. Introduction This paper is the sequel to a study [1] of existence and uniqueness of the weak solution to a time-fractional PDE of the form
~ u + (a∂t1−α + b)u = g ∂t u − ∇ · κ ∇∂t1−α − ~F ∂t1−α − G (1) for x ∈ Ω and 0 < t ≤ T , subject to the boundary and initial conditions u(x,t) = 0 for x ∈ ∂ Ω and 0 ≤ t ≤ T , (2) u(x, 0) = u0 (x) for x ∈ Ω. (3) ✩ The authors thank the University of New South Wales (Faculty Research Grant “Efficient numerical simulation of anomalous transport phenomena”), the King Fahd University of Petroleum and Minerals (project No. KAUST005) and the King Abdullah University of Science and Technology. Preprint submitted to Computers and Mathematics with Applications February 6, 2019 Various special cases of this problem occur in descriptions of subdiffusive transport, with the parameter α arising from a continuous-time, random walk model [2, section 3.4] in which the waiting-time distribution is a power law decaying like t −1−α as t → ∞ [3, 4, 5, 6, 7, 8]. For more details, see our related paper [1]. Our purpose here is to derive estimates for the derivatives of u, motivated by their crucial role in the error analysis of numerical methods [9, 10, 11, 12, 13, 14, 15, 16] for applications included in the class (1) of time-fractional problems. For the basic ~ = ~0 and a = b = 0, the fractional diffusion equation, given by the special case ~F = G solution admits a series representation via separation of variables, which, in combination with the asymptotics of the Mittag–Leffler function, yields bounds on the time derivatives of u in various spatial norms [17, 18]. One may also represent the solution in terms of a fractional resolvent [19, 20, 21]. These simple approaches no longer work in the general case, and the analysis that follows relies instead on the tools used in our study [1] of well-posedness: energy methods and a fractional Gronwall inequality. We assume that 0 < α < 1 and that the spatial domain Ω ⊆ Rd (d ≥ 1) is bounded ~ a and b, as well as the source term g, may depend and Lipschitz. The coefficients ~F, G, on x and t, but the generalized diffusivity κ = κ (x) may depend only on x. Our theory requires that for appropriate m ≥ 1,

~ ∈ Cm+1 [0, T ];W∞1 (Ω)d ~ (4) and a, b ∈ Cm [0, T ]; L∞ (Ω) , F, G where Wpk (Ω) denotes the Sobolev space of functions with all partial derivatives up to and including order k belonging to L p (Ω). The generalized diffusivity is permitted to be a bounded, d × d matrix-valued function, that is, κ ∈ L∞ (Ω; Rd×d ). In addition, we ensure that the spatial operator v 7→ −∇ · (κ ∇v) is uniformly elliptic on Ω by assuming κ (x) is symmetric and positive-definite with its minimal eigenvalue bounded away from zero. The fractional time derivative is understood in the Riemann–Liouville sense, that is ∂t1−α v(x,t) = ∂t I α v(x,t) with I α the fractional integral given by I α v(x,t) = Z t 0 ωα (t − s)v(x, s) ds where ωα (t) = t α −1 /Γ(α ). Let h·, ·i denote the inner product in L2 (Ω) or L2 (Ω)d . The weak solution u of (1) is defined by the condition hu(t), vi + (κ I α ∇u)(t), ∇v − (~B1 u)(t), ∇v + h(B2 u)(t), vi = h f (t), vi, (5) for all v ∈ H01 (Ω), where f (t) = u0 + I 1 g(t) and ~ φ )(t), ~B1 φ (t) = I 1 (~F ∂t1−α φ )(t) + I 1 (G B2 φ (t) = I 1 (a∂t1−α φ )(t) + I 1 (b φ )(t). (6) To see why, take the inner product of (1) with v, apply the first Green identity and integrate in time (making use of the initial condition). We proved in our previous paper [1, Theorems 4.1 and 4.2] that the above problem is well-posed in the following sense. 2 Theorem 1. Assume that the coefficients satisfy (4) for m = 1, that the source term g satisfies kg(t)k ≤ Mt η −1 for 0 < t ≤ T , where M and η are positive constants, and that the initial data
u0 ∈ L2 (Ω). Then, problem (1)–(3) has a unique weak solution u ∈ L2 (0, T ); L2 (Ω) that satisfies (5), and is such that 1. 2. 3. 4. 5. The restriction u : (0, T ] → L2 (Ω) is continuous.
If 0 < t ≤ T , then u(t) ∈ H01 (Ω) with ku(t)k + t α /2k∇u(t)k ≤ C ku0 k + Mt η . I α u, B2 u ∈ C([0, T ]; L2 (Ω)) and I α ∇u, ~B1 u ∈ C([0, T ]; L2 (Ω)d ). If t = 0, then I α u = B2 u = 0, I α ∇u = ~B1 u = 0 and u(0) = u0 . If t → 0, then hu(t), vi → hu(0), vi for each v ∈ L2 (Ω). Regarding part 1 and the weak continuity in part 5, we will show in Theorem 12 that u is continuous on the closed interval [0, T ] provided u0 ∈ Ḣ µ (Ω) for some µ > 0. Some partial results on the regularity of u are known. If the coefficients in (1) are ~ = ~0 and b = 0, then by applying I 1−α to both sides the independent of time, and if G fractional PDE may be written in the alternative form ∂tα u − ∇ · (κ ∇u − ~Fu) + au = I 1−α g. (7) Sakamoto and Yamamoto [22, Corollary 2.6] show, for example, that if g = 0 then the solution of (7) satisfies a bound of the form k∂tm uk ≤ Ct −m ku0 k, where k · k denotes the norm in L2 (Ω). Mu, Ahmad and Huang [23] obtain analogous estimates using ~ = ~0 and weighted Hölder norms. Recently, Le et al. [24] studied (1) for the case G ~ ~ a = b = 0, with F = F(x,t). One of their regularity results [24, Theorem 7.3] gives the bound k∂tm uk ≤ Ct −m+1/2 ku0kH 2 (Ω) when g = 0, subject to the restriction 1/2 < α < 1. The next section gathers together some technical preliminaries needed for our analysis, which uses delicate energy arguments, a fractional Gronwall inequality and several properties of fractional integrals to prove a priori estimates for the weak solution u of (1)–(3). In Section 3, we estimate the derivatives of u and ∇u with respect to time assuming u0 ∈ L2 (Ω). For example, Corollary 10 shows that if g(t) ≡ 0 then, with m ≥ 1 such that (4) holds, k∂tm uk + t α /2k∂tm ∇uk + t −α k∂tm−α uk + t −α /2k∂tm−α ∇uk ≤ Ct −m ku0 k for 0 < t ≤ T . Unlike a classical parabolic PDE, the fractional problem (1) exhibits only limited spatial smoothing as t increases [17], and in Section 4 we investigate the consequences of more regular initial data. For example, Theorems 12 and 13 show that when g(t) ≡ 0 and u0 ∈ Ḣ µ (Ω) for 0 ≤ µ ≤ 2, and under additional assumptions on κ and Ω, k∂tm uk + t −α k∂tm−α uk + t α k∂tm ukH 2 (Ω) + t −α /2k∂tm−α ∇uk ≤ Ct −m+α µ /2 ku0 k µ . The paper concludes with an Appendix containing three technical lemmas. 2. Preliminaries and notations This section introduces some notations and states some technical results that will be used in our subsequent regularity analysis. As in our recent paper [1], we define the 3 µ µ quadratic operators Q1 and Q2 , for µ ≥ 0 and 0 ≤ t ≤ T , by µ Q1 (φ ,t) = Z t hφ , I µ φ i ds 0 µ and Q2 (φ ,t) = Z t kI µ φ k2 ds. (8) 0 These operators coincide when µ = 0, so we write Q 0 = Q10 = Q20 . We recall the following positivity property [25, Theorem 2] µ Q1 (φ , T ) ≥ 0 for 0 ≤ µ ≤ 1. (9) µ µ The next four lemmas establish key inequalities satisfied by Q1 and Q2 . Lemma 2 ([13, Lemma 3.2]). If 0 < α < 1 and ε > 0, then Z t Q1α (φ ,t) + ε Q1α (ψ ,t), 4ε (1 − α )2 2t α Q α (φ ,t), Q2α (φ ,t) ≤ 1−α 1 Q1α (φ ,t) ≤ 2t α Q 0 (φ ,t). hφ , I α ψ i ds ≤ 0 (10) (11) (12)
Lemma 3 ([1, Lemmas 2.2 and 2.3]). If 0 < α ≤ 1, then for φ ∈ L2 (0,t), L2 (Ω) , Q2α (φ ,t) ≤ 2 Z t 0
ωα (t − s)Q1α (φ , s) ds and I 1−α kI α φ k2 ≤ 2Q1α (φ ,t).
Furthermore, if φ ∈ W11 (0,t); L2 (Ω) and φ (0) = I α φ ′ (0) = 0, then kφ (t)k2 ≤ 2ω2−α (t) Q1α (φ ′ ,t). µ Lemma 4 ([13, Lemma 3.1]). If 0 ≤ µ ≤ ν ≤ 1, then Q2ν (φ ,t) ≤ 2t 2(ν −µ ) Q2 (φ ,t). We will make essential use of the following fractional Gronwall inequality. Lemma 5 ([26, Theorem 3.1]). Let β > 0 and T > 0. Assume that a and b are nonnegative and non-decreasing functions on the interval [0, T ]. If q : [0, T ] → R is an integrable function satisfying 0 ≤ q(t) ≤ a(t) + b(t)I β q(t) for 0 ≤ t ≤ T , then q(t) ≤ a(t)Eβ b(t)t β
for 0 ≤ t ≤ T . Let M j denote the operator of pointwise multiplication by t j , that is, (M j φ )(t) = t j φ (t), and note the commutator properties (for any integer j ≥ 1 and any real µ ≥ 0) ∂t j M − M ∂t j = j∂t j−1 , ∂t M j − M j ∂t = jM j−1 , The following identities then follow by induction on m. 4 M I µ − I µ M = µ I µ +1 . (13) m,q Lemma 6. For 0 ≤ q ≤ m and µ ≥ 0, there exist constant coefficients am,q j , bj , µ m, cm,µ and d j such that q m− j q− j ∂tq M m = M m ∂tq + ∑ am,q ∂t , j M (14) q− j M m ∂tq = ∂tq M m + ∑ bm,q M m− j , j ∂t (15) j=1 q j=1 m m, µ I µ M m = M m I µ + ∑ c j M m− j I µ + j , j=1 m m, µ M mI µ = I µ M m + ∑ d j I µ + j M m− j . (16) (17) j=1 m, µ m,q For later reference, we set am,q 0 = b 0 = c0 m,q ãm,q j = aq− j , m, µ m,q b̃m,q j = bq− j , c̃ j m, µ = d0 = 1 and m, µ = cm− j , m, µ m, µ d˜j = dm− j . When µ = 0 the formulas involving I µ become redundant, and we see that cm,0 j = 0 = d m,0 for 1 ≤ j ≤ m. Likewise, ω− j (t) = 0 for 0 ≤ j ≤ m since Γ(z) has a pole j at z = − j. We conclude this section by noting that if m ≥ 1 and µ ≥ 0, then ∂tm I µ φ (t) = I µ ∂tm φ (t) + m−1 ∑ (∂t j φ )(0)ωµ − j (t) j=0
for φ ∈ W1m (0,t); L2 (Ω) , (18) which amounts to a restatement of the relation between the Riemann–Liouville and Caputo fractional derivatives. 3. Regularity of the weak solution Our aim in this section is to estimate higher-order time derivatives of u assuming appropriate bounds on the higher-order time derivatives of f (and hence, ultimately, of g), as well as sufficient smoothness of the coefficients in (1). We will not attempt to prove the existence of the higher-order derivatives of u, which could be done by estimating the corresponding derivatives of the projected solution uX from our earlier paper [1] corresponding to a finite dimensional subspace X = Xn ⊂ H01 (Ω), and then taking appropriate limits as n → ∞. For the remainder of the paper, we assume that (4) holds, and that kg( j−1)(t)k = O(t α − j ) as t → 0, for 1 ≤ j ≤ m. It follows that the existence and uniqueness of the weak solution u are guaranteed by Theorem 1. Henceforth, C will denote a generic constant that may depend on the coefficients in (1), the spatial domain Ω, the time interval [0, T ], the fractional exponent α , 5 the parameter η , and the integer m in (4). Also, we rescale the time variable if necessary so that the minimum eigenvalue of κ satisfies
λmin κ (x) ≥ 1 for x ∈ Ω. (19) For brevity, we introduce some more notations. Let µ (Bψ φ )(t) = ψ (t) I µ φ (t) − I 1 (ψ ′ I µ φ )(t), for 0 ≤ µ ≤ 1. (20) Integrating by parts and recalling (6), we find that ~B1 = Bα~ + B1~ F G and B2 = Bαa + B1b . (21) Generalizing (20), for j ∈ {0, 1, 2, . . .} we put   µ, j 1− µ µ Bψ φ (t) = ∂t j M j I 1 (ψ∂t φ ) (t) = (M j Bψ φ )( j) (t), and generalizing (8) we put
µ, j µ Qi (φ ,t) = Qi (M j φ )( j) ,t , for 0 ≤ t ≤ T and i ∈ {1, 2}, with Q 0, j = Q10, j = Q20, j . The next result relies on Lemma 15 from the Appendix. Lemma 7. For 0 < t ≤ T and for m ≥ 1, Q1α ,m (u,t) + Q2α ,m (∇u,t) ≤ Ct α m ∑ Q0, j ( f ,t), j=0 and m Q 0,m (u,t) + Q1α ,m (∇u,t) ≤ C ∑ Q 0, j ( f ,t). j=0 Proof. Since (I α ∇u)(0) = 0 by part 4 of Theorem 1,  Zt  Z t 1−α α ′ κ ∇∂s u(s), ∇v ds = κ (I ∇u) (s) ds, ∇v = κ (I α ∇u)(t), ∇v , 0 0 and by the identity in (17), M m I α ∇u = I α M m ∇u + m−1 ∑ d˜m,j α I α +m− j M j ∇u. j=0 Thus, multiplying both sides of (5) by t m yields m α α +m− j M j ∇u), ∇vi hM m u, vi + hκ I α M m ∇u, ∇vi + ∑ d˜m, j hκ I j=1 = m 1 m m α hM m Bα~F u + M m B1G ~ u, ∇vi − hM Ba u + M Bb u, vi + hM f , vi 6 for v ∈ H01 (Ω). We have ∂tm I α +m− j M j ∇u = ∂t j ∂tm− j I m− j I α M j ∇u = ∂t j I α M j ∇u = I α ∂t j M j ∇u, where the final step follows by (18) because ∂ti (M j ∇u)(0) = 0 for 0 ≤ i ≤ j − 1 ≤ m − 1. Likewise, ∂tm I α M m ∇u = I α ∂tm M m ∇u because ∂t j (M m ∇u)(0) = 0 for 0 ≤ j ≤ m − 1, and therefore h∂tm M m u, vi + hκ I α ∂tm M m ∇u, ∇vi = hB~Fα ,m u + B1,m ~ u, ∇vi G − hBaα ,m u + Bb1,mu, vi − m−1 ∑ d˜m,j α hκ I α ∂t j M j ∇u, ∇vi + h∂tm M m f , vi. (22) j=0 1,m 2 2 We let E (u) = 2kB~Fα ,m uk2 + kBaα ,m uk2 + 2kB1,m ~ uk + kBb uk , and conclude using the G Cauchy–Schwarz inequality that h∂tm M m u, vi + hκ I α ∂tm M m ∇u, ∇vi ≤ E (u) + C m−1 ∑ kI α ∂t j M j ∇uk2 j=0 + 12 k∇vk2 + 21 kvk2 + h∂tm M m f , vi. Choosing v = I α ∂tm M m u, integrating over the time interval (0,t) and using (19), we have Q1α ,m (u,t) + 12 Q2α ,m (∇u,t) ≤ Z t m−1 E (u) ds + C 0 ∑ Q2α , j (∇u,t) j=0 + Q2α ,m (u,t) + Zt 0 h∂tm M m f , I α ∂tm M m ui ds, and by the Cauchy-Schwarz inequality and Lemma 3,
I 1−α kI α (∂tm (M m u))k2 (t) ≤ 2Q1α (∂tm (M m u),t). Thus, Z t Z t k∂tm M m f k kI α ∂tm M m uk ds 0 1/2 Z t 1/2 Z t (t − s)−α kI α ∂tm M m uk2 ds ≤C (t − s)α k∂tm M m f k2 ds 0 0 0 h∂tm M m f , I α ∂tm M m ui ds ≤ 1/2  1/2  I 1−α (kI α (∂tm M m u)k2 )(t) ≤ C t α Q 0,m ( f ,t) ≤ Ct α Q 0,m ( f ,t) + 12 Q1α ,m (u,t), 7 implying that the function qm (t) = Q1α ,m (u,t) + Q2α ,m (∇u,t) satisfies qm (t) ≤ 2 Z t m−1 ∑ Q2α , j (∇u,t) + Ct α Q0,m ( f ,t). E (u) ds + C 0 j=0 By Lemma 15, m Q 0 (B~α ,m u,t) + Q 0(Baα ,m u,t) ≤ C ∑ Q2α , j (u,t) F (23) j=0 and, applying Lemma 4 with ν = 1 and µ = α , m m j=0 j=0 1, j 2(1−α ) 0 1,m Q 0 (B1,m ∑ Q2α , j (u,t) ~ u,t) + Q (Bb u,t) ≤ C ∑ Q2 (u,t) ≤ Ct G (24) for m ≥ 0. By combining the above estimates, qm (t) ≤ CQ2α ,m (u,t) + C m−1 ∑ q j (t) + Ct α Q0,m ( f ,t). j=0 Consequently, we conclude (recursively) that m qm (t) ≤ C ∑ Q2α , j (u,t) + Ct α j=0 m ∑ Q0, j ( f ,t), j=0 so, by applying the first inequality in Lemma 3 with φ = (M j u)( j) , qm (t) ≤ Ct α m m ∑ Q0, j ( f ,t) + C ∑ Z t j=0 0 j=0 ωα (t − s)q j (s) ds. Therefore, a repeated application of Lemma 5 yields the first desired estimate. To show the second estimate, choose v = ∂tm M m u in (22) and obtain k∂tm M m uk2 + hκ I α ∂tm M m ∇u, ∂tm M m ∇ui = −hEu, ∂tm M m ui m−1 − ∑ d˜m,j α hκ I α ∂t j M j ∇u, ∂tm M m ∇ui + h∂tm M m f , ∂tm M m ui, j=0 1,m where Eu = ∇ · B~Fα ,m u + Baα ,m u + ∇ · B1,m ~ u + Bb u. The first and the last terms on the G right-hand side are bounded by kEuk2 + k∂tm M m f k2 + 21 k∂tm M m uk2 so, after integrating in time, using (19) and applying (10) (for a sufficiently large ε ), 1 0,m (u,t) + Q1α ,m(∇u,t) 2Q ≤ Z t 0 kEu(s)k2 ds + Q 0,m( f ,t) + 21 Q1α ,m (∇u,t) m−1 +C ∑ Q1α , j (∇u,t). j=0 8 Since ∇ · (~ F ∂t1−α u) = (∇ · ~F)∂t1−α u + ~F · ∇∂t1−α u, we see that
α ,m α ,m u + B~F· ∇u, ∇ · B~Fα ,m u = ∂tm M m I 1 ∇ · (~F ∂t1−α u) = B∇· ~ F and therefore, applying Lemma 15 followed by Lemma 4, Z t 0  kEu(s)k2 ds ≤ 4 Q 0 (∇ · B~Fα ,m u,t) + Q 0(Baα ,m u,t)  0 1,m + Q 0 (∇ · B1,m u,t) + Q (B u,t) b ~ G  m  α, j α, j ≤ C ∑ Q2 (u,t) + Q2 (∇u,t) . j=0 Hence, the function qm (t) = Q 0,m (u,t) + Q1α ,m (∇u,t) satisfies qm (t) ≤ 2Q 0,m ( f ,t) + C m−1 m j=0 j=0 ∑ Q1α , j (∇u,t) + C ∑   Q2α , j (u,t) + Q2α , j (∇u,t) , and so, using (11) and (12), it follows that qm (t) ≤ 2Q 0,m ( f ,t) + C m−1 ∑ q j (t) + C j=0   Q2α ,m (u,t) + Q2α ,m (∇u,t) . By the first inequality in Lemma 3 and (12), Q2α ,m (u,t) + Q2α ,m (∇u,t) ≤ C Z t 0 ωα (t − s)qm (s) ds, and thus by Lemma 5, qm (t) ≤ CQ 0,m ( f ,t) + C m−1 ∑ q j (t). j=0 Applying this inequality recursively gives m qm (t) ≤ C ∑ Q 0, j ( f ,t), j=0 which completes the proof.  We can now show pointwise bounds for the norms in L2 (Ω) of the time derivatives of u and ∇u. Theorem 8. For m ≥ 1 and 0 < t ≤ T , k(∂tm u)(t)k2 + t α k(∂tm ∇u)(t)k2 ≤ Ct −1−2m m+1 ∑ Q0, j ( f ,t). j=0 9 m,m−1 Proof. Since M ∂tm = ∂tm M − m∂tm−1 , we see using (15) (and setting b̃m = 0) that M m+1 ∂tm = M m ∂tm M − mM m ∂tm−1 = m+1 ∑ j=1 and hence k(M m+1 ∂tm u)(t)k2 ≤ C m,m−1
j−1 b̃m,m ∂t M j , j−1 − mb̃ j−1 m+1 ∑ k(∂t j−1M j u)(t)k2. (25) j=1 Using the second inequality in Lemma 3 with φ = ∂t j−1 M j u and the first bound in Lemma 7, we get k(∂t j−1 j j M j u)(t)k2 ≤ Ct 1−α Q1α (∂t M j u,t) ≤ Ct ∑ Q 0,ℓ ( f ,t) ℓ=0 and so k(∂tm u)(t)k2 = t −2m−2 k(M m+1 ∂tm u)(t)k2 ≤ Ct −1−2m m+1 ∑ Q0, j ( f ,t). j=0 Applying the same argument to ∇u in place of u, and using the second bound in Lemma 7, the result follows.  Next, we estimate fractional time derivatives of u and ∇u. These bounds will later help in our study of spatial regularity, and reflect the presence of the fractional time derivative in (1). Theorem 9. For m ≥ 1 and 0 < t ≤ T , k(∂tm−α u)(t)k2 + t α k(∂tm−α ∇u)(t)k2 ≤ Ct −1−2(m−α ) m+1 ∑ Q0, j ( f ,t). j=0 Proof. Using the inequality (25), (M m ∂tm−α u)(t) 2 = (M m ∂tm−1 ∂t1−α u)(t) 2 m ≤ C ∑ (∂t j−1 2 M j ∂t1−α u)(t) , j=1 (26) and using (13) and (18) (with m = 1), M ∂t1−α u = M ∂t I α u = M I α ∂t u + u(0)ωα Thus, by (17),
= (I α M + α I α +1 )∂t u + u(0)M ωα
= I α M u′ + α I α u − u(0) + α u(0)ω1+α = I α (M u′ + α u). M j ∂t1−α u = M j−1 I α (M u′ + α u) = j−1 ∑ d˜ℓj−1,α I α + j−1−ℓM ℓ (M u′ + α u). ℓ=0 10 We have
∂t j−1 I α + j−1−ℓM ℓ (M u′ + α u) = ∂tℓ ∂t j−1−ℓ I j−1−ℓ I α M ℓ (M u′ + α u) = ∂tℓ I α M ℓ (M u′ + α u) = I α ∂tℓ M ℓ (M u′ + α u), where we used the identity (18) (with m = 1) and the fact that ∂ti M ℓ (M u′ + α u)(0) = 0 for 0 ≤ i ≤ ℓ − 1. Hence, (∂t j−1 where φℓ = j−1 j−1 ℓ=0 ℓ=0 ∑ d˜ℓj−1,α I α ∂tℓ M ℓ (M u′ + α u)(t) ≤ C ∑ kI α φℓ (t)k M j ∂t1−α u)(t) = ∂tℓ M ℓ (M u′ + α u). Using (14), ℓ ℓ ℓ,ℓ i i i ′ i ′ i−1 ′ u + α∂ti u φℓ = ∑ ãℓ,ℓ i M ∂t (M u + α u) = ∑ ãi M M ∂t u + i∂t i=0 i=0 and so, by Theorem 8, ℓ+1 ℓ+1 r+1 ℓ+2 r=0 r=0 i=0 r=0
kφℓ (t)k2 ≤ C ∑ k(M r ∂tr u)(t)k2 ≤ C ∑ t −1 ∑ Q 0,i ( f ,t) ≤ Ct −1 ∑ Q 0,r ( f ,t). (27) q 0,r Since kφℓ (t)k ≤ Cω1/2 (t)ψℓ (t) where ψℓ (t) = ∑ℓ+2 r=0 Q ( f ,t) is nondecreasing, we α see that kI φℓ (t)k ≤ Cωα +1/2 (t)ψℓ (t). Therefore, (∂t j−1 M j ∂t1−α u)(t) 2 j−1 j−1 ≤ C ∑ kI α φℓ (t)k2 ≤ C ∑ (t (α +1/2)−1)2 ψℓ (t)2 ℓ=0 ℓ=0 j+1 ≤ Ct 2α −1 ∑ Q 0,ℓ ( f ,t), ℓ=0 k(∂tm−α u)(t)k2 and the desired bound for follows at once from (26). Replacing u with ∇u in the preceding argument, we have (M m ∂tm−α ∇u)(t) 2 m ≤ C ∑ (∂t j−1 M j ∂t1−α ∇u)(t) 2 j=1 m j−1 ≤C ∑ ∑ kI α φℓ (t)k2 j=1 ℓ=0 where, this time, φℓ = ∂tℓ M ℓ (M ∇u′ + α ∇u) and hence kφℓ (t)k ≤ Cω(1−α )/2 (t)ψℓ (t). It follows that kI α φℓ k ≤ Cω(1+α )/2 (t)ψℓ (t) and therefore t α (M m ∂tm−α ∇u)(t) bounded by Ct α m−1 ∑ ℓ=0 kI α φℓ (t)k2 ≤ Ct α m−1 m+1 ℓ=0 ℓ=0 ∑ (t (1+α )/2−1)2 ψℓ (t)2 ≤ Ct 2α −1 is ∑ Q0,ℓ ( f ,t), as required. The following simplified bounds are perhaps more immediately useful. 11 2  Corollary 10. Let m ≥ 1 and suppose that g : (0, T ] → L2 (Ω) is Cm with kg( j) (t)k ≤ Mt η −1− j for 0 ≤ j ≤ m and some η > 0. (28) Then k(∂tm u)(t)k + t α /2k(∂tm ∇u)(t)k ≤ Ct −m (ku0 k + Mt η ) and k(∂tm−α u)(t)k + t α /2k(∂tm−α ∇u)(t)k ≤ Ct α −m (ku0 k + Mt η ). Proof. Since k f ( j) (t)k = kg( j−1) (t)k ≤ Mt η − j for 1 ≤ j ≤ m + 1, (14) implies that k(M j f )( j) (t)k ≤ CMt η for 1 ≤ j ≤ m + 1, with k f (t)k ≤ ku0 k + M η −1t η . Thus, Q 0, j ( f ,t) ≤ CM 2 t 2η +1 for 1 ≤ j ≤ m + 1, with Q 0 ( f ,t) ≤ Ct(ku0 k + Mt η )2 , so t −1−2m m+1 ∑ Q0, j ( f ,t) ≤ Ct −2m (ku0k + Mt η )2 j=0 and the result follows from Theorems 8 and 9.  4. More regular initial data We will now investigate further the relation between the regularity of u and that of the initial data u0 . In particular, Theorem 13 below extends Corollary 10 and proves a bound used in an error analysis of a finite element discretization of the fractional Fokker–Planck equation [13]. The fractional PDE (1) can be rewritten as u′ − ∇ · (κ∂t1−α ∇u) = h for x ∈ Ω and 0 < t < T ,
~ − (a∂t1−α u + bu). We can therefore apply known where h = g − ∇ · ~ F ∂t1−α u + Gu results for the fractional diffusion equation to establish the following bounds in the norm kvk µ = kAµ /2 vk of the fractional Sobolev space Ḣ µ (Ω), where Aµ /2 is defined via the spectral representation of Av = −∇ · (κ ∇v) using the Dirichlet eigenfunctions on Ω [17, 27]. The results of this section require H 2 -regularity for the Poisson problem, and to ensure this property we make the additional assumptions [28, Theorems 2.2.2.3 and 3.2.1.2] κ is Lipschitz on Ω and Ω is C1,1 or convex. (29) It follows that Ḣ 1 (Ω) = H01 (Ω) and Ḣ 2 (Ω) = H 2 (Ω) ∩ H01 (Ω). We also require that g satisfies (28). Our first result does not assume any additional smoothness of u0 . Theorem 11. Assume (28) and (29). If u0 ∈ L2 (Ω), then t m ku(m) (t)k µ ≤ Cku0 kt −µα /2 + CMt η −µα /2 12 for 0 ≤ µ ≤ 2 and 0 < t ≤ T . Proof. We have [17, Theorems 4.1 and 4.2, and the inequality stated after Theorem 5.4] m Z t m (m) − µα /2 t ku (t)kµ ≤ Ct ku0 k + C ∑ (t − s)−µα /2 s j kh( j) (s)k ds j=0 0 for m ≥ 0 and for 0 ≤ µ ≤ 2, with j  kh( j)(s)k ≤ k∂t j g(s)k + C ∑ k∂tℓ+1−α ∇u(s)k + k∂tℓ∇u(s)k ℓ=0  + k∂tℓ+1−α u(s)k + k∂tℓu(s)k . Corollary 10 shows that kh( j) (s)k is bounded by j   Msη −1− j + C ∑ sα /2−ℓ−1 + s−α /2−ℓ + sα −ℓ−1 + s−ℓ (ku0 k + Msη ) ℓ=0 so s j kh( j)(s)k ≤ Cku0 ksα /2−1 + Msη −1 and hence Z t 0 (t − s)−µα /2 s j kh( j) (s)k ds ≤ Cku0k(ω1−µα /2 ∗ ωα /2)(t) + CM(ω1−µα /2 ∗ ωη )(t)
≤ C ku0kt (1−µ )α /2 + Mt η −µα /2 , completing the proof.  Recall from part 5 of Theorem 1 that for any u0 ∈ L2 (Ω) the solution u(t) converges weakly to u(0) = u0 in L2 (Ω) as t → 0. The first estimate in our next result shows that u(t) → u0 in the norm of L2 (Ω) if we impose some additional spatial regularity on the initial data, namely if u0 ∈ Ḣ µ (Ω) for some µ > 0. The second and third estimates extend the results of Corollary 10. Theorem 12. Assume (28) and (29). If 0 ≤ µ ≤ 2 and u0 ∈ Ḣ µ (Ω), then ku(t) − u0k + t α /2 ∇ u(t) − u0 and, for m ≥ 1,
≤ Cku0 kµ t α µ /2 + Mt η , ku(m) (t)k + t α /2k∇u(m) (t)k ≤ Ct −m ku0 kµ t α µ /2 + Mt η with

k∂tm−α u(t)k + t α /2k(∂tm−α ∇u)(t)k ≤ Ct α −m ku0 kµ t α µ /2 + Mt η . Proof. Introduce the solution operator u(t) = U (u0 , g,t). By linearity, u = u1 + u2 where u1 (t) = U (u0 , 0,t) and u2 (t) = U (0, g,t). In view of Corollary 10, it suffices to consider u1 . Let w(t) = u1 (t) − u0 so that w(0) = 0, and suppose to begin with that u0 ∈ Ḣ 2 (Ω). Using (5), we find that hw(t), vi + κ (I α ∇w)(t), ∇v − (~B1 w)(t), ∇v + h(B2 w)(t), vi = hρ (t), vi, 13 where ρ (t) = I α ∇ · (κ ∇u0 ) − ∇ · ~B1 u0 − ~B2 u0 . Since (I α u0 )′ (t) = u0 ωα (t), and recalling the definitions (6), we have
~ ρ ′ (t) = ∇ · (κ ∇u0 ) − ∇ · (~F(t)u0 ) − a(t)u0 ωα (t) − ∇ · (G(t)u 0 ) − b(t)u0, so kρ ( j+1)(t)k ≤ Cku0 k2t α − j . Therefore, by Corollary 10,
kw(m) (t)k + t α /2k∇w(m) (t)k ≤ Ct −m kw(0)k + ku0k2t α = Cku0 k2t α −m , which proves the result for integer-order time derivatives in the case µ = 2. Similarly, for the fractional-order time derivatives,
k∂tm−α w(t)k + t α /2k∂tm−α ∇w(t)k ≤ Ct α −m kw(0)k + ku0k2t α = Ct 2α −m ku0k2 , completing the proof for µ = 2. Since Corollary 10 also implies the case µ = 0, the result follows for 0 < µ < 2 by interpolation.  With the help of Lemma 16 from the Appendix, we can generalize Theorem 11 as follows. Theorem 13. Assume (28) and (29). If 0 ≤ µ ≤ 2 and u0 ∈ Ḣ µ (Ω), then
ku(m) (t)k2 ≤ Ct −m ku0 k µ t −(2−µ )α /2 + Mt η −α for 0 < t ≤ T . Proof. We know from Theorem11 that
ku(m) (t)k2 ≤ Ct −α −m ku0 k + Mt η , (30) so there is nothing to prove if u0 = 0. Thus, by linearity, we may assume that g(t) ≡ 0 and so M = 0. Integrating (1) in time, we see that u − ∇ · (κ ∇I α u) + ∇ · ~B1u + B2u = u0 and, after applying the operator ∂t I 1−α to both sides, −∇ · (κ ∇u) = ρ where ρ = ∂t I 1−α u0 − u − ∇ · ~B1u − B2u). Since −∇ · (κ ∇u(m) ) = ρ (m) in Ω, with u(m) (t) = 0 on ∂ Ω for 0 ≤ t ≤ T , it follows by H 2 -regularity for the Poisson problem that ku(m) (t)k2 ≤ Ckρ (m) (t)k. The identity (18) (with m = 1) implies that ρ = I 1−α ∂t u0 − u − ∇ · ~B1u − B2u)
~ + a∂t1−α u + bu , = −I 1−α u′ − I 1−α ∇ · (~F ∂t1−α u + Gu) and Lemma 16 and Theorem 12 (with µ = 2) imply that t m+1 k∂tm I 1−α u′ (t)k ≤ C max m ∑ s1−α +1+ j ku( j+1)(s)k 0≤s≤t j=0
≤ C max s1−α ku0 k2 sα = Ctku0 k2 . 0≤s≤t 14 (31) Since ∇ · (~ F ∂t1−α u) = (∇ · ~ F)∂t1−α u + ~F · ∂t1−α ∇u, we see from Lemma 16 and Theorem 12 that
F ∂t1−α u) (t) t m+1 ∂tm I 1−α ∇ · (~   m ≤ C max ∑ s1−α +1+ j k∂sj (∂s1−α u)k + k∂sj (∂s1−α ∇u)k 0≤s≤t j=0 m ∑ s1−α +1+ j 0≤s≤t = C max j=0   k∂sj+1−α uk + k∂sj+1−α ∇uk
≤ C max s1−α (sα + sα /2 ) ku0 k2 sα ≤ Ct 1+α /2 ku0 k2 . 0≤s≤t Similarly, t m+1 ∂tm I 1−α (a∂t1−α u) ≤ Ct 1+α /2ku0 k2 , whereas
~ (t) t m+1 ∂tm I 1−α ∇ · (Gu) 
≤ C max s2−α ku(s)k + k∇(u − u0)(s)k + k∇u0k 0≤s≤t  m
( j) ( j) 1−α +1+ j ku (s)k + k∇u (s)k +∑s j=1    ≤ C max s2−α ku0 k + s−α /2(ku0 k2 sα ) + k∇u0k 0≤s≤t
 + s2−α (1 + s−α /2) ku0 k2 sα ≤ Ct 2−α ku0 k2 and t m+1 ∂tm I 1−α (bu) ≤ Ct 2−α ku0 k2 . Thus,
t m+1 kρ (m) (t)k ≤ C t + t 1+α /2 + t 2−α ku0 k2 , showing that t m kρ (m) (t)k ≤ Cku0 k2 and therefore, by (30) and (31), ku(m) (t)k2 ≤ Ct −m−α ku0 k and ku(m) (t)k2 ≤ t −m ku0 k2 , which proves the result in the cases µ = 0 and µ = 2. The case 0 < µ < 2 then follows by interpolation.  Appendix A. Further technical lemmas The following identity uses the notation from Lemma 6. Lemma 14. Let µ > 0 and 1 ≤ q ≤ m. If, for m − q + 1 ≤ j ≤ m, j−(m−q) M j φ ∈ W1 (0, T ) with (∂tk M j φ )(0) = 0 for 0 ≤ k ≤ j − (m − q) − 1, then ∂tq M m I µ φ = m−q m, µ ∑ d˜j I µ +m−q− j M j φ + j=0 m ∑ j=m−q+1 15 j−(m−q) m, µ d˜j I µ ∂t M jφ . Proof. By (17), M mI µ = m−q ∑ ∂tq I µ +m− j If 0 ≤ j ≤ m − q, then m − q − j ≥ 0 so Therefore, m−q m, µ ∑ d˜j m, µ d˜j I µ +m− j M j . ∑ j=m−q+1 j=0 ∂tq m m, µ d˜j I µ +m− j M j + m−q I µ +m− j M j φ = j=0 m, µ ∑ d˜j = ∂tq I q I µ +m−q− j = I µ +m−q− j . I µ +m−q− j M j φ for φ ∈ L1 (0, T ). j=0 If m − q + 1 ≤ j ≤ m then j − (m − q) ≥ 1 so q−(m− j) m− j m− j µ I ∂t I ∂tq I µ +m− j = ∂t = ∂t j−(m−q) Iµ and thus m ∂tq m m, µ d˜j I µ +m− j M j φ = ∑ m, µ j−(m−q) µ d˜j ∂t I M jφ . ∑ j=m−q+1 j=m−q+1 By (18), ∂t j−(m−q) I µ M j φ = I µ ∂t j−(m−q) j−(m−q)−1 M jφ + ∑ (∂tk M j φ )(0) ωµ −k , k=0  and our hypotheses on φ ensure that all terms in the sum over k vanish. The next lemma was used in the proof of Lemma 7.
Lemma 15. Let ψ ∈ W∞2m−1 (0, T ); L∞ (Ω)d for some m ≥ 1 and let µ ≥ 0. Then, m µ µ, j Q 0,m (Bψ φ ,t) ≤ C ∑ Q2 (φ ,t) for 0 ≤ t ≤ T and φ ∈ Cαm . j=0 Proof. We integrate by parts m times to obtain µ 1− µ Bψ φ = I 1 (ψ∂t m−1 φ) = ∑ (−1)i ψ (i)I µ +i φ + (−1)mI 1 i=0 and so µ ,m
ψ (m) I µ +m−1 φ , m µ Bψ φ = (M m Bψ φ )(m) = ∑ (−1)i Bim φ , i=0 where Bim φ = ( If 0 ≤ i ≤ m − 1, then Bim φ = ∂tm
∂tm M m ψ (i) I µ +i φ for 0 ≤ i ≤ m − 1,
m m 1 (m) µ +m−1 ∂t M I ψ I φ for i = m. (i) m ψ (M I µ +i m
φ) =   m ∑ q ψ (i+m−q)∂tq M m I µ +i φ q=0 16 (A.1) so our assumption on ψ implies that m k(Bim φ )(t)k ≤ C ∑ ∂tq (M m I µ +i φ )(t) . (A.2) q=0 By Lemma 14, ∂tq M m I µ +i φ = m−q m, µ +i ∑ d˜j I µ +i+m−q− j M j φ j=0 m ∑ + m, µ +i µ +i j−(m−q) I d˜j ∂t M jφ , j=m−q+1 and by (16), I µM j = j j, µ M j−k I µ +k ∑ ck for µ > 0 and j ≥ 0, k=0 with ∂tq M j = q ∑ arj,qM j−r ∂tq−r = r=0 q ∑ ãrj,q M j−q+r ∂tr for 1 ≤ q ≤ j. r=0 Thus, for 0 ≤ j ≤ m − q, j I µ +i+m−q− j M j φ = j, µ +i+m−q− j ∑ ck M j−k I µ +i+m−q− j+k φ k=0 and for m − q + 1 ≤ j ≤ m, I µ +i ∂t j−(m−q) j−(m−q) M jφ = ∑ j, j−(m−q) ãr I µ +i M m−q M r ∂tr φ r=0 j−(m−q) ∑ = j, j−(m−q) ãr r=0 m−q m−q, µ +i ∑ ck M m−q−k I µ +i+k M r ∂tr φ , k=0 so q (∂t M m I µ +i φ )(t) 2 m−q ≤C ∑ (I µ +i+m−q− j M j φ )(t)k2 j=0 m ∑ +C (I µ +i ∂t j−(m−q) M j φ )(t) 2 j=m−q+1 m−q j ≤C ∑∑ (M j−k I µ +i+m−q− j+k φ )(t) 2 j=0 k=0 m +C ∑ j=m−q+1 j−(m−q) m−q ∑ ∑ r=0 k=0 17 2 (M m−q−k I µ +i+k M r ∂tr φ )(t) . µ Integrating in time, since Q 0 (M j I µ φ ,t) ≤ t 2 j Q2 (φ ,t), we see that m−q j
µ +i+m−q− j+k Q 0 ∂tq M m I µ +i φ ,t ≤ C ∑ ∑ t 2( j−k) Q2 (φ ,t) j=0 k=0 j−(m−q) m−q m +C ∑ j=m−q+1 µ +i+k,r ∑ ∑ t 2(m−q−k)Q2 r=0 (φ ,t) k=0 and therefore, by Lemma 4, q
µ ,r Q 0 ∂tq M m I µ +i φ ,t ≤ Ct 2(i+m−q) ∑ Q2 (φ ,t). (A.3) r=0 Hence, recalling (A.2), m m
µ ,r Q 0 (Bim φ ,t) ≤ C ∑ Q 0 ∂tq M m I µ +i φ ,t ≤ Ct 2i ∑ Q2 (φ ,t), for 0 ≤ i ≤ m − 1. q=0 r=0 (A.4) It remains to estimate Bmm φ = ∂tm M m I 1 (ψ (m) I µ +m−1 φ ). Taking q = m and µ = 1 in Lemma 14 gives m 1 j j ∂tm M m I 1 = d˜0m,1 I 1 + ∑ d˜m,1 j I ∂t M , j=1 and so m j−1 Bmm φ = d˜01,m I 1 (ψ (m) I µ +m−1 φ ) + ∑ d˜1,m (ψ (m) M j I µ +m−1 φ ). j ∂t j=1 Thus, Q 0 (Bmm φ ,t) ≤ CQ 0 I 1 (ψ (m) I µ +m−1 φ ),t
m j−1 +C ∑ ∑ Q0 j=1 q=0
∂tq M j I µ +m−1 φ ,t , and since 1 I (ψ (m) I µ +m−1
φ (t) 2 ≤ Z t kψ (m) Z t  µ +m−1 2 (s)k ds φ )(s)k ds kI 2 0 0 µ +m−1 ≤ CtQ2 (φ ,t) we have, by Lemma 4,
µ +m−1 µ Q 0 I 1 (ψ (m) I µ +m−1 φ ),t ≤ Ct 2 Q2 (φ ,t) ≤ Ct 2m Q2 (φ ,t). 18 Finally, using (A.3) with m replaced by j and with i replaced by m − 1, m j−1 q j=1 q=0 r=0 µ Q 0 (Bmm φ ,t) ≤ Ct 2m Q2 (φ ,t) + C ∑ ≤ Ct 2m m−1 µ ,r ∑ Q2 µ ,r ∑ t 2(m−1+ j−q) ∑ Q2 (φ ,t) (φ ,t). q=0 The result now follows from (A.1) and (A.4). We used the following result in the proof of Theorem 13.
Lemma 16. If m ≥ 0, ψ ∈ W∞m (0, T ); L∞ (Ω)d and M k φ ∈ W1k (0, T ) with  for 0 ≤ k ≤ m + 1, q (∂t M k φ )(0) = 0 for 1 ≤ q ≤ k − 1 and 1 ≤ k ≤ m + 1, then t m+1 ∂tm I µ (ψφ )(t) ≤ C max m ∑ ksµ +1+ j φ ( j) (s)k 0≤s≤t for 0 < t ≤ T . j=0 Proof. By (15), M m+1 ∂tm I µ (ψφ )k = m j ∂t M 1+ j I µ (ψφ ) ∑ b̃m+1,m j j=0 m (A.5) j ≤ C ∑ ∂t M j+1 µ I (ψφ )(t) , j=0 and in turn, ∂t j M j+1 I µ (ψφ ) = j+1 j+1, µ ∑ d˜k ∂t j I µ + j+1−k M k (ψφ ). k=0 j j Since ∂t I µ + j+1−k = ∂t (∂t I 1 )I µ + j+1−k = ∂tk (∂t for 0 ≤ k ≤ j + 1, j+1 j j+1−k I j+1−k )I µ +1 = ∂tk I µ +1 j k∂t M j+1 I µ (ψφ )(t)k ≤ C ∑ k∂t I µ + j+1−k M k (ψφ )k k=0 j+1 j+1 = C ∑ k∂tk I µ +1 M k (ψφ )(t)k = C ∑ kI µ +1 ∂tk M k (ψφ )(t)k (A.6) k=0 k=0 q where, in the last step, we used the fact that ∂t M k (ψφ )(0) = 0 for 0 ≤ q ≤ k − 1. We have k   k ∂tk M k (ψφ ) = ∂tk (ψ M k φ ) = ∑ ψ (k−q) ∂tq M k φ q q=0 q k   k k−(q−r) r =∑ ψ (k−q) ∑ ãk,q ∂t φ , r M q q=0 r=0 19 and hence k kI µ ∂tk M k (ψφ )(t)k ≤ C ∑ q ∑ k(I µ +1M k−q+r ∂tr φ )(t)k q=0 r=0 k =C ∑ q ∑ Z t q=0 r=0 0 k ≤C ∑  ωµ +1 (t − s)sk−q−µ −1 ksr+µ +1 φ (r) (s)k ds max ksr+µ +1 φ (r) (s)k r=0 0≤s≤t  k ∑ (ωµ +1 ∗ ωk−q−µ )(t). q=r Since (ωµ +1 ∗ ωk−q−µ )(t) = ωk−q+1 (t) ≤ Ct k−q . the result now follows from (A.5) and (A.6).  References [1] W. McLean, K. Mustapha, R. Ali, O. Knio, Well-posedness of time-fractional advection-diffusion-reaction equations(In review). [2] J. Klafter, I. M. Sokolov, First Steps in Random Walks, Oxford University Press, 2011. [3] C. N. Angstmann, B. I. Henry, B. A. Jacobs, A. V. McGann, A time-fractional generalised advection equation from a stochastic process, Chaos, Solitons and Fractals 102 (2017) 175–183, future Directions in Fractional Calculus Research and Applications. doi:10.1016/j.chaos.2017.04.040. [4] B. I. Henry, T. A. M. Langlands, P. Straka, Fractional Fokker–Planck equations for subdiffusion with space- and time-dependent forces, Phys. Rev. Lett. 105 (2010) 170602. doi:10.1103/PhysRevLett.105.170602. [5] B. I. Henry, S. L. Wearne, Fractional reaction-diffusion, Physica A: Statistical Mechanics and its Applications 276 (2000) 448–455. doi:10.1016/S0378-4371(99)00469-0. [6] T. A. M. Langlands, B. I. Henry, S. L. Wearne, Fractional cable equation models for anomalous electrodiffusion in nerve cells, SIAM J. Appl. Math. 71 (2011) 1168–1203. doi:10.1137/090775920. [7] F. Liu, V. V. Anh, I. Turner, P. Zhuang, Time fractional advectiondispersion equation, J. Appl. Math. Computing 13 (2003) 233–245. doi:10.1007/BF02936089. [8] R. Metzler, E. Barkai, J. Klafter, Deriving fractional Fokker–Planck equations from a generalised master equation, Europhys. Lett. 46 (1999) 431–436. doi:10.1209/epl/i1999-00279-7. [9] E. Cuesta, C. Lubich, C. Palencia, Convolution quadrature time discretization of fractional diffusive-wave equations, Math. Comp. 75 (2006) 673–696. doi:10.1090/S0025-5718-06-01788-1. 20 [10] B. Jin, B. Li, Z. Zhou, Discrete maximal regularity of time-stepping schemes for fractional evolution equations, Numer. Math. 138 (2018) 101–131. doi:10.1007/s00211-017-0904-8. [11] S. Karaa, A. K. Pani, Error analysis of a FVEM for fractional order evolution equations with nonsmooth initial data, ESAIM: M2AN 52 (2) (2018) 773–801. doi:10.1051/m2an/2018029. [12] K. Le, W. McLean, K. Mustapha, Numerical solution of the time-fractional Fokker–Planck equation with general forcing, SIAM J. Numer. Anal. 54 (2016) 1763–1784. doi:10.1137/15M1031734. [13] K. Le, W. McLean, K. Mustapha, A semidiscrete finite element approximation of a time-fractional Fokker–Planck equation with non-smooth initial data, SIAM J. Sci. Computing 40 (2018) A3831–A3852. doi:10.1137/17M1125261. [14] H.-L. Liao, D. Li, J. Zhang, Sharp error estimate of the nonuniform L1 formula for linear reaction-subdiffusion equations, SIAM J. Numer. Anal. 66 (2) (2018) 1112–1133. doi:10.1137/17M1131829. [15] K. Mustapha, Time-stepping discontinuous Galerkin methods for fractional diffusion problems, Numer. Math. 130 (3) (2015) 497–516. doi:10.1007/s00211-014-0669-2. [16] M. Stynes, E. O’Riordan, J. L. Gracia, Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation, SIAM J. Numer. Anal. 55 (2) (2017) 1057–1079. doi:10.1137/16M1082329. [17] W. McLean, Regularity of solutions to a time-fractional diffusion equation, ANZIAM J. 52 (2010) 123–138. doi:10.1017/S1446181111000617. [18] B. Li, X. Xie, Regularity of solutions to time fractional diffusion equations, Discrete Cont. Dyn. Syst. Ser. B 22 (1531-3492 2017 11 237) (2017) 1. doi:10.3934/dcdsb.2018340. [19] W. McLean, V. Thomée, Numerical solution via Laplace transforms of a fractional order evolution equation, J. Integral Equations Appl. 22 (1) (2010) 57–94. doi:10.1216/JIE-2010-22-1-57. [20] Z. Fan, Existence and regularity of solutions for evolution equations with Riemann–Liouville fractional derivatives, Indagationes Mathematicae 25 (3) (2014) 516–524. doi:10.1016/j.indag.2014.01.002. [21] V. Keyantuo, C. Lizama, M. Warma, Existence, regularity and representation of solutions of time fractional diffusion equations, Adv. Differential Equations 21 (9/10) (2016) 837–886. URL https://projecteuclid.org:443/euclid.ade/1465912585 21 [22] K. Sakamoto, M. Yamamoto, Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems, Journal of Mathematical Analysis and Applications 382 (1) (2011) 426–447. doi:10.1016/j.jmaa.2011.04.058. [23] J. Mu, B. Ahmad, S. Huang, Existence and regularity of solutions to timefractional diffusion equations, Comput. Math. Appl. 73 (6) (2017) 985–996. doi:10.1016/j.camwa.2016.04.039. [24] K. Le, W. McLean, M. Stynes, Existence, uniqueness and regularity of the solution of the time-fractional Fokker–Planck equation with general forcing, Comm. Pure Appl. Math.(To appear). [25] J. A. Nohel, D. F. Shea, Frequency domain methods for Volterra equations, Advances in Mathematics 22 (1976) 278–304. doi:10.1016/0001-8708(76)90096-7. [26] J. Dixon, S. McKee, Weakly singular Gronwall inequalities, ZAMM Z. Angew. Math. Mech. 66 (1986) 535–544. [27] V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, 2nd Edition, Springer, 2006. [28] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Classics in Applied Mathematics, SIAM, 2011. doi:10.1137/1.9781611972030. 22