Strong convergence of an fractional exponential integrator scheme for the finite element discretization of time-fractional SPDE driven by standard and fractional Brownian motions

The aim of this work is to provide the first strong convergence result of numerical approximation of a general time-fractional second order stochastic partial differential equation involving a Caputo derivative in time of order $\alpha\in(\frac 12; 1)$ and driven simultaneously by a multiplicative standard Brownian motion and additive fBm with Hurst parameter $H\in(\frac 12, 1)$, more realistic to model the random effects on transport of particles in medium with thermal memory. We prove the existence and uniqueness results and perform the spatial discretization using the finite element and the temporal discretization using a fractional exponential integrator scheme. We provide the temporal and spatial convergence proofs for our fully discrete scheme and the result shows that the convergence orders depend on the regularity of the initial data, the power of the fractional derivative, and the Hurst parameter $H$.


Introduction
We analyse the strong numerical approximation of the following time fractional SPDE with initial value of the following type On the Hilbert space H = L 2 (Λ), where Λ ⊂ R d , d = 1, 2, 3 is bounded and has smooth boundary, T > 0 is the final time, A is a linear operator which is unbounded, not necessarily self-adjoint and is assumed to generate a semigroup S(t) := e −tA , C ∂ α t is the Caputo fractional derivative with α ∈ ( 1 2 , 1) and I 1−α t is the fractional integral operator which will be given in the next section. The . It is well known [2,22] that the noises can be represented as follows where q i , e i , i ∈ N are respectively the eigenvalues and eigenfunctions of the covariance operator Q, q 1 i , e 1 i , i ∈ N are respectively the eigenvalues and eigenfunctions of the covariance operator Q 1 , β i are mutually independent and identically distributed standard normal distributions and β H i are mutually independent and identically distributed fractional Brownian motion (fBm). The noises W and B H are supposed to be independent. Precise assumptions on the nonlinear mappings G and Φ to ensure the existence of the mild solution of (1) will be given in the following section.
Equation of type (1) with Φ = 0 might be used to model random effects on transport of particles in medium with thermal memory [31]. So due to the self-similar and long-range dependence properties of the fBm, when modelling such phenomena, it is recommended to incorporate the fBm process in order to obtain a more realistic model. During the last few decades, the theory of fractional partial differential equations has gained considerable interest over time. From the point of view of computations, several numerical methods have been proposed for solving time fractional partial differential equations (for details, see [13,15,6,25,8,5] and the reference therein). Note that the time stepping methods used in all the works mentioned until now are based on finite difference methods. However theses schemes are explicit, but unstable, unless the time stepsize is very small.
To solve that drawback, numerical method based on exponential integrators of Adams type have been proposed in [9]. The price to pay is the computation of Mittag-Leffler (ML) matrix functions.
As ML matrix function is the generalized form of the exponential of matrix function, works in [10,19,24] have extended some exponential computational techniques to ML. Note that up to now all the numerical algorithms presented are for time fractional deterministic PDEs with self adjoint linear operators.
Actually, few works have been done for numerical methods for Gaussian noise for time fractional stochastic partial differential equation (see [11,30,31,29]). Note that this above works have been done for self adjoint linear operator, so numerical study for (1) with Φ = 0 and non self adjoint operator is still an open problem in the field, to the best of our knowledge. However, it is important to mention that if H = 1 2 the process B H is not a semi-martingale and the standard stochastic calculus techniques are therefore obsolete while studying SPDEs of type (1). Alternative approaches to the standard Itô calculus are therefore required in order to build a stochastic calculus framework for such fBm. In recent years, there have been various developments of stochastic calculus and stochastic differential equations with respect to the fBm especially for H ∈ ( 1 2 , 1) (see, for example [1,2,18]) and theory of SPDEs driven by fractional Brownian motion has been also studied. For example, linear and semilinear stochastic equations in a Hilbert space with an infinite dimensional fractional Brownian motion are considered in [3,4]. However numerical scheme for time fractional SPDEs (1) driven both by fractional Brownian motion and standard Brownian motion have been lacked in the scientific literature to the best of our knowledge. Our goal in this work is to build the first numerical method to approximate the time fractional stochastic partial differential equation (1) driven simultaneously by a multiplicative standard Brownian motion and an additive fractional Brownian motion with Hurst parameter H ∈ ( 1 2 , 1) using finite element method for spatial approximation and a fractional version of exponential [16,21] Euler scheme for time discretization. Since the ML function is more challenging than the exponential function and our main result is based on novel preliminary results on ML functions. The analysis is complicated and is very different to that of a standard exponential integrator scheme [20](where α = 1) since the fractional derivative is not local and therefore numerical solution at given time depends to all previous numerical solutions up to that time. This is in contrast to the standard exponential numerical scheme where the numerical solution at a given time depends only of that of the previous nearest solution. We provided the strong convergence of our fully discrete scheme for (1). Our strong convergence results examine how the convergence orders depend on the regularity of the initial data, the power of the fractional derivative, and the Hurst parameter.
The rest of the paper is structured as follows. In Section 2, Mathematical settings for standard and fractional calculus are presented, along with the existence, uniqueness, and regularities results of the mild solution of SPDE (1). In Section 3, numerical schemes for SPDE (1) are presented, we discuss about space and time regularity of the mild solution X(t) of (1) given by (21). The spatial error analysis is done in Section 4. We end the paper in Section 5, by presenting the strong convergence proof of the our novel numerical scheme.

Mathematical setting, main assumptions and well posedness problem
In this section, we review briefly some useful results on standard and fractional calculus, introduce notations, definitions, preliminaries results which will be needed throughout this paper and the proof of existence and uniqueness of the mild solution of (1). • {B H (t), t ∈ [0, T ]} has continuous sample paths P a.s., where Cov(X, Y ) denotes the covariance operator for the Gaussian random variables X and Y and E stands for the mathematical expectation on (Ω, F , P, {F t } t≥0 ).
Notice that if H = 1 2 , the fractional Brownian motion coincides with the standard Brownian motion. Throughout this paper the Hurst parameter H is assumed to belong to ( 1 2 , 1).
Let (K, ., . K , . ) be a separable Hilbert space. For p ≥ 2 and for a Banach space U, we denote by L p (Ω, U) the Banach space of p-integrable U-valued random variables. We denote by L(U, K) the space of bounded linear mapping from U to K endowed with the usual operator norm . L(U,K) and L 2 (U, K) = HS(U, K) the space of Hilbert-Schmidt operators from U to K equipped with the following norm In the rest of this paper to simplify the presentation, we assume the SPDE (1) to be second order of the following type.
where f, g : Λ × R → R is globally Lipschitz continuous, φ : R → R is bounded. In the abstract framework (12), the linear operator A takes the form . We assume that there exists a positive constant c 1 > 0 such The functions F : H → H, G : H → L 0 2 and Φ ∈ L 0 2 are defined by for all x ∈ Λ, v ∈ H, u, w ∈ Q 1/2 (H). As in [7,16], we introduce two spaces H, and V such that H ⊂ V ; the two spaces depend on the boundary conditions of Λ and the domain of the operator A.
For Dirichlet (or first-type) boundary conditions, we take For Robin (third-type) boundary condition and Neumann (second-type) boundary condition, which is a special case of Robin boundary condition, we take V = H 1 (Ω) where ∂v/∂v A is the normal derivative of v and v A is the exterior pointing normal n = (n i ) to the boundary of A given by Using Gårding's inequality (see e.g. [26]), it holds that there exist two constants c 0 and λ 0 > 0 such that the bilinear form a(., .) associated to A satisfies By adding and subtracting c 0 Xdt in both sides of (12), we have a new linear operator still denoted by A, and the corresponding bilinear form is also still denoted by a. Therefore, the following coercivity property holds Note that we have create a new linear term −c 0 X in the right-hand side of (12). Thus we obtain a new equivalent form to (12) as we rewrite it in its contracted form as follows Note that the expression of nonlinear term F has changed as we included the term −c 0 X in the new nonlinear term that we still denote by F . The coercivity property (14) implies that A is the infinitesimal generator of a contraction semigroup S(t) = e −tA on L 2 (Λ) [7]. Note that this is due to the fact that the real part of the eigenvalues of A is positive. Note also that the coercivity property (14) also implies that A is a positive operator and its fractional powers are well defined and for any where Γ(α) is the Gamma function.
Remark 2.1. As we have mentioned, the linear operator A is the infinitesimal generator of a contraction semigroup S(t) = e −tA , and therefore Let be the fractional semigroups from Proposition 1 and (18), we also have Note that (18)- (20) are also hold is S(t) and S 1 (t) are replaced respectively by their semi discrete form S h (t) and S 1h (t) obtained after the finite element method. In the sequel of this paper, the contraction propriety (18)- (20) for S 1h (t) will play a key role in the analysis our numerical scheme.
Thanks to (11), the fractional semigroup operators S 1 (t) and S 2 (t) can be rewritten as and we obtain the following properties of that fractional semigroup (S i (t)) t∈(0,T ) , i = 1, 2. and Proof. See [21,Lemma 4] for the proof of (24), (27) and (26). The proof of (25) is similar to that of [21, (29)], we have hence Remark 2. Lemma 2 also holds with a uniform constant C (independent of h) when A and S i (t), i = 1, 2 are replaced respectively by their discrete versions A h and S ih (t) defined in Section 4.
In order to ensure the existence and the uniqueness of mild solution for SPDE (1) and for the purpose of convergence analysis we make the following assumptions.
Assumption 2 (Non linearity term F ). We assume the non-linear mapping F : H → H, to be linear growth and Lipschitz continuous ie, for κ ∈ [0, 1] there exists positive constant L > 0 such that Assumption 3 (Standard Noise term). We assume that the diffusion coefficient G : H → L 0 2 satisfies the global Lipschitz condition and the linear growth ie, for τ ∈ [0, 1], there exists a positive constant L > 0 such that Assumption 4 (Fractional Noise term). The deterministic mapping Φ : where β is defined as in Assumption 1.
We are no to present the following result of existence and uniqueness of mild solution of SPDE (1).
where Proof. We define the operator ξ : In order to obtain our result, we use the Banach fixed point to prove that the mapping ξ has a unique fixed point in L 2 (Ω × [0, T ], H). The proof will be splitted into two steps.
To see this, let x, y ∈ L 2 (Ω × [0, T ], H), then from (32) we get Cauchy-Schwarz inequality, the stability property of fractional semigroup S 2 (t) (20) and Assumption 2 yield Using the contraction argument of the semigroup, Ito isometry (6) and Assumption 3, we have Hence putting (39) and (40) in (38) holds The fact that 2L  (1), and x(t) is asymptotic stable in mean square. The proof is thus completed.
In all that follows, C denotes a positive constant that may change from line to line. In the Banach space D(A γ 2 ), γ ∈ R, we use the notation A γ 2 · = · γ and we now present the following regularity results.

Regularity of the mild solution
We discuss the space and time regularity of the mild solution X(t) of (1) given by (21) in this section. The following theorem presents the spatial and time regularity result.
and for 0 ≤ t 1 < t 2 ≤ T , the following optimal time regularity hold Moreover (42)  , taking the squared-norm and We bound II i , i = 1, 2, 3 one by one.

Space approximation and error estimates
We consider the discretization of the spatial domain by a finite element triangulation with maximal length h satisfying the usual regularity assumptions. Let V h ⊂ H denote the space of continuous functions that are piecewise linear over triangulation J h . To discretise in space, we introduce P h from L 2 (Ω) to V h define for u ∈ L 2 (Ω) by The discrete operator where a is the corresponding bilinear form of A. Like the operator A, the discrete operator A h is also the generator of a contraction semigroup S h (t) := e −tA h . The semidiscrete space version of Note that A h , P h G and P h Φ satisfy the same assumptions as A, G and Φ respectively. The mild solution of (60) can be represented as follows Where S 1h and S 2h are the semidiscrete version of S 1 and S 2 respectively defined by (22) and (23).
Let us define the error operators Then we have the following Lemma.
(ii) Let 0 ≤ γ ≤ 1, then there exists a constant C such that The following lemma provides an estimate in mean square sense for the error between the solution of SPDE (16) and the spatially semidiscrete approximation (61).
Lemma 4 (Space error). Let X and X h be the mild solution of (16) and (60), respectively. Let Assumptions 1 -4 be fulfilled then there exits a constant C independent of h, such that Proof. Define e(t) := X(t) − X h (t). By (21) and (61), taking the norm, using triangle inequality We will analyse the above terms IV i , i = 1, 2, 3, 4 one by one.
For the first term IV 1 , using Lemma 3 with r = γ = 2H + β − 1 and Assumption 1 yields For the second term IV 2 , by adding and subtracting a term, applying triangle inequality and the estimate (a + b) 2 ≤ 2a 2 + 2b 2 , we split it in two terms as follows Firstly, adding and subtracting a term, applying Cauchy-Schwartz inequality, Lemma 3 (i) with r = 2H + β − 1, ρ = 0, Lemma (2) more precisely (43) for the first term and Lemma 3 (ii) with γ = 0, Lemma 1 and Assumption 2 with κ = 0 yields Secondly, applying the Cauchy-Schwartz inequality, boundedness of P h , S 2h (t) and Assumption 2, it holds Putting (68) and (69) in (67), we obtain For the third term IV 3 , by adding and subtracting a term, using the triangle inequality and the estimate (a + b) 2 ≤ 2a 2 + 2b 2 , we have .

Fully discrete Euler scheme and its error estimate
In this section, we consider a fully discrete approximation scheme of SPDE (16), more precisely a fractional exponential scheme.

Fractional exponential scheme
Let ∆t the time step size and M ∈ N such that T = M∆t, hence for all m ∈ {0, 1, ..., M}, X h m or Y h m denotes the numerical approximation of X h (t m ) with t m = m∆t. Since we are not longer dealing with a semi group, the fractional exponential integrator scheme to (60) extending the standard exponential in [20]) is giving by.
As we can observe, the analysis will be more complicated and different to that of a standard exponential integrator scheme [20](where α = 1) since the fractional derivative is not local and therefore numerical solution at t m depends to all previous numerical solutions up to t m . This is in contrast to the standard exponential numerical scheme where the numerical solution at t m depends only of that at t m−1 .
Our main result, which is indeed our full convergence result is given in the following theorem.
Theorem 3. Let X(t m ) be the mild solution of (16) at time t m = m∆t given by (21). Let X h m be the numerical approximation through (73). Under Assumptions 1-4, the following estimation holds Before giving the proof, let us present a preparatory result.
where C is a positive constant independent of h.

Proof of Theorem 3
We are now ready to prove the first result of our main theorem. In fact using the standard technique in the error analysis, we split the fully discrete error in two terms as Note that the space error err 0 is estimate by Lemma 4. It remains to estimate the time error err 1 .
We recall that given the mild solution at time t m = m∆t of the semidiscrete problem (60) is given by Decomposing (77), we have and thanks to (73), we rewrite the numerical solution X h m in the integral form as where the notation s is defined as in [20, (89)] by s := t ∆t ∆t.
Subtracting (78) and (79), applying triangle inequality, taking the square and the estimate (a + b + c) 2 ≤ 3(a 2 + b 2 + c 2 ), it follows that By adding and subtracting a term, using triangle inequality, we recast V 2 as follows Using the martingale property of the stochastic integral, the Itô isometry (6), inserting an appropriate power of A h , (75), the semidiscrete version of (24) and (27) with σ = η = δ = 2H+β−1 2 and t 1 = t m − s, t 2 = t m − t j , Assumption 3 with τ = 2H+β−1 2 and the semidiscrete version of Theorem 2(more precisely (42)) yields To estimate the last two terms V 2 22 and V 2 23 , using the martingale property of the stochastic integral, applying Itô isometry, boundedness of P h and S 1h , Assumption 3 and the semidiscrete version of Theorem 2 (more precisely (43)) holds