A Time-Fractional Borel–Pompeiu Formula and a Related Hypercomplex Operator Calculus

In this paper we develop a time-fractional operator calculus in fractional Clifford analysis. Initially we study the $$L_p$$Lp-integrability of the fundamental solutions of the multi-dimensional time-fractional diffusion operator and the associated time-fractional parabolic Dirac operator. Then we introduce the time-fractional analogues of the Teodorescu and Cauchy–Bitsadze operators in a cylindrical domain, and we investigate their main mapping properties. As a main result, we prove a time-fractional version of the Borel–Pompeiu formula based on a time-fractional Stokes’ formula. This tool in hand allows us to present a Hodge-type decomposition for the forward time-fractional parabolic Dirac operator with left Caputo fractional derivative in the time coordinate. The obtained results exhibit an interesting duality relation between forward and backward parabolic Dirac operators and Caputo and Riemann–Liouville time-fractional derivatives. We round off this paper by giving a direct application of the obtained results for solving time-fractional boundary value problems.


Introduction
Nowadays, one of the most studied fractional partial differential equation is the time-fractional diffusion equation due to its wide range of applications (see [8,17,22,23,25,30] and the references therein indicated). In the context of physics-mathematics this equation is connected with the non-Markovian diffusion processes with memory (see [23]), while in probability theory it is related to jumping processes (see [8]). In the classical case, the diffusion equation describes the heat propagation in a homogeneous medium. The time-fractional diffusion equation models the anomalous diffusions exhibiting sub-diffusive behavior, due to particle sticking and trapping phenomena (see e.g. [24]).
The multi-dimensional time-fractional diffusion equation case was studied in several papers, e.g. [3,4,10,12,13,[19][20][21]32]. In [32] the fundamental solution of this equation was deduced in terms of H-functions. In [12,13]  fractional diffusion-wave equations. In [3,4,[19][20][21] Luchko and his collaborators used the representation in terms of Mellin-Barnes type-integrals to study some properties of the multi-dimensional space-time-fractional diffusionwave equation. Explicit series representations for the fundamental solution of the diffusion-wave operator and the so-called time-fractional parabolic Dirac operator were obtained in [10], for arbitrary dimension.
During the last decades, Clifford analysis proved to be a good tool to study partial differential equations of mathematical-physics. In particular, we have the work of Gürlebeck and Sprößig based on a Borel-Pompeiu formula and on an orthogonal decomposition of the underlying function space where one of the components is the kernel of the corresponding Dirac operator [11]. This theory was successfully applied to a large type of equations, e.g., Lamé equations, Maxwell equations, and Navier-Stokes equations. Fractional versions of these equations has been attracting recent interest (cf. [5,18,27,31]).
The aim of this paper is to continue ideas introduced in [10] and use the fundamental solution of the time-fractional diffusion-wave and parabolic Dirac operators to develop, in the context of fractional Clifford analysis, a time-fractional operator calculus related to the time-fractional parabolic Dirac operator defined via left Caputo time-fractional derivative. In order to do that we initially study the L p -integrability of the fundamental solutions of the time-fractional diffusion and the time-fractional parabolic Dirac operators deduced in [10]. Then, we introduce the time-fractional analogues of the Teodorescu and Cauchy-Bitsadze operators and we investigate some important mapping properties. Moreover, we present a Hodge-type decomposition for the L p -space with respect to the time-fractional parabolic Dirac operator. The results exhibit an interesting "double duality" between forward and backward time-fractional parabolic Dirac operators, and between Caputo and Riemann-Liouville time-fractional derivatives. This double duality appears in a non-trivial generalization of the Stokes' formula as well as in the time-fractional Borel-Pompeiu formula and in the Hodge-type decomposition. Throughout the paper we show that we can always recover the results of the classical function theory for the parabolic Dirac operator when considering the limit case of β = 1. Possible applications of our fractional integrodifferential hypercomplex operator calculus are the study of boundary value problems with time-fractional derivatives such as the time-fractional Navier-Stokes equation and the time-fractional Scrödinger equation. For integer time derivatives these equations had been already study in the context of Clifford analysis (see [6,7]).
The structure of the paper reads as follows: in the Prelimaries's section we recall some basic facts about Clifford analysis, fractional derivatives and their main properties. In Section 3, we recall the fundamental solution of the n-dimensional time-fractional diffusion operator deduced in [10] and some estimates for this function and its derivatives deduced in [17]. In Section 4 we study the conditions that ensure the L p -integrability of the fundamental solution of the time-fractional diffusion operator and of the fundamental solution of the time-fractional parabolic Dirac operator. In the following section, we study the time-fractional Teodorescu and Cauchy-Bitsadze operators in a cylindrical domain. Finally, in Section 6 we present a Hodge-type decomposition for the L q -space, where one of the components is the kernel of the time-fractional parabolic Dirac operator. This represents the main result of the paper, apart from the time-fractional Borel-Pompeiu formula. We round off this paper by giving an immediate application to the resolution of time-fractional boundary value problems.

Hypercomplex analysis
We consider the n-dimensional vector space R n endowed with an orthonormal basis {e 1 , · · · , e n }. The universal real Clifford algebra C 0,n is defined as the 2 n -dimensional associative algebra which obeys the multiplication rule e i e j + e j e i = −2δ ij , i, j = 1, . . . , n. (1) A vector space basis for C 0,n is generated by the elements e 0 = 1 and e A = e h1 · · · e h k , where A = {h 1 , . . . , h k } ⊆ M = {1, . . . , n}, for 1 ≤ h 1 < · · · < h k ≤ n. Each element x ∈ C 0,n can be represented by The Clifford conjugation is defined by x = A x A e A , where e A = e h k · · · e h1 , and e j = −e j , for j = 1, . . . , n, and e 0 = e 0 = 1. We introduce the complexified Clifford algebra C n as the tensor product where the imaginary unit i of C commutes with the basis elements, i.e., ie j = e j i for all j = 1, . . . , n. To avoid ambiguities with the Clifford conjugation, we denote the complex conjugation by , in the sense that for a complex scalar w A = a A + ib A we have that w A = a A − ib A . The complex conjugation can be extended linearly to whole of the Clifford algebra and leaves the elements e j invariant, i.e., e j = e j for all j = 1, . . . , n.
A C n -valued function defined on an open set U ⊆ R n has the representation f = A f A e A with C-valued components f A . Properties such as continuity, differentiability, and integrability of a C n -valued function need to be understood componentwise. For instance f ∈ L p (U, C n ), or shortly f ∈ L p (U ), means that {f A } ⊂ L p (U ) or, equivalent, that U |f (x)| p dx < +∞. L 2 (U ) can be turned into a C n -module, with the following Clifford inner product Next, we introduce the Euclidean Dirac operator D x = n j=1 e j ∂ xj , which factorizes the n-dimensional Euclidean Laplacian, i.e., D 2 In order to define the parabolic Dirac operator we need to introduce a Witt basis. We start considering the embedding of R n into R n+2 and two new elements e + and e − such that e 2 + = +1, e 2 − = −1, and e + e − = −e − e + . Moreover, e + and e − anticommute with all the basis elements e j , j = 1, . . . , n. Hence, {e 1 , . . . , e n , e + , e − } spans R n+1,1 . With the elements e + and e − we construct two nilpotent elements f and f † given by These elements satisfy the following relations The extended basis {e 1 , . . . , e n , f, f † } allow us to define the parabolic Dirac operator as D x,t := D x + f∂ t + f † , where D x stands for the Dirac operator in R n . The operator D x,t acts on C n -valued functions defined on time dependent domains Ω × I ⊆ R n × R + , i.e., functions in the variables (x 1 , x 2 , . . . , x n , t) where x j ∈ R for j = 1, . . . , n, and t ∈ R + . For the sake of readability, we abbreviate the space-time tuple (x 1 , x 2 , . . . , x n , t) simply by (x, t), assigning x = x 1 e 1 + · · · + x n e n . For additional details on Clifford analysis, we refer the interested reader for instance to [7,9,11].

Fractional derivatives and special functions
Now we recall the definitions of the fractional integrals and fractional derivatives that will be used in the paper. Let a, b ∈ R with a < b and let β > 0. The left and right Riemann-Liouville fractional integrals I β a + and I β b − of order β are given by (see [15]) By RL D β a + and RL D β b − we denote the left and right Riemann-Liouville fractional derivatives of order β > 0 on [a, b] ⊂ R, which are defined by (see [15]) where m = [β] + 1 and [β] means the integer part of β. Let C D β a + and C D β b − denote, respectively, the left and right Caputo fractional derivative of order β > 0 on [a, b] ⊂ R, which are defined by (see [15]) Throughout the paper, AC m ([a, b]) denotes the class of functions g which are continuously differentiable on [a, b] up to the order m − 1 and g (m−1) is supposed to be absolutely continuous on [a, b]. Now, we recall an important result about the boundedness of the fractional integrals I β a + and I β b − (see Theorem 3.5 in [29]).
Theorem 2.1 If 0 < β < 1 and 1 < p < 1 β , then the operators I β a + and The previous theorem will play an important role in the study of the mapping properties of the time-fractional integral operators introduced in Section 5.2.
The fundamental solution presented in [10,17] is represented in terms of a Fox H-function. The Fox Hfunction H m,n p,q is defined via a Mellin-Barnes type integral in the form (see [16]) where a i , b j ∈ C, α i , β j ∈ R + , for i = 1, . . . , p and j = 1, . . . , q, and L is a suitable contour in the complex plane separating the poles of the two factors in the numerator. A detailed study of the properties and convergence of the Fox H-function can be seen in [16].

Estimates of the fundamental solution of the time-fractional diffusionwave equation and its derivatives
We consider the multi-dimensional time-fractional diffusion-wave equation defined by where (x, t) ∈ R n × R + , β ∈]0, 2[, and c > 0. Here, C ∂ β 0 + is the left Caputo time-fractional derivative of order β (see (9)), and ∆ x is the Laplace operator in R n . The first fundamental solution of (12) satisfying the initial conditions u(x, 0) = δ(x) if 0 < β < 1, and u(x, 0) = δ(x) and ∂ t u(x, 0) = 0 if 1 < β < 2 was deduced by several authors (see e.g. [10,15,19,20]). In [10] the authors obtained the fundamental solution G β n (x, t) in the form where H m,n p,q (z) is the Fox H-function. Asymptotic behaviour of this fundamental solution and its derivatives were studied in [17]. Therein, the authors studied the multi-dimensional space-time-fractional diffusion-wave equation given by where β ∈]0, 2[, α ∈]0, +∞[, and ∆ α x is the fractional Laplacian. We remark that we changed the roles of α and β in [17], in order to have the same fractional parameter in the time-fractional derivative. In [17] was introduced the following auxiliar function where H σ,γ corresponds to the following H-function Moreover, the function p(x, t) := p 0,0 (x, t) is the first fundamental solution of (14). From now on we consider c = 1 in (12) and (13), and α = 1 in (14), (15) and (16), which implies that p(x, t) = G β n (x, t). Therefore, the results presented in [17] for p(x, t) are the same for G β n (x, t). Since we want to prove the L p -integrability of G β n and its derivatives, we consider the following particular cases of Theorems 5.1 and 5.5 presented in [17] (note again the change of the roles of α and β in [17]).
and for |x| 2 t −β ≤ 1 Before we present the particular case of Theorem 5.5, we introduce, for each multi-index a = (a 1 , . . . , a n ) and k ∈ N, the following notation: (20) and for We end this section with a short remark about the notation used in the previous theorems and in the remaining parts of the paper (see [17] for more details). We write f g for |z| ≤ (resp. |z| ≥ ) if there exists a positive constant κ independent of z such that f (z) ≤ κ g(z) for |z| ≤ (resp. |z| ≥ ).

L p -integrability of the fundamental solution and its derivatives 4.1 The case of the time-fractional diffusion operator
In this section we prove the L p -integrability of G β n and its derivatives only for β ∈]0, 1[, which corresponds to the diffusion case. Let us start with the L p -integrability of G β n .

Theorem 4.1
The fundamental solution G β n belongs to L p (R n ×]0, T ]), T ∈ R, whenever p and β satisfy the following conditions: .
Proof: The proof is quite long and technical, however the calculations are straightforward. Considering x = r ω, with r > 0 and ω ∈ S n−1 , we obtain where A(S n−1 ) denotes the surface area of S n−1 . Let us start with the analysis of I, which corresponds to the integral over the region |x| 2 t −β ≤ 1. Taking into account the estimates (18) in Theorem 3.1 we need to consider separately the cases when n = 1, n = 2, and n > 2. Hence, by (18) we have • Case n = 1: • Case n = 2: Considering the change of variable s = r 2 t −β in the integral with respect to r, we obtain where Γ(a, z) denotes the incomplete gamma function (see [2]). The previous result is valid only under the conditions • Case n > 2: Making s = r 2 t −β in the integral with respect to r, we obtain under the conditions Now, let us study II, which corresponds to the integral over the region |x| 2 t −β ≥ 1. Taking into account (17) in Theorem 3.1 and applying the same change of variables s = r 2 t −β , we have, for any n ∈ N The previous result is only valid under the conditions Finally, combining the conditions (22), (23), (24), and (25) we obtain (i) Case n = 1: Taking into account (22) and (25), we get (ii) Case n = 2: Taking into account (23) and (25), we get (iii) Case n > 2: Taking into account (24) and (25), we conclude that p ∈ ]1, p 1 [, with

The case of the time-fractional parabolic Dirac operator
In [10] the authors studied the so-called time-fractional parabolic Dirac operator, which is a first-order differential operator that factorizes the time-fractional diffusion operator, and obtained its first fundamental solution. This operator is defined in the Clifford algebra setting by where D x = n j=1 e j ∂ xj is the Euclidean Dirac operator, C ∂ β 0 + is the left Caputo fractional partial derivative of order β ∈]0, 1[ given by (9), and f, f † are the Witt basis elements defined in (3). This operator satisfies the following factorization property (see [10]) Applying C D β x,0 + to G β n we obtain the representation of the fundamental solution of C D β x,0 + in terms of Hfunctions (see [10]) Taking into account (27), (19), and the fact that G β n (x, t) = p 0,0 (x, t), we get the following relation between C G β + and the derivatives of p 0,0 : where the time-fractional derivative was calculated via the following relation proved in [17] Making use of (28), Theorem 3.1, Theorem 3.2, and proceeding in a very similar way as it was done in the proof of Theorem 4.1, we obtain to the following result.
, T ∈ R, whenever p and β satisfy the following conditions: and β ∈ 0, 1 2 ; (ii) If n = 2 then p ∈ ]1, p 2 [, with Proof: The steps of the proof are similar to those in Theorem 4.1. We have Let us start with the analysis of I, which corresponds to the integral over the region |x| 2 t −β ≤ 1. This integral is split in three integrals: Taking into account (21) in Theorem 3.2 and applying the change of variables s = r 2 t −β with r = |x|, we obtain under the conditions p ∈ 1, min 2 + βn β(n + 1) , n n − 1 , β ∈]0, 1], and n > 1.
When n = 1 we have We pass now to the analysis of I 2 . From (18) in Theorem 3.1 we conclude that it is necessary to consider three cases depending on the value of n. Hence, applying the same change of variables already considered, s = r 2 t −β with r = |x|, we have • Case n = 1: under the conditions • Case n = 2: where Γ(a, z) denotes the incomplete gamma function and the parameters p and β satisfy the conditions • Case n > 2: (2 − β (p (n + 2) − n)) (n − p (n + 2)) , where p and β are such that p ∈ 1, min 2 + βn β(n + 2) , n n − 2 and β ∈]0, 1[.
Concerning I 3 it corresponds to the integral I in the proof of Theorem 4.1. Hence, the convergence conditions are (22), (23), and (24) when n = 1, n = 2, and n > 2, respectively. We pass now to the analysis of J, which corresponds to the integral over the region |x| 2 t −β ≥ 1. Once again this integral is split in three integrals: Taking into account (20) in Theorem 3.2 and applying the change of variables s = r 2 t −β with r = |x|, we obtain where E ν (z) is the exponential integral function (see [2]), and the parameters p and β satisfy the following conditions p ∈ 1, 2 + βn β(n + 1) and β ∈]0, 1].
It remains to study the L p -integrability in R n ×]0, T ], T ∈ R, of the derivatives of C G β + , which will be useful in the next section to establish the mapping properties of the time-fractional Teodorescu operator. Taking into account (19), (28), and (29) we can write From (39) we have the following theorem. (ii) If n = 2 then p ∈ ]1, p 4 [, where (iii) If n > 2 then p ∈ ]1, p 5 [, with The proof of this theorem follows the same line of reasoning of the proof of Theorem 4.2 and it is omitted. However, we remark that for the analysis of the first term of (39) we use Theorem 3.2, and for the analysis of the second and third terms we use Theorem 3.1. Concerning the L p -boundedness of (40) we have the following result.
We omit the proof of this theorem due to the similarities with the previous ones. However, we would like to point out that during the proof are used the estimates in Theorem 3.2 with k = 2. For the case when |x| 2 t −β ≥ 1 we use the estimate (20) in Theorem 3.2. For the case of |x| 2 t −β ≤ 1, and in order to guarantee the convergence of the involved integrals, we need to consider the following improvement of the estimate (21) in Theorem 3.2 The new estimate (41) is obtained applying the inequality |x j | ≤ |x| 2 instead of |x j | ≤ |x|, for x ∈ R n with n ≥ 2, in the proof of Theorem 5.5 in [17]. We remark that in the previous estimates we need to consider k = 2 since ∂ xj C G β + = ∂ xj D x p 0,0 is a second order derivative of p 0,0 . When n = 1 we can not use either the estimates (21) or (41). Therefore, we need to study the L p -integrability of ∂ x C G β + when n = 1 using the explicit expression of C G β + . By [10] we have that for n = 1 it is given by .
Since | C G β + | is an even function we can assume x > 0, which implies that x |x| = 1. By straightforward calculations we obtain Finally, studying the L p -integrability of (42) in R×]0, T ], T ∈ R, we obtain the conditions given in (i) of Theorem 4.4.

Time-fractional Stokes' theorem
In this section we develop a Stokes' Theorem for the time-fractional parabolic Dirac operators of Caputo and Riemann-Liouville types given by and where D x = n j=1 e j ∂x j is the Dirac operator in R n , and C ∂ β 0 + , C ∂ β T − , RL ∂ β 0 + and RL ∂ β T − are the Caputo and Riemann-Liouville time-fractional derivatives with parameter β ∈]0, 1[ given by (9), (10), (7), and (8). Due to the properties of the Witt basis elements and since D 2 x = −∆ x we have the following factorizations of the time-fractional diffusion operators From now on until the end of the paper we consider a cylindrical space time domain C = Ω×]0, T ] ⊂ R n ×(0, +∞).
In the next theorem we present a time-fractional Stokes' formula involving the operators RL D β x,T − and C D β x,0 + . Theorem 5.1 (Time-fractional Stokes's theorem) For u, v ∈ AC 1 (C)∩AC(C) the following time-fractional Stokes' formulas hold where dσ x,t = dσ x dt, dσ x = D x dx = n j=1 (−1) j+1 e j d x j is the oriented surface element, and dx = dx 1 . . . dx n represents the n-dimensional (oriented) volume element.
Proof: In order to not overload the paper we present only the proof of (47), however, we remark that the proof for (48), (49), and (50) is analogous. Suppose that u, v ∈ AC 1 (C) ∩ AC(C). We start deducing the Stokes' theorem for the operators C D β x,T − and RL D β x,0 + , without the f † -component. From (43), (44), and (10) we obtain For the integral I we apply the classical Stokes' formula for the Dirac operator obtaining where Γ 1 = [0, T ]×∂Ω. Concerning integral II, taking into account the definition of I 1−β T − (see (6)) and changing the order of integration we obtain Now, making w = t and applying integration by parts with respect to the time variable, we get dx.

Merging (52) and (53) in (51) leads to
Rearranging the terms in the previous identity we obtain Adding and subtracting u f † v in the right-hand side of (54) leads finally to the time-fractional Stokes' formula: We observe that the Stokes's formula for the time-dependent case in the classical Clifford analysis is given by (see [7]) where D ± x,t := D x + f∂ t ± f † are the forward/backward parabolic Dirac operators. However, in the fractional Clifford analysis setting we obtain a more complicatedly "double duality" relation. On the one hand the formula involves forward and backward time-fractional parabolic Dirac operators, and on the other hand it also involves both the Caputo and Riemann-Liouville derivatives. Moreover, in the fractional case the integral over ∂C is split into two integrals over Γ 1 and Γ 2 , which does not occur in the classical case. In the next section we deduce a time-fractional Borel-Pompeiu formula and we introduce the time-fractional analogous of the Teodorescu and Cauchy-Bitsadze operator. Before we do that we need to understand the behaviour of the time-fractional Dirac operator C D β y,T − when the arguments of the function u in (47) are translated and reflected in space and time. Denoting the translation operator by T θ1,θ2 u(y, w) := u(θ 1 +y, θ 2 +w) and the reflection operator by R y,w u(y, w) := u(−y, −w), and taking into account the definitions of the timefractional parabolic Dirac operators presented in (43), and the definitions of the right and left Caputo fractional derivatives presented in (10) and (9), we can deduce by straightforward calculations the following relation (where the derivative is with respect to the variable (y, w) ∈ C):

Time-fractional Teodorescu and Cauchy-Bitsadze operators
In this section we introduce the time-fractional Teodorescu and Cauchy-Bitsadze operators and study some of their properties. Replacing u by C G β + in (47) leads to the following time-fractional Borel-Pompeiu formula and also to the time-fractional Cauchy's integral formula.
Thus, replacing in (47) u by C G β + (x − y, t − w) = T x,t R y,w C G β + (y, w) and using (57) with θ 1 = x and θ 2 = t, we obtain which leads to the time-fractional Borel-Pompeiu formula (58). Additionally, if v is in the kernel of RL D β y,0 + then the second integral of the left-hand side of the preceding expression is equal to zero. Therefore, we arrive at the time-fractional Cauchy's integral formula stated in (59).
We observe that, as a consequence of the "double duality" mentioned previously, we can also deduce another three alternative versions of the time-fractional Borel-Pompeiu formula from (48), (49), and (50).

Remark 5.4
When β = 1 the time-fractional Borel-Pompeiu formula and the time-fractional Cauchy's integral formula reduce to the corresponding formulae presented in [7].
Based on (58) we introduce the time-fractional Teodorescu and Cauchy-Bitsadze operators and study their properties.
Definition 5.5 Let g ∈ AC 1 (C). Then the linear integral operators of convolution type are called the time-fractional Teodorescu and Cauchy-Bitsadze operators, respectively.

Remark 5.6
In the case β = 1, the time-fractional operators C T β + and C F β + coincide with the correspondent ones defined in [7].
The previous definition allows us to rewrite (58) in the alternative form Before we deduce two operating properties of the fractional integral operators (60) and (61), we need to understand the behaviour of the time-fractional Dirac operator C D β x,T − when the arguments of the function over which we apply the derivatives is only translated. Taking into account the definition of C D β x,T − and the definition of the left Caputo fractional derivative presented in (9), we can deduce by straightforward calculations the following relation (where the derivative is with respect to the variable (x, t) ∈ C): Theorem 5.7 The time-fractional operator C T β + is the right inverse of C D β x,0 + , i.e., for g ∈ L p (Ω), with p ∈ 1, 1 1−β and β ∈]0, 1], we have Proof: In view of the definition of C T β + given in (60), and the relation deduced in (62) with θ 1 = −y, θ 2 = −w, and the fact that C G β + is the fundamental solution of C D β x,0 + satisfying C D β In a similar way as in [14], we introduce the fractional Sobolev space W 1,β,p a + (C), specifically adapted to our problem and considering the left Caputo fractional derivative (9), with the norm · W 1,β,p a + (C) given by: where · Lp(C) is the usual L p -norm in C and β ∈]0, 1].
Proof: In view of the definition of C F β + given in (61), and the relation deduced in (62) with θ 1 = −y, θ 2 = −w, and the fact that C G β + is the fundamental solution of C D β x,0 + satisfying C D β where the integral over Γ 1 is equal to 0 because (x, t) ∈ C and (y, w) ∈ Γ 1 ∪ Γ 2 and, therefore, (x − y, t − w) = (0, 0). Now, we present some mapping properties of the fractional operators C T β + and C F β + .
Theorem 5.9 The operator C T β + is bounded from L q (C) to L r (C) with 1 + 1 r = 1 p + 1 q , such that q ∈ 1, 1 1−β and the parameters p and β are in the conditions of Theorem 4.2.

Proof:
Let q ∈ 1, 1 1−β , and the parameters p and β be in the conditions of Theorem 4.2. Using the Young's inequality for convolutions (see Theorem 1.4 in [29]) and taking into account the L p −integrability of C G β + studied in Theorem 4.2, we obtain Concerning the derivatives of C T β + we have the following theorem.
Theorem 5.10 Let g ∈ L q (C) with q = 1, 1 1−β . The spatial derivatives of C T β + with respect to x j , j = 1, 2, . . . , n and the time-fractional derivative of C T β + are bounded and satisfy the mapping properties Proof: To show the boundedness of the operator C T β + , it suffices to study the convolution terms (see [26]) where on the right-hand side of the last equality we used (62) with θ 1 = y and θ 2 = w. In Theorems 4.4 and 4.3 were studied the conditions over p and β such that the kernels of the previous convolutions are L p -functions. Hence, making use of the Young's inequality for convolutions (see Theorem 1.4 in [29]), we get for the partial derivatives with respect to x j g Lq(C) , j = 1, . . . , n, and for the time-fractional derivative where 1 + 1 r = 1 p + 1 q , with q ∈ 1, 1 1−β , and the parameters p and β are in the conditions of Theorems 4.4 and 4.3, respectively.
The preceding last two theorems allow us to prove the continuity of C T β + .
This result is a direct consequence of Theorem 5.9 and Theorem 5.10 and, therefore, we omit the proof. We observe that when β = 1 the mapping results for the Teodorescu operator correspond to the classical ones (see [7]). Now, we study the mapping properties of C F β + .

Remark 5.13
When β = 1 we obtain the classical results obtained in [7] for the time-dependent T and F operators.

Hodge-type decomposition and Boundary Value Problems
The aim of this section is to obtain a Hodge-type decomposition and to present an immediate application of this decomposition in the resolution of boundary value problems involving the time-fractional diffusion operator.
Theorem 6.1 Let q ∈ 1, p 1−(1−β)p , p ∈ 1, 1 1−β and β ∈]0, 1]. The space L q (C) admits the following direct decomposition Proof: we denote the unique operator solution for the problem (see [1,28] where the existence and uniqueness of solutions of this type of boundary value problems is studied) where u, v ∈ L p (C). As a first step we take a look at the intersection of the two spaces that appear in the decomposition. Let We directly see that C D β x,0 + h = 0 in C. Moreover, since h ∈ C D β we obtain that g = 0. Consequently, h = 0. Hence, the intersection of these subspaces only contains the zero function, which implies that the sum is direct. Now, let h ∈ L q (C) and h 2 such that .
Therefore, h 1 ∈ ker C D β x,0 + . Since h ∈ L p (C) was arbitrarily chosen our decomposition is a direct decomposition of the space L q (C). Corollary 6.2 From the decomposition (65) we have the following projectors Remark 6.3 As a consequence of the "double duality" mentioned previously we can also deduce another Hodgetype decompositions and P and Q-type projectors for the operators C D β x,T − , RL D β x,0 + and RL D β x,T − .
As in the previous results, when β = 1 the Hodge-type decomposition presented in Theorem 6.1 coincides with the decomposition presented in [7]. Moreover, for β = 1 we have that q = p and p ∈]1, +∞[. For the particular case of p = 2 the decomposition is orthogonal with respect to the inner product (2) (see Theorem 3.3 in [7]). We end this section by presenting an application of our results solving boundary value problems.
Theorem 6.4 Let g ∈ W 1,β,p 0 + (C) with p ∈ 1, 1 1−β and β ∈]0, 1]. Then, the solution of the problem is given by h = C T β The proof is based on the properties of the operator C T β + and of the projector C Q β + . Since C T β + is the right inverse of C D β x,0 + , we get Corollary 6.5 Let g ∈ W 1+ 1 p ,β+ 1 p ,p 0 + (∂C) with p ∈ 1, 1 1−β and β ∈]0, 1]. Then, the solution of the problem is given by where g is the W 2,β+1,p 0 + -extension of g.
Proof: Let p and β in the conditions given. Since g ∈ W 1+ 1 p ,β+ 1 p ,p 0 + (∂C) there exists a W 2,β+1,p 0 + -extension g with tr( g) = g. If we consider h = f + g, then the problem (67) will be transformed into From Theorem 6.4 we conclude that, the solution f of the previous problem has the form Using the time-fractional Borel-Pompeiu formula (58), and taking into account that C P β + = I − C Q β + and C D β Since h = f + g we get (68).
From Theorem 6.4 and Corollary 6.5 we get Theorem 6.6 Let f ∈ W 1,β,p 0 + (C) and g ∈ W is given by where g is the W 1+p,β+p,p 0 + -extension of g.
Proof: Let h 1 and h 2 be the solutions of the problems (66) and (67), respectively. Then h = h 1 + h 2 solves the boundary value problem (69).
We observe that as consequence of the "double duality" mentioned across the paper, the previous results can be also deduced with the correspondent rearrangements in the definition of the Teodorescu operators and Q projectors, according to the correspondent Borel-Pompeiu formula and Hodge-type decomposition, for the

Conclusions
In this work we presented a generalization of several results studied in Clifford analysis in the context of nonstationary problems (see e.g., [7]) to the context of fractional Clifford analysis. Possible applications of our fractional integro-differential hypercomplex operator calculus are the study of boundary value problems with time-fractional derivatives such as the time-fractional Navier-Stokes equation and the time-fractional Scrödinger equation. This study is out of the scope of this paper and it will be subject of future research.