{'fraction_non_alphanumeric': 0.07515384879843229, 'fraction_numerical': 0.026884862653419054, 'mean_word_length': 3.5300576234231427, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 40, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Time fractional gradient flows: Theory and numerics

We develop the theory of fractional gradient flows: an evolution aimed at the minimization of a convex, l.s.c.~energy, with memory effects. This memory is characterized by the fact that the negative of the (sub)gradient of the energy equals the so-called Caputo derivative of the state. We introduce the notion of energy solutions, for which we provide existence, uniqueness and certain regularizing effects. We also consider Lipschitz perturbations of this energy. For these problems we provide an a posteriori error estimate and show its reliability. This estimate depends only on the problem data, and imposes no constraints between consecutive time-steps. On the basis of this estimate we provide an a priori error analysis that makes no assumptions on the smoothness of the solution.


Introduction
In recent times problems involving fractional derivatives have garnered considerable attention, as it is claimed that they better describe certain fundamental relations between the processes of interest; see, for instance [29,15,46]. In this, and many other references the models considered are linear. However, it is well known that real world phenomena are not linear, not even smooth. It is only natural then to consider nonlinear/nonsmooth models with fractional derivatives.
The purpose of this work is to develop the theory and numerical analysis of so-called timefractional gradient flows: an evolution equation aimed at the minimization of a convex and lower semicontinuous (l.s.c.) energy, but where the evolution has memory effects. This memory is characterized by the fact that the negative of the (sub)gradient of the energy equals the so-called Caputo derivative of the state.
The Caputo derivative, introduced in [11], is one of the existing models of fractional derivatives. It is defined, for α ∈ (0, 1), by (1.1) D α c w(t) = 1 Γ(1 − α)ˆt 0ẇ (r) (t − r) α dr, where Γ denotes the Gamma function. This definition, from the onset, seems unnatural. To define a derivative of a fractional order, it seems necessary for the function to be at least differentiable. Below we briefly describe several attempts at circumventing this issue. We focus, in particular, on the results developed in a series of papers by Li and Liu, see [25,28,26,27], where they developed a distributional theory for this derivative; see also [16]. The authors of these works also constructed, in [26], so-called deconvolution schemes that aim at discretizing this derivative. With the help of this definition and the schemes that they develop the authors were able to study several classes of equations, in particular time fractional gradient flows.
Let us be precise in what we mean by this term. Let T > 0 be a final time, H be a separable Hilbert space, Φ : H → R ∪ {+∞} be a convex and l.s.c. functional, which we will call energy. Given u 0 ∈ H, and f : (0, T ] → H we seek for a function u : [0, T ] → H that satisfies (1.2) D α c u(t) + ∂Φ(u(t)) ∋ f (t), t ∈ (0, T ], u(0) = u 0 , where by ∂Φ we denote the subdifferential of Φ. Our objectives in this work can be stated as follows: We will introduce the notion of "energy solutions" of (1.2), and we will refine the results regarding existence, uniqueness, and regularizing effects provided in [28]. This will be done by generalizing, to non-uniform time steps the "deconvolution" schemes of [26,28], and developing a sort of "fractional minimizing movements" scheme. We will also provide an a priori error estimate that seems optimal in light of the regularizing effects proved above. We also develop an a posteriori error estimate, in the spirit of [30] and show its reliability.
We comment, in passing, that nonlinear evolution problems with fractional time derivative have been considered in other works. From a modeling point of view, their advantages have been observed in [15,12]. Some other types of nonlinear problems have been studied in [8,40,2,24,23,39,45] and [31,38] where, for a particular type of nonlinear problem other "energy dissipation inequalities" than those we obtain are derived. Regularity properties for nonlinear problems with fractional time derivatives have been obtained in [22,14,21,1,44,43,42,41]. Of particular interest to us are [28] which we described above and [3] which also considers time fractional gradient flows. The assumptions on the data, however, are slightly different than ours. As such, some of the results in [3] are stronger, and some weaker than ours; in particular, we conduct a numerical analysis of this problem. Nevertheless, we refer to this reference for a nice historical account and particular applications to PDEs.
Our presentation will be organized as follows. We will establish notation and the framework we will adopt in Section 2. Here, in particular, we will study several properties of a particular space, which we denote by L p α (0, T ; H), and that will be used to characterize the requirements on the right hand side f of (1.2). In addition, we also review the various proposed generalizations of the classical definition of Caputo derivatives, with particular attention to that of [25,28,27]; since this is the one we shall adopt. In Section 3 we generalize the deconvolution schemes of [26,28] and their properties, to the case of nonuniform time stepping. Many of the simple properties of these schemes are lost in this case, but we retain enough of them for our purposes. Section 4 introduces the notion of energy solutions for (1.2) and shows existence and uniqueness of these. This is accomplished by introducing, on the basis of our generalized deconvolution formulas, a fractional minimizing movements scheme; and showing that the discrete solutions have enough compactness to pass to the limit in the size of the partition. In Section 5 we provide an error analysis of the fractional minimizing movements scheme. First, we show how an error estimate follows as a side result from the existence proof. Then, in the spirit of [30], we provide an a posteriori error estimator for our scheme and show its reliability. This estimator is then used to independently show rates of convergence. This section is concluded with some particular instances in which the rate of convergence can be improved. Section 6 is dedicated to the case in which we allow a Lipschitz perturbation of the subdifferential. We extend the existence, uniqueness, a priori, and a posteriori approximation results of the fractional gradient flow. Finally, Section 7 presents some simple numerical experiments that illustrate, explore, and expand our theory.

Notation and preliminaries
Let us begin by presenting the main notation and assumptions we shall operate under. We will denote by T ∈ (0, ∞) our final (positive) time. By H we will always denote a separable Hilbert space with scalar product ·, · and norm · . As it is by now customary, by C we will denote a nonessential constant whose value may change at each occurrence.
We will always assume that our energy is bounded from below, that is As we are not assuming smoothness in our energy beyond convexity, a useful substitute for its derivative is the subdifferential, that is, The effective domain of the subdifferential is D(∂Φ) = {w ∈ H : ∂Φ(w) = ∅} . Recall that, in our setting, we always have that D(∂Φ) = D(Φ). We refer the reader to [13,33] for basic facts on convex analysis.
In applications, it is sometimes useful to obtain error estimates on (semi)norms stronger than those of the ambient space, and that are dictated by the structure of the energy. For this reason, we introduce the following coercivity modulus of Φ, see [30,Definition 2.3].

2.2.
Vector valued time dependent functions. We will follow standard notation regarding Bochner spaces of vector valued functions, see [32,Section 1.5]. For any w ∈ L 1 (0, T ; H) and E ⊂ [0, T ] that is measurable, we define the average by where |E| denotes the Lebesgue measure of E.
Since eventually we will have to deal with time discretization, we also introduce notation for time-discrete vector valued functions. Let P be a partition of the time interval [0, T ] (2.2) P = {0 = t 0 < t 1 < . . . < t N −1 < t N = T }, with variable steps τ n = t n − t n−1 and τ = max{τ n : n ∈ {1, . . . , N }}. We will always denote by N the size of a partition. For t ∈ [0, T ] we define ⌊t⌋ P = max{r ∈ P : r < t}, ⌈t⌉ P = min{r ∈ P : t ≤ r}, and n(t) to be the index of ⌈t⌉ P , so that t ∈ (⌊t⌋ P , ⌈t⌉ P ] = (t n(t)−1 , t n(t) ]. Given a partition P, To quantify the assumptions we need on the right hand side f of (1.2) we introduce the following space.
The proof of the first embedding is a simple application of Hölder inequality. Indeed, we have and hence as we intended to show.
When dealing with discretization we will approximate the right hand side f of (1.2) by its local averages over a partition P. Thus, we must provide a bound on this operation that is independent of the partition.
Then, there exists a constant C > 0 only depending on p and α such that . We first, for n ∈ {1, . . . , N }, bound the integral tn 0 (t n − s) α−1 F P (s) p ds.
To achieve this, we decompose this integral as We use Hölder inequality in the definition of F k to obtain that Since, for every p ∈ (1, ∞) the function s → s α−1 belongs to the Muckenhoupt class A p (R + ), see [20, Example 7.1.7], there exists a constant C p,α that only depends on p and α such that Therefore, for any k, we have (2.8) Substituting (2.7) and (2.8) into (2.6) we get Now consider t ∈ [0, T ]. Taking advantage of the estimate we obtained above we write Therefore by taking supremum over t ∈ [0, T ] and C = (2C p,α ) 1/p , we finish the proof of this lemma.
For p = 1, the proof proceeds almost the same way as before. The only difference worth noting is that, instead of (2.7), we have Next, we observe that, since α−1 ∈ (−1, 0), then the function s → s α−1 belongs to the Muckenhoupt class A 1 (R + ). Thus, With this information, the proof proceeds without change.
It turns out that averaging is not only continuous, but possesses suitable approximation properties in this space. Namely, we have a control on the difference between fractional integrals of f ∈ L p α (0, T ; H) and its averages.
where the constants C, C ′ depend only on p and α. In addition, for any β ∈ (0, 1) we also have where the constants C 1 , C ′ 1 depend on p, α, and β. As usual, when p = 1, we have p ′ = ∞ and 1/p ′ is treated as 0.
Proof. We first notice that the second inequalities in both (2.10) and (2.11) follow directly from Lemma 2.4 and the triangle inequality.
To show the first inequality in (2.10), given P we consider t ∈ [0, T ]. Using that f − F P has zero mean on each subinterval of the partition, we can write For the first term, denoted I 1 (t), we have where C 1 only depends on p and α. For the second term, noticing that t − t k−1 + τ > t − s for s ∈ (t k−1 , t k ) we have we obtain , and (2.10) follows after combining the bounds for I 1 (t) and I 2 (t) that we have obtained.
To prove (2.11) we apply the Hölder inequality to (2.12) with α replaced by β to get Arguing as in the bound for I 2 (t) Thus, to obtain (2.11) it suffices to show that, for every r ∈ [0, T ], with some constant C 2 only depending on p, α, and β. To estimate the fractional integral of II 2 by Fubini's theorem we have where we set a ∧ b = min{a, b}. We claim that there exists a constant C 3 depending on α and β such that On the one hand, for r − s ≤ 2τ , we simply havê ⌈s⌉P ∧r On the other hand, if r − s > 2τ , then Therefore (2.14) is proved, and thus (2.14) implies that To conclude, we claim that for a constant C 4 depending on α and β. Indeed, if this is the case, we havê α (0,T ;H) , and we combine the estimates for II 2 (t) and II 3 (t) together and conclude the proof of (2.11).
Let us now turn to the proof of (2.15). First, if r − s ≤ τ then it suffices to observe that The first term can be bounded using that r − s > τ as follows On the other hand, since for t ∈ (0, τ ) we have that r − s + τ − t ≥ r − s, the second term can be estimated asˆτ This concludes the proof.
We refer the reader to [27, section 4] for further results concerning the space L p α (0, T ; H).
2.3. The Caputo derivative. As we mentioned in the Introduction, the definition of the Caputo derivative, given in (1.1) seems unnatural. Smoothness of higher order is needed to define a fractional derivative. Several attempts at resolving this discrepancy have been proposed in the literature and we here quickly describe a few of them. First, one of the main reasons that motivate practitioners to use, among the many possible definitions, the Caputo derivative (1.1) is, first, that D α c 1 = 0 and second that this derivative allows one to pose initial value problems like (1.2). However, it is by now known that even in the linear case solutions of problems involving the Caputo derivative possess a weak singularity in time [37,36,35]. This singular behavior of the solution forces one to wonder: If fractional derivatives describe processes with memory, why is it sufficient to know the state at one particular point (initial condition) to uniquely describe the state at all future times? Is it possible that the singularity is precisely caused by the fact that we are ignoring the past states of the system? This motivates the following: Set w(t) = w 0 for t ≤ 0. Therefore, where, in the last step, we integrated by parts. The expression D α m w(t) is known as the Marchaud derivative of order α of the function w. This is the way that the Caputo derivative has been understood, for instance, in [6,5,7,4]. We comment, in passing, that owing to [9] this fractional derivative satisfies an extension problem similar to the (by now) classical Caffarelli Silvestre extension [10,34] for the fractional Laplacian.
Another approach, and the one we shall adopt here, is to notice that (1.1) can be converted, for sufficiently smooth functions, into a Volterra type equation This identity is the beginning of the theory developed in [25] to extend the notion of Caputo derivative. To be more specific, [25] considers the set of distributions for a fixed time T > 0. Then the modified Riemann Liouville derivative for any distribution w ∈ E T is defined, following classical references like [18, Section 1.5.5], as , with θ being the Heaviside function, is a distribution supported in [0, ∞) and the convolution is understood as the generalized definition between distributions. Here D denotes the distributional derivative. Reference [25] then uses this to define the generalized Caputo derivative of w ∈ L 1 loc ([0, T ); H) associated with w 0 by D α c w = D α rl (w − w 0 ). If there exists w(0) ∈ H such that lim t↓0 ffl t 0 w(s) − w(0) ds = 0, then we always impose w 0 = w(0) in this definition. It is shown in [25,Theorem 3.7] that for such a function w, (2.17) holds for Lebesgue a.e. t ∈ (0, T ) provided that the generalized Caputo derivative D α c w ∈ L 1 loc ([0, T ); H). We also comment that [25, Proposition 3.11(ii)] implies that for every function w ∈ L 2 (0, Finally, we recall that the Mittag-Leffler function of order α ∈ (0, 1) is defined via .
We refer the reader to [19] for an extensive treatise on this function. Here we just mention that this function satisfies, for any λ ∈ R, the identity 2.3.1. An auxiliary estimate. Having defined the Caputo derivative of a function, we present an auxiliary result. Namely, an estimate on functions that have piecewise constant, over some partition P, Caputo derivative.
Lemma 2.6 (continuity). Let p ∈ [1, ∞); P be a partition, as in (2.2), of [0, T ]; and w ∈ L 1 (0, T ; H) be such that its generalized Caputo derivative D α c w ∈ L p α (0, T ; H), and it is piecewise constant over P. Then we have (2.20) sup where the constant C depends only on α.
Proof. The representation (2.17) allows us to write Therefore by Hölder inequality, we have where the constants C, C 1 , and C 2 depend only on p and α. For I 2 (t), we simply havê . Now to bound the integral for I 1 (t), we use Fubini's theorem to get We claim that where C 3 only depends α. If this is true, then we havê . The proof of (2.21) proceeds as the one for (2.15). For brevity we skip the details.
2.4. Some comparison estimates. As a final preparatory step we present some auxiliary results that shall be repeatedly used and are related to differential inequalities involving the Caputo derivative, and a Grönwall-like lemma.
First, we present a comparison principle which is similar to [17,Proposition 4.2]. The proof can be done easily by contradiction, and therefore it is omitted here. Lemma 2.7 (comparison). Let g 1 , g 2 : [0, T ]×R → R be both nondecreasing in their second argument and g 2 be measurable. Assume that v, w ∈ C([0, T ]; R) satisfy v(0) < w(0), and there is some We now present a result that can be interpreted as an extension of [30,Lemma 3.7] to the fractional case. However, unlike the classical case, here we have the restriction that λ ≥ 0 because we have to argue from a fractional integral inequality. Nevertheless, this is sufficient for our purposes.
+∞] be measurable functions, and λ ≥ 0. If the following differential inequality is satisfied Proof. From (2.22) we obtain that (2.24) where a(t) = max 0≤s≤t a(s) and the functions C, D are defined in (2.23). This immediately implies that In order to bound a, we construct a barrier function and hence for every t ∈ (0, T ) provided that Applying Lemma 2.7 we obtain that Plugging this back into (2.24) and noticing that this holds for any K satisfying (2.25) we obtain that which is the desired result.

Deconvolutional discretization of the Caputo derivative
To discretize the Caputo fractional derivative, references [26, 28] consider a so-called deconvolutional scheme on uniform time grids and prove some properties of this discretization. In this section, we generalize this deconvolutional scheme to the variable time step setting, and prove properties that will be useful in deriving a posteriori error estimates later, in Section 5.2.
3.1. The discrete Caputo derivative. Let P be a partition as in (2.2). To motivate this discretization, let us assume that w : . Then formally by (2.17), we have Let K P ∈ R N ×N be the matrix induced by the partition P, which is defined as Then we can rewrite (3.1) in matrix form as where V, W, W 0 ∈ H N with V n = V n , W n = w(t n ), and (W 0 ) n = w(0). Notice that K P is lower triangular and all the elements on and below the main diagonal are positive. Therefore K P is invertible and its inverse is also lower triangular. Thus, the previous identity is equivalent to where we set K −1 P,n0 = − n j=1 K −1 P,nj . This motivates the following approximation of the Caputo derivative provided W ∈ H N and W 0 ∈ H are given. For n ∈ {1, . . . , N } we set 3.2. Properties of K −1 P . We note that, when the partition is uniform, both K P and its inverse will be Toeplitz matrices, and hence the product K P V can be interpreted as the convolution of sequences. Consequently, multiplication by K −1 P is equivalent to taking a sequence deconvolution. This motivates the name of this scheme and enables [28] to apply techniques for the deconvolution of a completely monotone sequence and prove properties of K −1 P . We were not successful in extending, to a general partition P, all the properties of K −1 P presented in [28] for the case when the partition is uniform. This is mainly because their techniques are based on ideas that rely on completely monotone sequences, which do not easily extend to a general P. Nevertheless we have obtained sufficient, for our purposes, properties. The following result is the counterpart to [28, Proposition 3.2(1)]. Proposition 3.1 (properties of K −1 P ). Let P be a partition as in (2.2), and K P be defined in (3.2). The matrix K P is invertible, and its inverse satisfies: Proof. We already showed that K P is nonsingular. We prove (3.4) and (3.5) separately.
First, to prove that K −1 P,n0 < 0. For this, it suffices to show that for a vector W ∈ R N such that W i = 1 for any i ≥ 1, then the vector We prove this by induction on n. For n = 1, clearly Suppose that F j > 0 for all 1 ≤ j ≤ k, now we want to show that F k+1 > 0 as well. Notice that then taking the difference we have We claim that K P,k+1,j − K P,k,j < 0 for any j. In fact, this can be seen through the definition of the entries of K P Using K P,k+1,j −K P,k,j < 0 and F j > 0 for all j ∈ {1, . . . , k} in (3.6), we see that K P,k+1,k+1 F k+1 > 0 and thus F k+1 > 0. Therefore by induction we proved that K −1 P,n0 < 0 for n ≥ 1. Next, we prove that K −1 P,ii > 0 and K −1 P,ni < 0. Consider a vector W ∈ R N that is such that W i = 1 and W j = 0 for j = i. It suffices to prove that for, and thus F i = 1/K P,ii > 0. Now we prove by induction that (3.7) holds. First, when n = i + 1, we have and hence This shows that (3.7) is true for n = i + 1. Now suppose that we have already shown that F n < 0 for n satisfying n ∈ {i + 1, . . . , k}, we want to prove F k+1 < 0. To this aim, notice that and thus, since K P,k,i > 0, we can get Since by the induction hypothesis F j < 0 for j ∈ {i + 1, . . . , k}, it only remains to show that Due to j > i, we have ξ < η and hence Therefore from the arguments above we see that F k+1 < 0, and by induction K −1 P,ni < 0 for n > i.
Remark 3.2 (generalization). The discretization of the Caputo derivative, described in (3.3), and its properties presented in Proposition 3.1 can be extended to more general kernels. Indeed, for a general convolutional kernel g ∈ L 1 (0, T ; R) the entries of the matrix K P will be The proof of (3.4) follows verbatim provided g ′ (t) < 0, as the reader can readily verify. The proof of (3.5) only requires that the function   .2), and K P be defined as in (3.2). Then, its inverse satisfies: 1. For n ∈ {1, . . . , N − 1}, Proof. To prove (3.8) it suffices to show that for a vector W ∈ R N such that W i = 1 for any i ≥ 1, then the vector We prove this by induction on n. For n = 1, Clearly, Hence we have where we set F 0 = 0 in the equations above. Therefore to show F k+1 < F k , we only need to prove that Since we also have Taking the difference between the equation above and the one for W k , we obtain that In light of this identity, we claim that to obtain (3.10) it suffices to show that To finish the proof, we focus on (3.11), fix i and define Applying Cauchy's mean-value theorem to the two fractions above, we know there exists η ∈ (x + c 2 , x + c 3 ) and ξ ∈ (x + c 1 , x + c 2 ) such that where the last inequality holds because α < 1 and ξ < x + c 2 < η. This shows the monotonicity of function h and confirms (3.11). This concludes the inductive step and proves (3.8).
The proof of (3.9) is obtained similarly. For convenience we only write the proof for i = 1, but the extension to general i is straightforward. Consider a vector W ∈ R N such that W j = 1 if j = 1 and W j = 0 if j = 1, then it suffices to prove that vector F = K −1 P W satisfies (3.12) F n < F n+1 for n ∈ {2, . . . , N − 1}. We prove (3.12) by induction on n. For n = 2, observe that from the proof of (3.8) with F 0 = 0, we have From the first and second equation above, we see that F 1 > 0 and F 2 − F 1 < 0. Combining the second and the third equation we deduce that It also remains to prove that when (3.12) holds for n ∈ {2, . . . , k − 1}, then it also holds for n = k, i.e. F k < F k+1 , provided that k < N . To this aim, we first see that Therefore in order to prove F k < F k+1 , we only need to show that (3.14) Similar to (3.13) we also have Thanks to the inductive hypothesis, we know that F j+1 − F j < 0 for j = 2 and F j+1 − F j > 0 for j ∈ {3, . . . , k − 1}, Therefore using a similar argument used in the proof for (3.8), to prove (3.14) we only need to show which is similar to (3.11). We rewrite the inequality above as and define the function This implies that h ′ 1 (x) < 0 for 0 < x < t k−1 and finishes inductive step of the induction. Hence (3.9) is proved.
Remark 3.4 (generalization). Notice that, for a general kernel g, property (3.8) remains valid provided G(t) = ln(g(t)) satisfies G ′′ (t) > 0. Figure 1. Given a partition P, the figure shows the nonlocal basis functions {ϕ P,i } N i=0 for different values of α. Every function whose Caputo derivative is piecewise constant can be written as a linear combination of these functions. Notice that, for any partition point ϕ P,i (t j ) = δ ij . In addition, Proposition 3.5 shows that these functions form a partition of unity.

A continuous interpolant. Given a partition
By definition, we have that W P (t n ) = W n . Moreover, where we defined (3.19) ϕ P,0 (t) = 1 + 1 Γ(α + 1) The functions {ϕ P,i } N i=0 play the role, in this context, of the standard "hat" basis functions used for piecewise linear interpolation over a partition P. Indeed, they are such that any function with piecewise constant (Caputo) derivative can be written as a linear combination of them. Figure 1 illustrates the behavior of these functions. As expected, and in contrast to the hat basis functions, these functions are nonlocal, in the sense that they have global support. Something worth noticing is also that the figure seems to indicate that, as α ↓ 0, the functions resemble piecewise constants and, in contrast, when α ↑ 1 they tend to the classical hat basis functions.
An important feature of the hat basis functions is that they form a partition of unity. It is easy to check that, for any t ∈ [0, T ] we have n(t) i=0 ϕ P,i (t) = 1. The following result shows that ϕ P,i (t) ≥ 0. Thus, for any t ∈ [0, T ], W P (t) is a convex combination of its nodal values {W j } N j=0 . This observation will be crucial to derive an a posteriori error estimate in Section 5.2.
Proof. By definition, for t = t n , we have ϕ P,n (t n ) = 1 and ϕ P,i (t n ) = 0 for any i = n. Also, for i > n(t), we see that ϕ P,i (t) = 0, and hence it only remains to show that ϕ P,i (t) > 0 for i ≤ n(t).
To show this, consider W i = 1 and W j = 0 for j = i, a piecewise constant V P and its interpolation W P defined in (3.16) and (3.17). Then our goal is to show that W P (t) > 0. If i = n(t) > 0, then it is easy to check by definition that (D α P W) n > 0 and (D α P W) j = 0 for j ∈ {1, . . . , i − 1}. Therefore we obtain If i < n(t), the proof is not that straightforward. The trick is to insert the time t, which is not on the partition P, to get a new partition P ′ = P ∪ {t} and then apply Propositions 3.1 and 3.3 in an appropriate way. Let us now work out the details. Let P ′ = {t ′ k } N +1 k=0 and notice that t ′ n(t) = t, t ′ n(t)+1 = t n(t) . On the basis of this partition we define the vector Since the only possible nonzero components of W ′ are W ′ i = W i = 1 and W ′ n(t) = W P (t), therefore we deduce from the equality above that , which can be rearranged as From Proposition 3.3 we see that K −1 P ′ ,n(t)+1,i − K −1 P ′ ,n(t)i > 0 and from Proposition 3.1 we see that K −1 P ′ ,n(t)n(t) − K −1 P ′ ,n(t)+1,n(t) > 0 as a consequence of K −1 P ′ ,n(t)n(t) > 0 and K −1 P ′ ,n(t)+1,n(t) < 0. This leads to the fact that W P (t) > 0 and finishes our proof.

Time fractional gradient flow: Theory
We have now set the stage for the study of time fractional gradient flows, which were formally described in (1.2). Throughout the remaining of our discussion we shall assume that the initial condition satisfies u 0 ∈ D(Φ) and that f ∈ L 2 α (0, T ; H). We begin by commenting that the case f = 0 was already studied in [28, Section 5] where they studied so-called strong solutions, see [28, Definition 5.4]. Here we trivially extend their definition to the case f = 0.

Energy solutions.
Since H is a Hilbert space, we will mimic the theory for classical gradient flows and introduce the notion of energy solutions for (1.2). To motivate it, suppose that at some t ∈ (0, T ) , then, by definition of the subdifferential, this is equivalent to the evolution variational inequality (EVI) Notice that, provided u 0 ∈ D(Φ) we can set w(t) = u 0 in (4.2) and obtain that´T 0 Φ(u(t))dt < ∞, which motivates the name for this notion of solution. In addition, as the following result shows, any energy solution is a strong solution. Proof. Evidently, it suffices to prove that that f (t) − D α c u(t) ∈ ∂Φ(u(t)) for almost every t ∈ (0, T ). Let w 0 ∈ H, t 0 ∈ (0, T ), and choose h > 0 sufficiently small so that (t 0 − h, t 0 + h) ⊂ (0, T ). Define where by χ S we denote the characteristic function of the set S. This choice of test function on (4.2) gives The assumptions of an energy solution guarantee that all terms inside the integrals belong to L 1 (0, T ; R) so that for almost every t 0 we have, as h ↓ 0, that which is (4.1) and, as we intended to show, is equivalent to the claim.
Remark 4.4 (coercivity). By introducing the coercivity modulus of Definition 2.1 one realizes that an energy solution u satisfies, instead of (4.1) and (4.2), the stronger inequalities and, for any w ∈ L 2 (0, T ; H),

4.2.
Existence and uniqueness. In this section, we will prove the following theorem on the existence and uniqueness of energy solutions to (1.2) in the sense of Definition 4.2. The main result that we will prove reads as follows.
Theorem 4.5 (well posedness). Assume that the energy Φ is convex, l.s.c., and with nonempty effective domain. Let u 0 ∈ D(Φ) and f ∈ L 2 α (0, T ; H). In this setting, the fractional gradient flow problem (1.2) has a unique energy solution u, in the sense of Definition 4.2. For almost every t ∈ (0, T ), the solution u satisfies that f (t) − D α c u(t) ∈ ∂Φ(u(t)) and for any t ∈ [0, T ] we have In addition, u ∈ C 0,α/2 ([0, T ]; H) with modulus of continuity where the constant C depends only on α.
We point out that our assumptions are weaker than those in [28, Theorem 5.10]. First, we allow for a nonzero right hand side. In addition, we do not require [28, Assumption 5.9], which is a sort of weak-strong continuity of subdifferentials.
The remainder of this section will be dedicated to the proof of Theorem 4.5. To accomplish this, we follow a similar approach to [28, Section 5]. To show existence of solutions, we consider a sort of fractional minimizing movements scheme. We introduce a partition P with maximal time step τ and compute the sequence U = {U n } N n=0 ⊂ H as follows. Assume U 0 ∈ D(Φ) is given, the n-th iterate, for n ∈ {1, . . . , N }, is defined recursively via We will usually choose U 0 = u 0 , but other choices of U 0 ∈ D(Φ) are also allowed. From the approximation scheme (4.7) and the expression of the discrete Caputo derivative (D α P U) n given in (3.3), it is clear that Thanks to Proposition 3.1, for i = 0, . . . , n − 1, we have that K −1 P,ni < 0 and as a consequence the functional on the right hand side of (4.9) is uniformly convex. Combining with the fact that Φ is lower semicontinuous, the functional on the right hand side has a unique minimizer, and hence U n is well-defined. Now, in order to define a continuous in time function from U, we use the interpolation introduced in (3.16). Let V P (t) = (D α P U) n(t) . Then we have Recall that F P can be defined from {F n } N n=1 using (2.3) and that Lemma 2.4 showed that F P ∈ L 2 α (0, T ; H) with a norm bounded independently of P. We now obtain some suitable bounds for U P and V P . Lemma 4.6 (a priori bounds). Let P be any partition. The functions U P and V P satisfy where the constant C only depends on α.
Therefore noticing that K −1 P,ni < 0 for i ∈ {0, . . . , n − 1}, we get where we denoted Φ(U) = {Φ(U n )} N n=0 . We can now proceed to obtain the claimed estimates. To prove the first one, we use that to obtain that for any n, where the constant C depends only on α. Now, since Proposition 3.5 has shown that U P is a convex combination of the values U n , we have which finishes the proof of the first claim.
We now proceed to prove the second claim. Using (4.12) we get for any t ∈ [0, T ]. This implies that which, using Lemma 2.4, implies the result.
Remark 4.7 (the function Φ). Notice that, during the course of the proof of the first estimate in (4.11) we also showed that, if we define Φ P (t) = N i=0 ϕ P,i (t)Φ(U i ), then Φ(t) is the interpolation of Φ P (U) with piecewise constant Caputo derivative. Moreover, These estimates immediately yield a modulus of continuity estimate on the interpolant U P which is independent of the partition P.
Lemma 4.8 (Hölder continuity). Let P be any partition and U ∈ H N be the solution to (4.7) associated to this partition. For t 1 , t 2 ∈ [0, T ] the interpolant U P , defined in (3.16), satisfies where the constant C depends only on α.
Proof. As proved in [28, Lemma 5.8], D α c w ∈ L 2 α (0, T ; H) guarantees w ∈ C 0,α/2 ([0, T ]; H). Therefore using D α c U = V α ∈ L 2 α (0, T ; H) and the estimate from Lemma 4.6, we obtain the result. Next we control the difference between discrete solutions corresponding to different partitions. Lemma 4.9 (equicontinuity). Let, for i = 1, 2, P i be partitions of [0, T ] with maximal step size τ i , respectively, and denote by U (i) the associated solutions to (4.7). Let U i be their interpolations, defined by (4.10), and U i be their piecewise constant interpolations as in (2.3). Assuming that U where the constant C only depends on α.
Proof. For almost every t ∈ [0, T ], we have that where to bound I we used that F i (t) − D α c U i (t) ∈ ∂Φ(U i (t)) and Definition 2.1. Define now so that D α c G(t) = F 1 (t) − F 2 (t) and by (2.10) of Lemma 2.5 one further has (4.16) where C is a constant that depends only on α. Using these estimates, from (4.15) we deduce that and, using (2.17) and (4.16), we then conclude . It remains then to estimate the fractional integral on the right hand side. We estimate each term separately.
Recalling that U i (t) = U i (⌈t⌉ i ) we can invoke Lemma 2.6 and, again, Lemma 4.6 to arrive at Combining all estimates proves the desired result.
We are finally able to prove Theorem 4.5. We will follow the same approach as in [28, Theorem 5.10]; we will pass to the limit τ i ↓ 0 and study the limit of discrete solutions U i .
Proof of Theorem 4.5. Let us first prove uniqueness of energy solutions. Suppose that we have two energy solutions u 1 , u 2 to (1.2). Let t ∈ (0, T ) be arbitrary and h > 0 be sufficiently small so Setting as test function, in the EVI that characterizes u 1 , the function w = u 1 − χ (t−h,t+h) (u 1 − u 2 ) and vice versa, and adding the ensuing inequalities we obtain Combining with the facts that d ≥ 0 and ffl t 0 |d(s)|ds → 0 we obtain, by [25,Corollary 3.8], d(t) = 0. This proves the uniqueness.
We now turn our attention to existence. Let {P k } ∞ k=1 be a sequence of partitions such that τ k ↓ 0 as k → ∞. We denote by U (k) the discrete solution, on partition P k , given by (4.7) with U (k) 0 = u 0 . The symbols U k , V k and F k carry analogous meaning. Owing to Lemma 4.9 there exists u ∈ C([0, T ]; H) such that U k converges to u in C([0, T ]; H).
The embedding of Proposition 2.3 and an application of Lemma 4.6 shows that there is a subsequence for which V kj ⇀ v in L 2 (0, T ; H) as j → ∞. Moreover, we can again appeal to Lemma 4.6 to see that, for every t ∈ [0, T ], the sequence is uniformly bounded in L 2 (0, t; H) so that by passing to a further, not retagged, subsequence This, in addition, shows that v ∈ L 2 α (0, T ; H) so that if we define Recall that for any j ∈ N and any t ∈ [0, T ] we have that Since, for an arbitrary w ∈ H we have that (t − ·) α−1 2 w is in L 2 (0, t; H) , we can use (4.18) to obtain that lim j→∞ U kj (t), w = lim The statement above holds for any w ∈ H and all t ∈ [0, T ]. Thus, in H. However, this implies that u = u, as U kj converges to u in C([0, T ]; H). Therefore D α c u = v ∈ L 2 α (0, T ; H) and, by Lemma 4.6, we have the estimate for some constant C depending on α. As in the proof of Lemma 4.8 this implies that (4.6) holds. From this, we also see that the initial condition is attained in the required sense. It remains to show that the EVI (4.2) holds for u. From the construction of discrete solutions, one derives that for any w ∈ L 2 (0, T ; H) We will pass to the limit in this inequality. For the right hand side, it suffices to observe that U kj → u in C([0, T ]; H), V kj ⇀ v in L 2 (0, T ; H) and F kj → f in L 2 (0, T ; H). Thus, For the left hand side, the uniform convergence of U kj and the lower semicontinuity of Φ, give It remains to recall that D α c u = v ∈ L 2 (0, T ; H) to conclude that, according to Definition 4.2, u is an energy solution.
Theorem 4.5 also holds for this new definition. However, at least with our techniques, the requirements on the data u 0 ∈ D(Φ) and f ∈ L 2 α (0, T ; H) do not change.

Fractional gradient flows: Numerics
Since the existence of an energy solution was proved by a rather constructive approach, namely a fractional minimizing movements scheme, it makes sense to provide error analyses for this scheme. We will provide an a priori error estimate which, in light of the smoothness u ∈ C 0,α/2 ([0, T ]; H) proved in Theorem 4.5, is optimal. In addition, in the spirit of [30] we will provide an a posteriori error analysis.

5.1.
A priori error analysis. The a priori error estimate reads as follows. We comment that this result gives us a better rate compared to [28, Theorem 5.10].
Theorem 5.1 (a priori I). Let u be the energy solution of (1.2). Given a partition P, of maximal step size τ , let U ∈ H N be the discrete solution defined by (4.7) starting from U 0 ∈ D(Φ). Let U P and U P be defined as in (4.10) and (2.3), respectively. Then we have, where Φ 0 = max{Φ(U 0 ), Φ(u 0 )}, and the constant C depends only on α.
Proof. The proof can be obtained by following the same procedure employed in the proof of Lemma 4.9.
In the current situation, however, instead of comparing two discrete solutions we compare the exact and discrete ones. The only difference is that we allow U 0 = u 0 here, but this presents no essential difficulty. For brevity, we skip the details.

5.2.
A posteriori error analysis. Let us now provide an a posteriori error estimate between the discretization in (4.7) and the solution of (1.2). We will also show how, from this a posteriori error estimator, an a priori error estimate can be derived. Let us first introduce the a posteriori error estimator.
Definition 5.2 (error estimator). Let P be a partition of [0, T ] as in (2.2), and U ∈ H N denote the discrete solution given by (4.7). We define the error estimator function as Notice that the quantity E P (t) is nonnegative because . It is also, in principle, computable since it only depends on data, and the discrete solution U. It is then a suitable candidate for an a posteriori error estimator.
The derivation of an a posteriori error estimate begins with the observation that, for any w ∈ H, we have In other words, the function U P solves an EVI similar to (4.3) but with additional terms on the right hand side. We can then compare the EVIs by a now standard approach, that is, set w = u(t) in (5.4) and w = U P (t) in (4.3), respectively, to see that for almost every t ∈ [0, T ]. Consider the following notions of error: (t − s) α−1 σ(u(s); U P (s)) + σ(U P (s); u(s)) ds We have the following error estimate for E.
Theorem 5.3 (a posteriori). Let u be the energy solution of (1.2). Let P be a partition of [0, T ] defined as in (2.2) and let U ∈ H N be the discrete solution given by (4.7) starting from U 0 ∈ D(Φ). Let E and E P be defined in (5.6) and (5.3), respectively, The following a posteriori error estimate holds Proof. From (2.18) we infer ; u(t)) − σ(u(t); U P (t)).
The claimed a posteriori error estimate (5.7) follows from Lemma 2.8 by setting λ = 0, a(t) = ( U P − u)(t) , b(t) = 2 σ(U P (t); u(t)) + σ(u(t); U P (t)) , 5.3. Rate of convergence. Although we have already established an optimal a priori rate of convergence for our scheme in Theorem 5.1, in this section we study the sharpness of the a posteriori error estimator E P by obtaining the same convergence rates through it. We comment that neither in Theorem 5.1 nor in our discussion here, we require any relation between time steps. We will also consider some cases when the rate of convergence can be improved.

5.3.1.
Rate of convergence for energy solutions. Let us now use the estimator E P to derive a convergence rate or order O(τ α/2 ) for the error E, defined in (5.6), when f ∈ L 2 α (0, T ; H). Notice that such regularity a priori does not give any order of convergence for f − F P L 1 α (0,T ;H) in (5.7). Observe also that the rate that we obtain is consistent with classical gradient flow theories, where an order O(τ 1/2 ) is proved provided that u 0 ∈ D(Φ) and f ∈ L 2 (0, T ; H); see [30,Sec 3.2].
We first bound E P L 1 α (0,T ;H) . Theorem 5.4 (bound on E P L 1 α (0,T ;H) ). Under the assumption that U 0 ∈ D(Φ), the estimator E P , defined in (5.3), satisfies Proof. We bound the contributions E P,1 and E P,2 separately. The bound of E P,1 follows without change that of the term II of (4.15) in Lemma 4.9. Thus, To bound E P,2 , we recall the function Φ P , defined in Remark 4.7, and its properties. Define also Φ P (t) = Φ(U P (t)). We have On the one hand, proceeding as in the proof of Lemma 2.6 we obtain sup r∈[0,T ]ˆr 0 On the other hand, using we have for any r ∈ [0, T ] that Therefore combining the estimates for I 1 and I 2 we have proved that which together with (5.9) proves (5.8) because E P is nonnegative.
We next take advantage of Lemma 2.5, and derive a rate for E without additional smoothness assumptions on the right hand side f .
Theorem 5.5 (a priori II). Let u be the energy solution of (1.2). Let P be a partition of [0, T ] defined as in (2.2) and U ∈ H N be the discrete solution given by (4.7) starting from U 0 ∈ D(Φ). Let E be defined in (5.6). Then we have where the constant C depends only on α.

5.3.2.
Rate of convergence for smooth energies. Let us show that, at least for smoother energies, it is possible to obtain a better rate of convergence. We will, essentially, assume that the energy is locally C 1+β for β ∈ (0, 1]. More specifically in this section we consider energies that satisfy the following. There exists β ∈ (0, 1] such that for every R > 0, there is a constant C β,R > 0 for which where B R denotes the ball of radius R in H. Notice that, by Lemma 4.8, all the discrete solutions U P are uniformly bounded in C([0, T ]; H). Thus, we can fixR > 0 depending only on the data such that, for any partition P and all t ∈ [0, T ], U P (t) ∈ BR. Therefore, (5.11) implies that for some constant C β = C β,R .
A particular example to which this situation applies is the following. Let H = R d and Φ(w) = 1 p |w| p with p > 1. In this case, (5.12) holds with β = 1 for p ≥ 2 and β = p − 1 for p ∈ (1, 2). For p < 2, to reach β = 1, we must assume that u and U P stay uniformly away from zero. This example can, of course, be generalized. In this setting, we have the following improved estimate for E P L 1 α (0,T ;H) . Theorem 5.6 (improved bound). Assume that the energy Φ satisfies (5.12). Let u be the energy solution to (1.2), and denote by P a partition of [0, T ] defined as in (2.2). Denote by U P the solution of (4.7) starting from U 0 ∈ D(Φ). In this setting, the estimator E P defined in (5.3) satisfies for some constant C that depends on α, β, and the problem data.
Theorem 5.7 (improved rate: smooth energies). Assume that the energy Φ satisfies (5.12). Let u be the energy solution to (1.2), and denote by P a partition of [0, T ] defined as in (2.2). Denote by U P the solution of (4.7) starting from U 0 ∈ D(Φ). In this setting, if there is q > 1/α for which f ∈ W α(β+1)/2,q (0, T ; H) then the error E, defined in (5.6), satisfies where the constant C depends on α, β, q, T , and the problem data.

5.3.3.
Rate of convergence for linear problems. Let us now show how for certain classes of linear problems an improved rate of convergence can be obtained. We first assume that we have a Gelfand triple, V ֒→ H ֒→ V ′ and that where a : V × V → R is a nonnegative, symmetric, bounded, and semicoercive bilinear form. In this setting, (4.1) becomes D α c u, w + a(u, w) = f, w , ∀w ∈ V. Notice that the bilinear form induces an operator A : V → V ′ given by which implies that, for almost every t ∈ (0, T ), we have a problem in V ′ which reads We further assume that f ∈ L 2 α (0, T ; [·] V ). More essentially we also require u 0 ∈ D(∂Φ). The motivation for an improved rate of convergence is then the following, at this stage formal, calculation. From (2.18) we have which says that D α c u is uniformly bounded in L 2 α (0, T ; [·] V ). To make these considerations rigorous, we consider the discrete problem (4.7), which in this case reduces to (D α P U) n + AU n = F n , Then the computations can be followed verbatim to obtain that Similar to Lemma 2.4, we know that and hence D α c U P is uniformly bounded L 2 α (0, T ; [·] V ). With this additional regularity, we can obtain an improved rate of convergence. To see this, we will use that Φ is, essentially, quadratic to observe that in this case the error estimator, defined in (5.3) reduces to (5.17) These ingredients together give us the following improved estimate.
Theorem 5.8 (improved rate: linear problems). Assume that the energy Φ is given by (5.15), that the initial data satisfies Au 0 ∈ H, and that f ∈ L 2 α (0, T ; [·] V ). Let u be the energy solution to (1.2), and denote by P a partition of [0, T ] defined as in (2.2). Denote by U P the solution to (4.7) starting from U 0 ∈ H, such that AU 0 ∈ H. In this setting, we have that where the constant C depends only on α. This, immediately, implies that , so that if, in addition, we further have f ∈ W α,q (0, T ; H) for some q > 1/α, then where the constant C depends only on α, q and T .
Proof. Owing to Theorem 5.3 and equation (5.14), the convergence rate (5.19) follows directly from (5.18) in the same way as Theorem 5.7. We only need to prove (5.18) and bound E P L 1 α (0,T ;H) . Using (5.17), for every r ∈ (0, T ] we have Now, we invoke Lemma 2.6 with p = 2 and the semi-norm [·] V to obtain that . By (5.16), we have that D α c U P ∈ L 2 α (0, T ; [·] V ) uniformly in P and thus arrive at . This implies the desired bound for E P L 1 α (0,T ;H) and finishes the proof.

Lipschitz perturbations
In this section, inspired by the results of [3], we consider the analysis and approximation of a fractional gradient flow with a Lipschitz perturbation. Namely, we consider the following problem We assume that the perturbation function Ψ : (0, T ] × H → H satisfies 1. (Carathéodory) For every w ∈ H the mapping t → Ψ(t, w) is strongly measurable on (0, T ) with values in H. Moreover, there exists L > 0 such that for almost every t ∈ (0, T ) and every 2. (Integrability) There is w 0 ∈ L 2 α (0, T ; H) for which t → Ψ(t, w 0 (t)) ∈ L 2 α (0, T ; H). We immediately comment that our assumptions can fit the case where Φ is merely λ-convex. Moreover, these assumptions also guarantee the existence of ψ ∈ L 2 α (0, T ; R) for which Ψ(t, w) ≤ ψ(t) + L w , ∀w ∈ H.
Evidently, an energy solution to (6.1) satisfies, for almost every t ∈ (0, T ) and all w ∈ H, the EVI 6.1. Existence, uniqueness, and stability. Our main result in this direction is the following. Proof. We begin by proving existence. We essentially follow the idea used for the classical ODEs. A similar argument was also used in the proof of [25,Theorem 4.4].
For w ∈ L 2 α (0, T ; H) we denote by S(w) ∈ L 2 α (0, T ; H) the energy solution to D α c u(t) + ∂Φ(u(t)) ∋ f (t) − Ψ(t, w(t)), a.e. t ∈ (0, T ], u(0) = u 0 . Our assumptions and the results of Theorem 4.5 guarantee that this mapping is well defined, and moreover, S(w) ∈ L ∞ (0, T ; H). We want to show that there exists a fixed point w such that S(w) = w. If u i = S(w i ) for i = 1, 2, then for almost every t we have This readily implies that which as a consequence yields that, for every t ∈ [0, T ], . We claim that by induction, we can further obtain the following stability result for any t ∈ [0, T ] and positive integer n. In fact, for n = 1, we simply have Furthermore, if (6.3) holds for n = k, then for n = k + 1 which proves (6.3). Now consider w 0 ∈ L 2 α (0, T ; H) and the sequence of functions defined via w n = S n (w 0 ). It is easy to see that, for n ≥ 1, we have w n ∈ L ∞ (0, T ; H), and ∞ n=1 w n −w n+1 L ∞ (0,T ;H) converges because ∞ n=0 L n t αn Γ(αn + 1) = E α (Lt α ).
This shows that w n → u in L ∞ (0, T ; H) for some u. Since w n+1 = S(w n ), it follows immediately that u = S(u). This proves the existence of solutions. As for uniqueness, assume that we have two solutions u 1 and u 2 , for almost every t, we have Combining with the fact that u 1 (0) = u 2 (0) = u 0 , one obtains that u 1 (t) − u 2 (t) 2 = 0 for almost every t, which proves uniqueness. Finally, the estimate on the Caputo derivative trivially follows from the iteration scheme. We skip the details.
For diversity in our arguments, we present an alternative proof. The arguments here are inspired by those of [3, Theorem 5.1].
Alternative proof of Theorem 6.2. Let us, for µ > L 1/α , define which by the obvious inequalities e −µT ≤ e −µt ≤ 1, defines an equivalent norm in L 2 α (0, T ; H). Let S : L 2 α (0, T ; H) → L 2 α (0, T ; H) be as before. As shown, if u i = S(w i ) for i = 1, 2, then for every t we have which as a consequence yields that, for every r ∈ [0, T ], Obvious manipulations then yield which implies so that S is a contraction with respect to the norm · µ . We conclude then by invoking the contraction mapping principle. This unique fixed point, evidently, is a energy solution in the sense of Definition 6.1. Uniqueness and stability follow as before.
6.2. Discretization. Let us now present the numerical scheme for problem (6.1). We follow the previous notations and conventions regarding discretization so that, for any partition P of [0, T ] defined as in (2.2), we can also consider the discrete solution defined recursively via where F n is defined in (4.8) and Ψ n : H → H is defined by Clearly, for every n, Ψ n is Lipschitz continuous with Lipschitz constant L. Using the definition of D α P in (3.3) and K −1 P,nn = (K P,nn ) −1 = Γ(α + 1)τ −α n , we can rewrite (6.4) as Hence the discrete scheme can be recursively well-defined provided Lτ α < Γ(α + 1). For this reason, moving forward, we will implicitly operate under this assumption. It is possible to show that the discrete solutions in (6.4) satisfy (6.5) D α c U P L 2 α (0,T ;H) ≤ C, with a constant that depends on problem data but is independent of the partition P. To see this, we follow the arguments of either proof of Theorem 6.2, and realize that while the operator S may depend on P, the estimates that we obtain do not. 6.3. Error estimates. Let us now show how to derive error estimates for the problem with Lipschitz perturbation (6.1). We recall that the energy solution u to this problem satisfies (6.2). In addition, for simplicity, we will operate under the assumption that the perturbation does not depend explicitly on time, i.e., Ψ(t, w) = Ψ(w) for all w ∈ H. The general case only lengthens the discussion but brings nothing substantive to it, as the additional terms that appear can be controlled via arguments used to control terms of the form f (t) − F P (t).
Similar to the discussion before, we define the error estimator (6.6) E P,L (t) = E P (t) + Ψ(U P (t)), U P (t) − U P (t) , which, as before, is nonnegative. In addition, for any w ∈ H we have Setting w = u(t) in the inequality above and setting w = U (t) in (6.2) leads to for almost every t ∈ (0, T ). This implies the following error estimates. Theorem 6.3 (a posteriori: Lipschitz perturbations). Let u be the unique energy solution of (6.1).
Let P be a partition of [0, T ] defined as in (2.2) and let U ∈ H N be the discrete solution given by (6.4) starting from U 0 ∈ D(Φ). Let E and E P,L be defined in (5.6) and (6.6), respectively, The following a posteriori error estimate holds Proof. We argue as in the proof of (5.3). To make formulas shorter we omit the coercivity terms. From (2.18) and (6.7) we infer Then the error estimate (6.8) follows from Lemma 2.8 with We also comment here that by Lemma 2.6 where the constant C only depends on α. In addition, the norm on the right hand side is bounded independently of the partition P; see (6.5). Hence the convergence rates proved in Theorems 5.5 and 5.7 also hold for problems with a Lipschitz perturbation. Since the proofs are almost identical, we only state the theorems below without proofs. Theorem 6.4 (convergence rate: Lipschitz perturbations). Let u be the energy solution of (6.1). Let P be a partition of [0, T ] defined as in (2.2) and U ∈ H N be the discrete solution given by (6.4) starting from U 0 ∈ D(Φ). Let E be defined in (5.6). Then we have E ≤ u 0 − U 0 (E α (2LT α )) 1/2 + Cτ α/2 f L 2 α (0,T ;H) + D α c U P L 2 α (0,T ;H) , where the constant C depends only on α, L and T , but not on P. Theorem 6.5 (improved rate: smooth energies and Lipschitz perturbations). Assume that the energy Φ satisfies (5.12). Let u be the energy solution to (6.1), and denote by P a partition of [0, T ] defined as in (2.2). Denote by U P the solution of (6.4) starting from U 0 ∈ D(Φ). In this setting, if there is q > 1/α for which f ∈ W α(β+1)/2,q (0, T ; H) then the error E, defined in (5.6), satisfies where the constants C 1 and C 2 depend only on α, β, q, L, T , and the problem data, but are independent of P.
Finally we consider the setting of Section 5.3.3 with a Lipschitz perturbation. Similar to (6.5), we can show that D α c U P L 2 α (0,T ;[·]V ) is bounded uniformly with respect to the partition P. For this reason, an improved error estimate analogous to Theorem 5.8 can be proved in this case. Theorem 6.6 (improved rate: quadratic energies and Lipschitz perturbations). Assume that the energy Φ is given by (5.15), that the initial data satisfies Au 0 ∈ H, and that f ∈ L 2 α (0, T ; [·] V ). Let u be the energy solution to (6.1), and denote by P a partition of [0, T ] defined as in (2.2). Denote by U P the solution to (6.4) starting from U 0 ∈ H, such that AU 0 ∈ H. In this setting, we have that where the constant C depends only on α, L and T .

Numerical illustrations
In this section we present some simple numerical examples aimed at illustrating, and extending, our theory. All the computations were done with an in-house code that was written in MATLAB © . 7.1. Practical a posteriori estimators. We begin by commenting that, unlike the a posteriori estimators for the classical gradient flow proposed in [30], our a posteriori estimator E P is not constant on each subinterval of our partition P; see (5.3). Here we mention more computationally friendly alternatives, and their properties.
First, we define an estimator that is piecewise constant in time via D P (t) = max s∈[⌊t⌋P ,⌈t⌉P ] D α c U P (s) − F (s), U P (s) − U P (s) + Φ( U P (s)) − Φ(U P (s)) This is clearly an upper bound for E P (t).
One may also consider the simpler indicator (7.1) E P,n = (D α P U) n − F n , U n−1 − U n + Φ(U n−1 ) − Φ(U n ), n = 1, . . . , N. Although it is not always true that E P (t) ≤ E P,n(t) , this indicator is convenient to use in practice and gives reasonable results. In fact, this is the one that we implemented in the numerical examples of Section 7.3 below. 7.2. A linear one dimensional example. As a first simple example we consider the one dimensonal fractional ODE (7.2) D α c u + λu = 0, u(0) = 1, with λ > 0. From (2.19) we have u(t) = E α (−λt α ). This, obviously, fits our framework with H = R, and Φ(w) = λ 2 |w| 2 . Notice also that all the assumptions of Section 5.3.2 are also satisfied with β = 1. Thus, we expect a rate of order O(τ α ) when using (4.7) to approximate the solution over a uniform partition with time step τ .  Table 1. Convergence rate for the approximation of (7.2) using scheme (4.7) over a uniform partition of size τ . As predicted by Section 5.3.2, the rate is O(τ α ).  Table 1 shows, for λ = 0.001 and different values of α, the difference |u(1) − U N | which we use as a proxy for the error E H of (5.6). The rate of convergence is verified. 7.3. Adaptive time stepping. We now illustrate the use of the a posteriori error estimator E P given in (5.3) to drive the selection of the size of the time step. For a given tolerance ε we, at every step, choose the local time step τ n to guarantee that 2T α Γ(α + 1) where E P,n is given in (7.1). Then, by Theorem 5.3, we expect that u − U P L ∞ (0,T ;H) ≤ ε, provided the approximation error f − F P L 1 α (0,T ;H) is negligible. Notice that to drive the process we are using the simpler estimator E P ; see the discussion in Section 7.1.
Notice that this example fits the framework of Section 5.3.2 with β = p − 1. However, as mentioned there, it is not expected that the solution reaches zero in finite time, so we do not expect a reduced rate.
To compute the discrete solution, at every time step, we need to solve a nonlinear equation of the form U n + c U n |U n | p−2 − W n = 0, c = λτ α Γ(α + 1) where W n is known. We found the solution to this problem using Newton's method, which works for small values of τ . Table 2 shows the results for α = 0.5, p = 1.5, and λ = 1. These clearly indicate a rate of O(τ ).
At each time step one needs to solve a problem of the form U n + c ln(U n ) − W n = 0, c = λτ α Γ(α + 1) , and W n is known. This is solved with a Newton scheme, which runs into difficulties at the initial time step. We go around this issue by using as initial value for the iteration a very small positive number.  Table 4. Convergence rate for α = 0.5 and λ = 10 −6 in Example 3 of Section 7.4. The rate seems to be of order O(τ ), which is better than what the theory predicts. Table 3 presents the results for α = 0.5 and λ = 10 −6 . These indicate that the convergence rate is O(τ α ). Similar results for other choices of α and λ were obtained. 7.4.3. Example 3. As a final example we consider Φ(w) = −λ 1 − (1 − u) 2 + , u 0 = 0.
Acknowledgement AJS is partially supported by NSF grant DMS-1720213.