Yavar Kian
email: yavar.kian@univ-amu.fr

Masahiro Yamamoto

RECONSTRUCTION AND STABLE RECOVERY OF SOURCE TERMS APPEARING IN DIFFUSION EQUATIONS

Keywords: Inverse source problems, fractional diffusion equation, reconstruction, well-posedness, stability estimate Mathematics subject classification 2010 : 35R30, 35R11

come

Reconstruction and stable recovery of source terms appearing in diffusion equations

∂ xi a ij (x)∂ xj u , x ∈ Ω, where a ij = a ji ∈ C 2 (ω × [-, ]), 1 ≤ i, j ≤ d, satisfy d i,j=1
a ij (x)ξ i ξ j c|ξ| 2 , x ∈ Ω, ξ = (ξ 1 , . . . , ξ d ) ∈ R d , a jd (x , ± ) = 0, a dd (x , ± ) > 0, x ∈ ω, j = 1, . . . , d -1.

(1.1)

For α ∈ (0, 2) \ {1}, we denote by ∂ α t the Caputo fractional derivative with respect to t given by

∂ α t u(t, x) := 1 Γ([α] + 1 -α) t 0 (t -s) [α]-α ∂ [α]+1 s u(s, x)ds, (t, x) ∈ Q,
and by ∂ 1 t the usual derivative in t. We set

m α =    0, 0 < α 1, 1, 1 < α < 2, 1
T ∈ (0, +∞), Σ = (0, T ) × ∂Ω, 0 < α < 2 and we consider the following problem    ∂ α t u + Au = F (t, x), (t, x , x d ) = (t, x) ∈ Q, ∂ k t u(0, x) = 0, x ∈ Ω, k = 0, . . . , m α .

(1.2)

We associate with this problem the following boundary conditions ∂ m x d u(t, x , ± ) = 0, (t, x ) ∈ (0, T ) × ω, m = 0, 1,

∂ n ν u(t, x) = 0, (t, x) ∈ (0, T ) × ∂ω × (-, ), n = 0, 1.

(1.4)

In the same way, we consider the problem

   ∂ α t u + Au + q(t, x )u = F (t, x), (t, x , x d ) = (t, x) ∈ Q, ∂ k t u(0, x) = 0, x ∈ Ω, k = 0, m α .
(1.5)

∂ x d u(t, x , 0) = u(t, x , ) = 0, (t, x ) ∈ (0, T ) × ω, (1.6)
u(t, x) = 0, (t, x) ∈ (0, T ) × ∂ω × (0, ).

(

In the present paper, we treat the inverse problem of determining the source term F or the coefficient q from measurements of the solution of (1.2)-(1.4) or (1.5)-(1.7) on a subboundary of the cylindrical domain Q or Q.

1.2. Obstruction against the uniqueness. We recall that there is an obstruction against the recovery of general source terms F from any type of measurements of the solution of (1.2)-(1.4) (resp. (1.5)-

(1.7)) on (0, T ) × ∂Ω (resp. (0, T ) × ∂ Ω). Indeed, choose χ = 0, ∈ C ∞ 0 (Q) and consider F := ∂ α t χ + Aχ. From the uniqueness of the weak solution of problem (1.2)-(1.4) (see Section 1. [START_REF] Cheng | Uniqueness in an inverse problem for a one dimensional fractional diffusion equation[END_REF] for more details and see Definition 1.1 below for the definition of weak solutions), we knows that u = χ and, since χ = 0, we deduce that F = 0. However, we have ∂ k ν u(t, x) = ∂ k ν χ(t, x) = 0, (t, x) ∈ (0, T ) × ∂Ω, k = 0, 1, . . . , with ν the outward unit normal vector of ∂Ω.

Facing this obstruction against the uniqueness, we will consider source terms of the form

F (t, x , x d ) := f (t, x )R(t, x , x d ), (t, x , x d ) ∈ (0, T ) × ω × (-, ) (1.8)
and, assuming that R is known, we will consider the problem of determining f . 1.3. Motivations. We recall that the problem (1.2)-(1.4) (resp. (1.5)-(1.7)) is associated with different models of diffusion. The non-integer value of the power α, is frequently used for describing anomalous diffusion derived from continuous-time random walk models (see e.g., [START_REF] Metzler | The random walk's guide to anomalous diffusion:a fractional dynamics approach[END_REF]). In this context the recovery of the source term f can be seen as the recovery of a time evolving source of diffusion. For instance, in the case α = 1 and d = 3, our problem can be associated with the recovery of a source moving in the subset ω of the plan R 2 from a single measurement of the heat flux at the boundary. Such a problem can be associated with the determination of different properties of materials such as metal (see e.g., [START_REF] Kim | An analytical solution to heat conduction with a moving heat source[END_REF] for the heat equation). We refer to [START_REF] Hu | Inverse moving source problems in electrodynamics[END_REF] for applications of recovery of source terms of the form (1.8) to the recovery of moving sources in the electrodynamics. For non-integer value of the power α, in the spirit of [START_REF] Nakagawa | Overview to mathematical analysis for fractional diffusion equations-new mathematical aspects motivated by industrial collaboration[END_REF] (see also [START_REF] Jin | A tutorial on inverse problems for anomalous diffusion processes[END_REF]), our inverse problem can be seen as the recovery of a moving source of diffusion of a contaminant under the ground. As unknown sources, we assume the form of (1.8), which can be interpreted for example that an uknown source f (x , t) depends only on the depth variable and t in the case of d = 2, which corresponds to a layer structure, and on the planar locations (x 1 , x 2 ) and t but not on the depth in the case of d = 3, which may be a good approximation if Ω is a very thin domain in the direction of x 3 .

1.4. Known results. Inverse source problems have received a lot of attention these last decades among the mathematical community (see [START_REF] Isakov | Inverse source problems[END_REF] for an overview). For diffusion equations corresponding to the case α = 1 with time independent source terms, several authors investigated the conditional stability (e.g. [START_REF] Choulli | Some stability estimates in determining sources and coefficients[END_REF][START_REF] Yamamoto | Conditional stability in the determination of force terms of heat equations in a rectangle[END_REF][START_REF] Yamamoto | Conditional stability in the determination of densities of heat sources in a bounded domain, Estimation and Control of Distributed Parameter Systems[END_REF]). Following the Bukhgeim-Klibanov approach introduced in [3], [START_REF] Yu | Lipschitz stability in inverse parabolic problems by the Carleman estimate[END_REF] established Lipschitz stable recovery of the source from one Neumann boundary measurement. In [START_REF] Choulli | Some stability estimates in determining sources and coefficients[END_REF], the authors derived a stability estimate for this problem from a single Neumann observation of the solution on an arbitrary portion of the boundary. This last result has been extended to an infinite cylindrical domain by [START_REF] Kian | Logarithmic stability inequality in an inverse source problem for the heat equation on a waveguide[END_REF]. For fractional diffusion equations corresponding to the case α = 1, in spite of the phyiscal backgrounds related to various anomalous diffusion phenomena stated above, to our best knowledge, there is no result in the mathematical literature dealing with the recovery of source terms, depending on both time and space variables, of the form (1.8). We mention also the work of [START_REF] Jiang | Weak unique continuation property and a related inverse source problem for time-fractional diffusion-advection equations[END_REF][START_REF] Jin | An inverse problem for a one-dimensional time-fractional diffusion problem[END_REF][START_REF] Jin | A tutorial on inverse problems for anomalous diffusion processes[END_REF][START_REF] Kian | Global uniqueness in an inverse problem for timefractional diffusion equations[END_REF][START_REF] Kian | On time-fractional diffusion equations with space-dependent variable order[END_REF][START_REF] Li | Initial-boundary value problem for distributed order time-fractional diffusion equations[END_REF][START_REF] Rundell | Recovering an unknown source in a fractional diffusion problem[END_REF] where some inverse coefficient problems and some related results have been considered.

The above mentioned results are all concerned with the determination of time-independent source terms f (x). Several authors considered also the recovery of time-dependent source terms. In [START_REF] Fujishiro | Determination of time dependent factors of coefficients in fractional diffusion equations[END_REF][START_REF] Sakamoto | Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems[END_REF] the authors proved the stable recovery of a source term f (t) depending only on the time variable from measurements of solutions at one spatial point over time interval. As long as the classical partial differential equations with natural number α are concerned, some papers have also been devoted to the unique existence and the stability in finding source terms of the form (1.8) (see e.g. [START_REF] Ya | Determination of coefficients of higher terms in a parabolic equation[END_REF][START_REF] Hu | Uniqueness and stability for the recovery of a time-dependent source and initial conditions in elastodynamics[END_REF][START_REF] Isakov | Inverse source problems[END_REF]). In particular, in [START_REF] Isakov | Inverse source problems[END_REF]Section 6.3] the author proved the reconstruction and the unique recovery of source terms of the form (1.8) appearing in a parabolic equation on the half space.

Our main results stated below, seem the first achievements for the inverse source problem of determining f (t, x ) in (1.8) for fractional partial differential equations.

1.5. Preliminary properties. In the present paper, following [START_REF] Kian | On existence and uniqueness of solutions for semilinear fractional wave equations[END_REF][START_REF] Sakamoto | Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems[END_REF], we consider solutions of problem (1.2)-(1.6) (resp. (1.5)-(1.7)) in the following weak sense.

Definition 1.1. Let F ∈ L 2 (Q) (resp. F ∈ L 2 ((0, T ) × Ω)).
We say that problem (1.2)-(1.4) (resp.

(1.5)-(1.7) with q = 0) admits a weak solution u if there exists v ∈ L 2 loc (R + ; L 2 (Ω)) such that:

1) v |Q = u and inf{ε > 0 : e -εt v ∈ L 1 (R + ; L 2 (Ω))} = 0,
2) for all p > 0 the Laplace transform V (p) = +∞ 0 e -pt v(t, .)dt with respect to t of v, satisfies (1.3)-(1.4) (resp. (1.6)-(1.7)) and solves

(A + p α )V (p) = F (p) in Ω, where F (p) = L[F (t, .)1 (0,T ) (t)](p) =
T 0 e -pt F (t, .)dt and 1 (0,T ) is the characteristic function of (0, T ).

By the results in [START_REF] Lions | Non-homogeneous Boundary Value Problems and Applications[END_REF][START_REF] Sakamoto | Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems[END_REF], we can prove that for

F ∈ L 2 (Q) problem (1.2)-(1.4) (resp. (1.5)-(1.7)) admits a unique solution u ∈ L 2 (0, T ; H 1 (Ω)) satisfying ∂ α t u, Au ∈ L 2 (Q).
For sake of completeness we recall this result in the Appendix.

In (1.5), we assume that a ij ∈ C 3 (Ω), i, j ∈ {1, . . . , d}, satisfy the following condition

a ij (x , -x d ) = a ij (x , x d ), (x , x d ) ∈ Ω, i, j ∈ {1, . . . , d}.
(1.9)

Then we consider the problem

                   ∂ α t u + Au = F (t, x), (t, x , x d ) ∈ (0, T ) × Ω, u(t, x) = 0, (t, x) ∈ (0, T ) × [∂ω × (0, ) ∪ ω × { }], -∂ x d u(t, x , 0) = 0 (t, x ) ∈ (0, T ) × ω, ∂ k t u(0, x) = 0, x ∈ Ω, k = 0, m α .
(1.10)

Here we consider weak solutions in the sense of Definition 1.1 and, following [START_REF] Kian | On existence and uniqueness of solutions for semilinear fractional wave equations[END_REF][START_REF] Sakamoto | Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems[END_REF], we can prove that there exists an operator valued function J(t) ∈ B(L 2 ( Ω)) such that the solution of (1.10) takes the form

u(t, •) = t 0 J(t -s)F (s, •)ds.
Using this definition, for q ∈ L ∞ ((0, T ) × ω) we can define the solution of

                   ∂ α t u + Au + q(t, x )u = F (t, x), (t, x , x n ) ∈ (0, T ) × Ω, u(t, x) = 0, (t, x) ∈ (0, T ) × [∂ω × (0, ) ∪ ω × { }], -∂ x d u(t, x , 0) = 0 (t, x ) ∈ (0, T ) × ω, ∂ k t u(0, x) = 0, x ∈ Ω, k = 0, . . . , m α , (1.11)
in the mild sense which is defined by the integral equation

u(t, •) = - t 0 J(t -s)qu(s, •)ds + t 0 J(t -s)F (s, •)ds.
1.6. Main results. From now on, we assume that F takes the form (1.8). For our first result we need an assumption on ω and A that guarantees the elliptic regularity of the operator A. Indeed, due to the fact that the domain Ω is only Lipschitz, some extra assumptions will be required for guaranteeing the elliptic regularity of A with the boundary conditions (1.3)- (1.4). For this purpose, for m, n = 0, 1, we introduce the condition:

(Hmn) For all v ∈ H 1 (Ω) satisfying Av ∈ L 2 (Ω) and (1.3)-(1.4) with these value of m, n, we have v ∈ H 2 (Ω) and there exists C > 0 depending only on A and Ω such that

v H 2 (Ω)) C( Av L 2 (Ω)) + v L 2 (Ω)) ).
Note that conditions (H00) and (H11) will be fulfilled if, for instance, we assume that ω is convex.

Indeed, in that case Ω will also be convex and, in virtue of [9, Theorem 3.2.1.2] and [9, Theorem 3.2.1.3],

(H00) and (H11) will be fulfilled. In the same way, assuming that

a id = 0, ∂ xj a dd = 0, ∂ x d a ij = 0, i, j ∈ {1, . . . , d -1}, (1.12)
we deduce from a separation of variable argument similar to [START_REF] Gaitan | A stability result for a time-dependent potential in a cylindrical domain[END_REF]Lemma 2.4] that, for all m, n = 0, 1, (Hmn) is fulfilled.

Using the conditions (Hmn), we obtain the following.

Theorem 1.2. Let (H00), (H10) be fulfilled and assume that R, ∂ x d R ∈ L ∞ (Q) and there exists a constant c > 0 such that

|R(t, x , )| c, (t, x ) ∈ (0, T ) × ω. (1.13)
Assume also that the condition

∂ x d a ij = 0, i, j = 1, . . . , d -1, (1.14) is fulfilled. Then, for f ∈ L 2 ((0, T ) × ω), the solution u of (1.2)-(1.4) with m = 1 and n = 0 satisfies u ∈ H 3 (-, ; L 2 ((0, T ) × ω)), ∂ α t u, Au ∈ H 1 (-, ; L 2 ((0, T ) × ω))
. Therefore, we can define

h := (t, x ) → [∂ α t u + (A + a dd ∂ 2 x d )u](t, x , ) R(t, x , ) ∈ L 2 ((0, T ) × ω). (1.15)
Moreover, for every f ∈ L 2 ((0, T ) × ω), we can define an operator valued function

K(t, s) ∈ B(L 2 (ω)),
t, s ∈ (0, T ), such that f solves the integral equation

f (t, •) = h(t, •) + t 0 K(t, s)f (s, •)ds, t ∈ (0, T ), (1.16)
which is well-posed. Finally, for every [START_REF] Jin | An inverse problem for a one-dimensional time-fractional diffusion problem[END_REF]. In the same way, assuming that (H01) and (H11) are fulfilled, the same results hold true for the problem (1.2)-(1.4) with m = 1 and n = 1.

(h, f ) ∈ L 2 ((0, T ) × ω) × L 2 ((0, T ) × ω) satisfying (1.16) the solution u of (1.2)-(1.4) satisfies (1.
For problem (1.11), we consider first the following condition:

( H) For all v ∈ H max(1,s) ( Ω) satisfying Av ∈ H s ( Ω), s ∈ [0, 2], and (1.6)-(1.7), we have v ∈ H 2+s (Ω)
and there exists C > 0 depending only on A, s and Ω such that

v H 2+s (Ω)) C( Av H s (Ω)) + v H s (Ω)) ).
Assuming that ω is of class C 4 and using a separation of variable argument similar to [8, Lemma 2.4],

one can check that (1.12) implies ( H).

Using ( H) we obtain the following well-posedness result.

Proposition 1.3. Assume that ( H) is fulfilled. Let γ ∈ (0, 1) be such that for α ∈ (0, 1], γ ∈ (1/2, 1
)

and for α ∈ (1, 2), γ ∈ (1/2, 1/α). Fix p ∈ (1, +∞) such that 1 p < min(1 -αγ, α(1 -γ)) and let F ∈ W 1,p (0, T ; L 2 ( Ω)) ∩ C([0, T ]; H 2γ ( Ω)) satisfy F |t=0 = 0. Let q ∈ C 1 ([0, T ]; L ∞ (ω)) ∩ C([0, T ]; W 2,∞ (ω)). Then problem (1.11) admits a unique weak solution u ∈ C([0, T ]; H 2(1+γ) ( Ω)) ∩ C 1 ([0, T ]; H γ ( Ω)).
Applying this well-posedness result, we can state our second main result as follows.

Theorem 1.4. Assume that ( H), (H00), (1.9) and (1.14) are fulfilled, α ∈ (0, 1], ω is of C 4 and

γ ∈ (3/4, 1). Fix p ∈ (1, +∞) such that 1 p < min(1 -αγ, α(1 -γ)). Let q ∈ C 1 ([0, T ]; L ∞ (ω)) ∩ C([0, T ]; W 2,∞ (ω)), f ∈ W 1,p (0, T ; L ∞ (ω))∩C([0, T ]; W 2,∞ (ω)) and R ∈ W 1,∞ (0, T ; L 2 ( Ω))∩C([0, T ]; H 2γ ( Ω)) satisfy R, ∂ x d R ∈ L ∞ ((0, T ) × Ω) and (1.13). Assume also that f (0, x) = 0, x ∈ Ω.
Then, for δ ∈ 0, 2γ - 3 4 , there exists a constant C depending on q, R, Ω, T , α, A, δ, such that

f L ∞ (0,T ;L 2 (ω)) C( ∂ x d u(•, •, ) L ∞ (0,T ;H 3 2 (ω)) + ∂ x d u(•, •, ) W 1,∞ (0,T ;H δ (ω)) ).
(1.17)

Applying this result, we can also prove the stable recovery of the coefficient q appearing in the problem

                         ∂ α t v + Av + q(t, x )v = 0, (t, x) ∈ (0, T ) × Ω, v(t, x , ) = h 0 (t, x ) (t, x ) ∈ (0, T ) × ω, v(t, x) = h 1 (t, x), (t, x) ∈ (0, T ) × ∂ω × (-, ), ∂ x d v(t, x , -) = 0 (t, x ) ∈ (0, T ) × ω, v(0, x) = w 0 , x ∈ Ω (1.18) with α ∈ (0, 1], h k , k = 0, 1, w 0 such that there exists H ∈ C 1 ([0, T ]; H 2+γ ( Ω)) ∩ W 2,p (0, T ; H 2 ( Ω)) satisfying                          H(t, x , ) = h 0 (t, x ) (t, x ) ∈ (0, T ) × ω, H(t, x) = h 1 (t, x), (t, x) ∈ (0, T ) × ∂ω × (0, ), ∂ x d H(t, x , 0) = H(t, x , 0) = 0 (t, x ) ∈ (0, T ) × ω, H(0, x) = w 0 (x), x ∈ Ω, ∂ α t H(0, x) = -AH(0, x) -q(0, x )H(0, x), x = (x , x d ) ∈ Ω. (1.19)
The result for the determination of q can be stated as follows.

Corollary 1.5. Let the condition of Theorem 1.4 be fulfilled with f = 0, d 3 and let h k , k = 0, 1, w 0 be given by (1.19) and the condition

|h 0 (t, x )| c, (t, x ) ∈ (0, T ) × ω (1.20) be fulfilled. Fix q j ∈ C 1 ([0, T ]; L ∞ (ω)) ∩ C([0, T ]; W 2,∞ (ω)), j = 1, 2, such that q j W 1,∞ (0,T ;L ∞ (ω)) + q j L ∞ (0,T ;W 2,∞ (ω)) M, q 1 (0, x ) = q 2 (0, x ) = q(0, x ), x ∈ ω
and consider the solution v j of (1.18) with q = q j . Then, for δ ∈ 0, 2γ -3 4 , we have

q 1 -q 2 L ∞ (0,T ;L 2 (ω)) C ∂ x d (v 1 -v 2 ) |x d = L ∞ (0,T ;H 3 2 (ω)) + ∂ x d (v 1 -v 2 ) |x d = W 1,∞ (0,T ;H δ (ω)) , (1.21)
where the constant C depends on α, c, T , Ω, a, A, α, M , h 0 , h 1 , w 0 , δ. with Ω corresponding to the half space and see also [START_REF] Ya | Determination of coefficients of higher terms in a parabolic equation[END_REF]. In contrast to [START_REF] Isakov | Inverse source problems[END_REF], we state our result on a bounded cylindrical domain and we restrict our analysis to solutions lying in Sobolev spaces while [13, Section 6.3] is stated with Hölder continuous functions. Moreover our approach admits a natural extension to fractional diffusion equations (α = 1).

Let us remark that Theorem 1.2 gives a reconstruction algorithm for the recovery of the source term f under consideration. It is actually stated as a well-posedness result for the pair of functions (u, f ) appearing in (1.2) and (1.8). In contrast to Theorem 1.2, Theorem 1.4 provides only a stability estimate. However, Theorem 1.4 can be applied to more general boundary conditions and it can also be applied to the stable recovery of a coefficient depending on both time and space variables (see Corollary 1.5).

Applying Theorem 1.4, we prove in Corollary 1.5 the stable recovery of the coefficient of order zero q provided ∂ x d q = 0. It seems that this result is the first result of stable recovery of a coefficient depending on both time and space variables for a fractional diffusion equation. In [START_REF] Kian | On existence and uniqueness of solutions for semilinear fractional wave equations[END_REF]Theorem 3.6], the authors derived a similar result for the heat equation (α = 1) stated with stronger regularity conditions and measurements on both γ × (0, ) and ω × { 1 }, with γ ⊂ ∂ω and 1 ∈ (0, ]. Even, for α = 1, our result improves the one of [START_REF] Kian | On existence and uniqueness of solutions for semilinear fractional wave equations[END_REF]Theorem 3.6] in terms of regularity conditions and restriction of the data.

Theorem 1.4 is stated only for α ∈ (0, 1] but it can be extended to α ∈ (0, 2) without any difficulty. Indeed, one can easily extend our argumentation to the case α ∈ (1, 2). In Theorem 1.4

we have restricted our analysis to α ∈ (0, 1] in order to simplify the statement of this theorem and its proof.

1.8. Outline. This paper is organized as follows. In Section 2 we prove Theorem 1.2. Section 3 is devoted to the proof of Theorem 1.4 and in Section 4 we consider the application of Theorem 1.4 stated in Corollary 1.5. Finally, in the Appendix we recall and prove a result related to regularity of solutions of (1.2).

Proof of Theorem 1.2

We start with considering the first part of Theorem 1.2. For this purpose, we assume that (H00), (H10), (1.13)-(1.14) are fulfilled and we will show (1.16). We denote by A (resp. Ã) the operator A acting on L 2 (Ω) with domain

D(A) := {g ∈ L 2 (Ω) : Ag ∈ L 2 (Ω), u |∂Ω = 0} (resp. D( Ã) := {g ∈ L 2 (Ω) : Ag ∈ L 2 (Ω), u |∂ω×(0, ) = 0, ∂ x d u |x d =± = 0}
Thanks to (1.1) we know that A (resp. Ã) are selfadjoint operators with a spectrum consisting of a non-decreasing sequence of non-negative eigenvalues (λ n ) n 1 (resp. ( λn ) n 1 ). Moreover, conditions (H00) and (H10) imply that D(A) and D( Ã) embedded continuously into H 2 (Ω). Let us also introduce an orthonormal basis in the Hilbert space L 2 (Ω) of eigenfunctions (ϕ n ) n 1 (resp. ( φn ) n 1 ) of A (resp.

Ã) associated to the non-decreasing sequence of eigenvalues (λ n ) n 1 (resp. ( λn ) n 1 ). We consider also the operator valued function S(t) (resp. S(t)) defined by

S(t)h = ∞ n=1 t α-1 E α,α (-λ n t α ) h, ϕ n L 2 (Ω) ϕ n , h ∈ L 2 (Ω), S(t)h = ∞ n=1 t α-1 E α,α (-λn t α ) h, φn L 2 (Ω) φn , h ∈ L 2 (Ω),
where E α,α corresponds to the Mittag-Leffler function given by

E α,α (z) = ∞ n=0 z n Γ(α(n + 1)) , z ∈ C.
Following [START_REF] Kian | On existence and uniqueness of solutions for semilinear fractional wave equations[END_REF][START_REF] Sakamoto | Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems[END_REF] (see also Lemma 3.1), one can check that problem (1.2) admits a unique weak solution u ∈ L 2 (0, T ; D( Ã)) taking the form

u(t, •) = t 0 S(t -s)F (s, •)ds, t ∈ (0, T ).
Recall that the function v given by Definition 1.1 takes the form

v(t, •) = min(t,T ) 0 S(t -s)F (s, •)ds, t ∈ (0, +∞).
Moreover, using the fact that D( Ã) embedded continuously into H 2 (Ω), we deduce by interpolation that for s 1 ∈ 3 4 , 1 , D( Ãs1 ) embedded continuously into H 2s1 (Ω). Therefore, applying formula (1.148) of [START_REF] Podlubny | Fractional differential equations[END_REF]Theorem 1.6], one can check that inf{ε > 0 :

e -εt v ∈ L 1 (R + ; H 2s1 (Ω))} = 0. Fixing w = ∂ x d u
we deduce that, for all p > 0, the Laplace transform in time W (p) of the extension of w to R + × Ω given by w = ∂ x d v, with v defined in Definition 1.1, is lying in H 2s1-1 (Ω) and it satisfies

   (A + p α )W (p) = -(∂ x d A)V (p) + ∂ x d F (p), in Ω, W (p) = 0, on ∂Ω.
Note that here we use the fact that W (p) ∈ H 2s1-1 (Ω), with 2s 1 -1 > 1 2 , for defining its trace on ∂Ω. Moreover, applying formula (1.148) of [START_REF] Podlubny | Fractional differential equations[END_REF]Theorem 1.6] and the fact that D(A 1/2 ) = H 1 0 (Ω), we can extend S(t), t > 0, to a bounded operator from H -1 (Ω) to L 2 (Ω) satisfying

S(t) B(H -1 (Ω),L 2 (Ω)) Ct α/2-1 .
In addition, using the fact that v ∈ L 2 loc (R

+ ; D(A)) ⊂ L 2 loc (R + ; H 2 (Ω)), we deduce that (∂ x d A)v ∈ L 2 loc (R + ; L 2 (Ω))
. Thus, extending F by zero to (0, +∞) × Ω and using the fact that

inf{ε > 0 : e -εt (∂ x d A)v ∈ L 1 (R + ; H -1 (Ω))} = 0,
we deduce that the function

w 1 (t, •) := t 0 S(t -s)[(∂ x d A)v(s, •)) + ∂ x d F (s, •)]ds = t 0 S(t -s 1 )(∂ x d A) s1 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 ds 1 + t 0 S(t -s)f (s 2 , •)∂ x d R(s 2 , •)ds.
is well defined and the Laplace transform in time of w 1 coincide with the one of w. This proves that

∂ x d u = - t 0 S(t -s 1 )(∂ x d A) s1 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 ds 1 + t 0 S(t -s)f (s 2 , •)∂ x d R(s 2 , •)ds.
(2.1)

With the additional assumptions (H10), (H11) and the fact that ω is C 2 one can extend these arguments to problem (1.2)-(1.4), with m = n = 1, by using the fact that, for any

s 1 ∈ (0, 1/2), C ∞ 0 (Ω) is dense in H s 1 (Ω).
In view of (1.14), for w = ∂ x d u, we have

-(∂ x d A)u = ∂ x d (∂ x d a dd ∂ x d u) + d-1 j=1 [∂ xj (∂ x d a jd ∂ x d u) + ∂ x d (∂ x d a jd ∂ xj u)] = B 1 w + B 2 u, (2.2)
where

B 1 w := ∂ x d a dd ∂ x d w + 2 d-1 j=1 ∂ x d a jd ∂ xj w +   d j=1 ∂ xj ∂ x d a jd   w, (2.3)
B 2 u = d j=1 (∂ 2 x d a jd )∂ xj u. (2.4)
Then, from (2.1) and the above arguments, we deduce that w solves the integral equation

w(t, •) = t 0 S(t -s)B 1 w(s, •)ds + t 0 S(t -s 1 )B 2 s1 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 ds 1 + t 0 S(t -s)f (s, •)∂ x d R(s, •)ds.
(2.5)

Now let us consider the following.

Lemma 2.1. The integral equation (2.5) admits a unique solution w ∈ L 2 (0, T ; H 2 (Ω)) satisfying

w(t, •) H 1 (Ω) C t 0 (t -s) α/2-1 f (s, •) L 2 (ω) , t ∈ (0, T ), (2.6)
with C > 0 depending on R, A, T . Moreover, we have

∂ α t w ∈ L 2 (Q).
Proof. We introduce the maps K 1 , M : L 2 (0, T ;

H 1 (Ω)) -→ L 2 (0, T ; H 1 (Ω)), K 2 : L 2 ((0, T ) × ω)) -→ L 2 (0, T ; H 1 (Ω)) defined by K 1 v(t, •) := t 0 S(t -s)B 1 v(s, •)ds, K 2 f (t, •) := t 0 S(t -s 1 )B 2 s1 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 ds 1 + t 0 S(t -s)f (s, •)∂ x d R(s, •)ds
and Mv = K 1 v + K 2 f . Applying formula (1.148) of [29, Theorem 1.6] and the fact that both D(A

)

and D( Ã

2 ) embedded continuously into H 1 (Ω), we find

K 1 v(t, •) H 1 (Ω) C K 1 v(t, •) D(A 1 2 ) C t 0 (t -s) α 2 -1 Γ(α/2) v(s, •) H 1 (Ω) ds, t ∈ (0, T ) (2.7) K 2 f (t, •) H 1 (Ω) C t 0 (t -s) α 2 -1 f (s, •) L 2 (ω) ds + t 0 s1 0 (t -s 1 ) α 2 -1 (s 1 -s 2 ) α 2 -1 f (s 2 , •) L 2 (ω) ds 2 ds 1 C t 0 (t -s) α 2 -1 f (s, •) L 2 (ω) ds, t ∈ (0, T ).
(2.8)

Following the proof of [7, Proposition 1], for n ∈ N, we find by iteration

K n 1 v(t, •) H 1 (Ω) C n t 0 (t -s) nα 2 -1 Γ(nα/2) v(s, •) H 1 (Ω) ds, t ∈ (0, T ) (2.9) K n-1 1 K 2 f (t, •) H 1 (Ω) C n t 0 (t -s) nα 2 -1 Γ(nα/2) f (s, •) L 2 (ω) ds, t ∈ (0, T ).
(2.10)

It follows that for n ∈ N sufficiently large we have

M n v 1 (t, •) -M n v 2 (t, •) H 1 (Ω) C n t 0 (t -s) nα 2 -1 Γ(nα/2) v 1 (s, •) -v 2 (s, •) H 1 (Ω) ds, t ∈ (0, T )
and an application of the Young inequality for convolution product implies that

M n v 1 -M n v 2 L 2 (0,T ;H 1 (Ω)) C n T 0 t nα 2 -1 Γ(nα/2) dt v 1 -v 2 L 2 (0,T ;H 1 (Ω)) C n T nα 2 nα 2 Γ(nα/2) v 1 -v 2 L 2 (0,T ;H 1 (Ω)) C n T nα 2 +1 Γ(nα/2 + 1) v 1 -v 2 L 2 (0,T ;H 1 (Ω)) .
Then, using the fact that

lim n→+∞ C n T nα 2 +1
Γ(nα/2 + 1) = 0, we deduce that there exists n 0 ∈ N such that M n0 is a contraction. Moreover, conditions (2.9)-(2.10) imply

M n v L 2 (0,T ;H 1 (Ω)) C n T nα 2 +1 Γ(nα/2 + 1) v L 2 (0,T ;H 1 (Ω)) + n k=1 C k T kα 2 +1 Γ(kα/2 + 1) f L 2 ((0,T )×ω) C n T nα 2 +1 Γ(nα/2 + 1) v L 2 (0,T ;H 1 (Ω)) + ∞ k=1 C k T kα 2 +1
Γ(kα/2 + 1) f L 2 ((0,T )×ω) .

Therefore, by eventually increasing the size of n 0 , we deduce that M n0 admits a unique fixed point w ∈ L 2 (0, T ; H 1 (Ω)) which by uniqueness of this fixed point is also a fixed point of M. Moreover, in view of (2.9)-(2.10), for n 1 ∈ N satisfying n 1 4 α and for a.e t ∈ (0, T ), we have

w(t, •) H 1 (Ω)) = M n1 w(t, •) H 1 (Ω)) C t 0 (t -s) n 1 α 2 -1 w(s, •) H 1 (Ω) ds + t 0 (t -s) α 2 -1 f (s, •) L 2 (ω) ds C T n 1 α 2 -1 t 0 w(s, •) H 1 (Ω) ds + t 0 (t -s) α 2 -1 f (s, •) L 2 (ω) ds .
Therefore, applying Gronwall inequality for function lying in L 2 (0, T ) (see e.g [1, Lemma 6.3]), we find

w(t, •) H 1 (Ω)) C t 0 (t -s) α 2 -1 f (s, •) L 2 (ω) ds e Ct , t ∈ (0, T ),
which clearly implies (2.6). Finally, using the fact that

B 1 w, s 1 → B 2 s1 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 , ∂ x d F ∈ L 2 (Q),
we deduce from Lemma 3.1 and assumption (H00), that M takes values in L 2 (0, T ; H 2 (Ω)) and therefore w ∈ L 2 (0, T ; H 2 (Ω)). In the same way, we prove that

∂ α t w ∈ L 2 (Q).
According to Lemma 2.1, we have u ∈ L 2 (0, T ;

H 1 (Ω)), Au, ∂ α t u ∈ L 2 (Q) and w = ∂ x d u ∈ L 2 (0, T ; H 2 (Ω)), ∂ α t ∂ x d u, A∂ x d u ∈ L 2 (Q).
Combining this with (2.2), we deduce that

x d → u(•, •, x d ) ∈ H 3 ((-, ); L 2 ((0, T ) × ω)) ∩ H 1 (-, ; L 2 (0, T ; H 2 (ω))), x d → ∂ α t u(•, •, x d ) ∈ H 1 ((-, ); L 2 ((0, T ) × ω)).
Then, (1.2) implies

a dd ∂ x d w(t, x , ) = a dd ∂ 2 x d u(t, x , ) = [∂ α t u + (A + a dd ∂ 2 x d )u](t, x , ) -R(t, x , )f (t, x ) = R(t, x , )[h(t, x ) -f (t, x )], (t, x ) ∈ (0, T ) × ω, (2.11)
with h given by (1.15). Fixing r ∈ (3/4, 1) and τ 1 :

L 2 (0, T ; H 2r (Ω)) y -→ a dd ∂x d y(t,x , ) R(t,x , ) ∈ L 2 ((0, T )× ω), we obtain f (t, •) = h(t, •) -τ 1 w(t, •), t ∈ (0, T ).
Moreover, applying (2.2), (2.5) and using the fact that u, w ∈ L 2 (0, T ; H 1 (Ω)), we get

(∂ x d A) s1 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 ds 1 = (∂ x d A)u ∈ L 2 (Q), (2.12)
which combined with (H00) implies

f (t, •) =h(t, •) -τ 1 t 0 S(t -s 1 )(∂ x d A) s1 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 ds 1 -τ 1 t 0 S(t -s)f (s, •)∂ x d R(s, •)ds , t ∈ (0, T ).
(2.13)

In view of (H00) by interpolation D(A r ) embedded continuously into H 2r (Ω). Moreover, in view of formula (1.148) of [29, Theorem 1.6], we have

S(t -s 1 )(∂ x d A) s1 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 H 2r (Ω) C S(t -s 1 )(∂ x d A) s1 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 D(A r ) C(t -s 1 ) α(1-r)-1 (∂ x d A) s1 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 L 2 (Ω)
, s 1 ∈ (0, t).

On the other hand, applying Young inequality for convolution product, we obtain

T 0 t 0 S(t -s 1 )(∂ x d A) s1 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 H 2r (Ω) ds 1 dt C T 0 t 0 (t -s 1 ) α(1-r)-1 (∂ x d A) t 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 L 2 (Ω) ds 1 dt CT 1-r (∂ x d A) t 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 L 1 (0,T ;L 2 (Ω)) C (∂ x d A) t 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 L 2 (Q) < ∞.
Thus, by Fubini theorem for a.e. t ∈ (0, T ), we have

s 1 -→ S(t -s 1 )(∂ x d A) s1 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 ∈ L 1 (0, t; H 2r (Ω)).
In the same way, for a.e. t ∈ (0, T ), we get

s -→ S(t -s)f (s, •)∂ x d R(s, •) ∈ L 1 (0, t; H 2r (Ω)).
Therefore, from (2.13), for a.e. t ∈ (0, T ), we find

f (t, •) =h(t, •) - t 0 τ 1 S(t -s 1 )(∂ x d A) s1 0 S(s 1 -s 2 )f (s 2 , •)R(s 2 , •)ds 2 ds 1 - t 0 τ 1 S(t -s)f (s, •)∂ x d R(s, •)ds.
(2.14)

This proves (1.16), let us prove that this problem is well-posed. We fix the maps G, H :

L 2 ((0, T )×ω) -→ L 2 ((0, T ) × ω), with
Hg(t, •) := - t 0 τ 1 S(t -s 1 )(∂ x d A) s1 0 S(s 1 -s 2 )g(s 2 , •)R(s 2 , •)ds 2 ds 1 - t 0 τ 1 S(t -s)g(s, •)∂ x d R(s, •)ds, g ∈ L 2 ((0, T ) × ω)
and

Gg(t, •) := h(t, •) + Hg(t, •), g ∈ L 2 ((0, T ) × ω).

Note that

Hg(t, •) L 2 (ω) t 0 τ 1 S(t -s 1 )(∂ x d A) s1 0 S(s 1 -s 2 )g(s 2 , •)R(s 2 , •)ds 2 L 2 (ω) ds 1 + t 0 τ 1 S(t -s)g(s, •)∂ x d R(s, •) L 2 (ω) ds C t 0 S(t -s 1 )(∂ x d A) s1 0 S(s 1 -s 2 )g(s 2 , •)R(s 2 , •)ds 2 H 2r (Ω) ds 1 + C t 0 S(t -s)g(s, •)∂ x d R(s, •) H 2r (Ω) ds C(t α(1-r)-1 1 (0,+∞) ) * (∂ x d A) 1 (0,T ) (t) t 0 S(t -s)g(s, •)R(s, •)ds L 2 (Ω) + C ∂ x d R L ∞ (Q) (t α(1-r)-1 1 (0,+∞) ) * ( g(•, ) L 2 (ω) 1 (0,T ) ), t ∈ (0, T ). (2.15)
On the other hand, in view of (2.2), fixing w 1 ∈ L 2 (0, T ; H 1 (Ω)) the solution of (2.5) with f = g we obtain

(∂ x d A) t 0 S(t -s)g(s, •)R(s, •)ds = B 1 w 1 (t, •) + B 2 t 0 S(t -s)gR(s, •)ds,
where B 1 and B 2 are defined in formula (2.3)-(2.4). Thus, applying Lemma 2.1 and formula (1.148) of [29, Theorem 1.6], we obtain

(∂ x d A) t 0 S(t -s)g(s, •)R(s, •)ds L 2 (Ω) C t 0 (t -s) α 2 -1 g(s, •) L 2 (ω) ds, t ∈ (0, T ).
Combining this with (2.15), we get

Hg(t, •) L 2 (ω) C t 0 (t -s) α(1-r)-1 Γ(α(1 -r)) g(s, •) L 2 (ω) ds, t ∈ (0, T ).
By iteration, for all n ∈ N, we deduce that

H n g(t, •) L 2 (ω) C n t 0 (t -s) nα(1-r)-1 Γ(nα(1 -r)) g(s, •) L 2 (ω) ds, t ∈ (0, T )
and in a similar way to Lemma 2.1, we deduce that there exists n 2 ∈ N such that G n2 is a contraction and G admits a unique fixed point f ∈ L 2 ((0, T ) × ω) satisfying

f (t, •) L 2 (ω) = Gf (t, •) L 2 (ω) h(t, •) L 2 (ω) + C t 0 (t -s) α(1-r)-1 f (s, •) L 2 (ω) ds, t ∈ (0, T ).
Therefore, in view of [7, Lemma 3] (see also [33, Theorem A]), we have

f (t, •) L 2 (ω) C h(t, •) L 2 (ω) + t 0 (t -s) α(1-r)-1 h(s, •) L 2 (ω) ds , t ∈ (0, T ),
and an application of the Young inequality yields

f L 2 ((0,T )×ω) C h L 2 ((0,T )×ω) .
This proves the well-posedness of (1.14) and the reconstruction of f from the data h. Now let us consider the proof of the last part of Theorem 1.2. For this purpose, we fix

(h 1 , f ) ∈ L 2 ((0, T )×ω)×L 2 ((0, T )×ω)
satisfying (1.16) with h = h 1 and we consider u ∈ L 2 (0, T ; H 1 (Ω)) solving (1.2). Following the above argumentation, we can define h ∈ L 2 ((0, T ) × ω) given by (1.15) and f solves (1.16). This implies that

h 1 (t, •) = f - t 0 K(t, s)f (s, •)ds = h(t, •), t ∈ (0, T ).
Therefore we have h = h 1 and the proof of Theorem 1.2 for (1.2)-(1.4), with m = 1 and n = 0, is completed. Using similar arguments, one can check that this result is still true for the problem (1.2)-(1.4), with m = 1 and n = 1.

Proof of Theorem 1.4

This section is devoted to the proof of Theorem 1.4. In contrast to the preceding section, for u solving (1.11), ∂ x d u will not be a solution of an initial boundary value problem with homogeneous boundary condition. However, with suitable regularity conditions on F we can consider the trace of u at {x d = }. We will start by proving Proposition 1.3

Proof of Proposition 1.3. We consider first the case q = 0. Let A be the operator A acting on

L 2 ( Ω) with domain D(A) := {g ∈ L 2 ( Ω) : Ag ∈ L 2 ( Ω), u |∂ω×(0, ) = 0, u |x d = = 0, ∂ x d u |x d =0 = 0}.
The spectrum of A consists of a non-decreasing sequence of strictly positive eigenvalues (λ n ) n 1 . Let us also introduce an orthonormal basis in the Hilbert space L 2 ( Ω) of eigenfunctions (ϕ n ) n 1 of A associated with the non-decreasing sequence of eigenvalues (λ n ) n 1 . Then, for n ∈ N, the solution of (1.11) is given by

u n (t) := u(t, •), ϕ n L 2 ( Ω) = t 0 (t -s) α-1 E α,α (-λ n (t -s) α )F n (s)ds,
where F n (t) := F (t, •), ϕ n L 2 ( Ω) . Since F n ∈ W 1,p (0, T ), F n (0) = 0, applying [30, Lemma 3.2] and integrating by parts we find

u n (t) = -λ -1 n F n (t) + λ -1 n t 0 E α,1 (-λ n (t -s) α )F n (s)ds.
Thus, we have u = v + w where

v(t, •) = A -1 F (t, •), w n (t) := w(t, •), ϕ n L 2 ( Ω) = λ -1 n t 0 E α,1 (-λ n (t -s) α )F n (s)ds.
Note that, for all t ∈ [0, T ], v(t, •) solves the boundary value problem

                   Av(t, •) = F (t, •), in Ω, v(t, x) = 0, x ∈ ∂ω × (0, ), v(t, x , ) = 0, x ∈ ω, ∂ x d v(t, x , 0) = 0, x ∈ ω.
Therefore, applying assumption ( H) and the fact that F ∈ C([0, T ]; H 2γ ( Ω)), we deduce that v ∈ C([0, T ]; H 2(1+γ) ( Ω)). Thus, in order to complete the proof, we only need to check that w ∈ C([0, T ]; H 2(1+γ) ( Ω)).

For this purpose, using the fact that F n ∈ L p (0, T ) and t → t α-1 ∈ L p (0, T ), with 1/p = 1 -1/p, one can check that w n ∈ C([0, T ]), n ∈ N. Moreover, fixing m, n ∈ N, with m < n, and applying formula (1.148) of [START_REF] Podlubny | Fractional differential equations[END_REF]Theorem 1.6], we find

n k=m λ 1+γ k w k (t)ϕ k L 2 ( Ω) t 0 n k=m λ γ k E α,1 (-λ k (t -s) α )F k (s)ϕ k L 2 ( Ω) ds C t 0 (t -s) -αγ n k=m F k (s)ϕ k L 2 ( Ω)
ds.

On the other hand, since 1 p < 1 -αγ, we have 1 p = 1 -1 p > αγ and we deduce T 0 s -p αγ ds < ∞.

Therefore, applying the Hölder inequality, we obtain

n k=m λ 1+γ k w k (t)ϕ k L 2 ( Ω) C n k=m F k (s)ϕ k L p (0,T ;L 2 ( Ω))
.

Therefore, we deduce that ũ ∈ C([0, T ]; H 3 (Ω)) ∩ C 1 ([0, T ]; H 2γ (Ω)) and, thanks to (1.9), ũ solves

                   ∂ α t ũ + Aũ + q(t, x )ũ = f (t, x ) R(t, x), (t, x , x n ) ∈ (0, T ) × ω × (-, ), ũ(t, x) = 0, (t, x) ∈ (0, T ) × ∂ω × (-, ), ũ(t, x , ± ) = 0 (t, x ) ∈ (0, T ) × ω, ∂ k t ũ(0, x) = 0, x ∈ ω × (-, ), k = 0, . . . , m α ,
where R is the extension of R to (0, T ) × ω × (-, ) given by

R(t, x , x d ) := R(t, x , -x d ), (t, x , x d ) ∈ (0, T ) × ω × (-, 0).
Fixing v = ∂ x d ũ and using the fact that ũ ∈ C([0, T ];

H 3 (Ω)) ∩ C 1 ([0, T ]; H 2γ (Ω)), we deduce that v solves the problem                    ∂ α t v + Av + q(t, x )v = H(t, x , x d ) + f (t, x )∂ x d R(t, x), (t, x , x n ) ∈ (0, T ) × ω × (-, ), v(t, x) = 0, (t, x) ∈ (0, T ) × ∂ω × (-, ), v(t, x , ± ) = ∂ x d ũ(t, x , ± ) = ±∂ x d u(t, x , ), (t, x ) ∈ (0, T ) × ω, ∂ k t v(0, x) = 0, x ∈ ω × (-, ), k = 0, . . . , m α , (3.1)
with

H(t, x , x d ) := d i,j=1 ∂ xi ∂ x d a ij (x)∂ xj ũ .
Note that here since R, ∂

x d R ∈ L ∞ ((0, T ) × Ω), we have R, ∂ x d R ∈ L ∞ ((0, T ) × Ω).
We are now in position to prove Theorem 1.4.

Proof of Theorem 1.4. In all this proof C denotes a generic constant depending on α, c, T , Ω, A, q, d. According to Lemma 1.3 the solution u of (1.11) is lying in

C([0, T ]; H 2(1+γ) ( Ω))∩C 1 ([0, T ]; H γ ( Ω))
and v = ∂ x d ũ, with ũ the even extension of u to (0, T ) × ω × (-, ), which solves (3.1). Note first that projecting the equation (1.11) in (0, T ) × ω × { } and using the fact that, for all (t, x ) ∈ (0, T ) × ω, we have u(t, x , ) = 0, we deduce that

f (t, x )R(t, x , ) = [a d,d ∂ 2 x d u + 2 d-1 j=1 a jd ∂ xj ∂ x d u + d j=1 ∂ xj a jd ∂ x d u](t, x , ), (t, x ) ∈ (0, T ) × ω.
Combining this with (1.13) and the fact that v = ∂ x d u, for all t ∈ (0, T ), we deduce that

f (t, •) L 2 (ω) c -1 C( ∂ x d v(t, •, ) L 2 (ω) + ∂ x d u(t, •, ) H 1 (ω) ) C( ∂ x d v(t, •, ) L 2 (ω) + ∂ x d u(t, •, ) H 3 2 (ω)
).

(3.2)

Thus, the proof of (1.17) will be completed if we can derive a suitable estimate of ∂

x d v(t, •, ) L 2 (ω)
For this purpose, we decompose v into

v = v 1 + v 2 , where v 1 solves                    ∂ α t v 1 + Av 1 = 0, in (0, T ) × ω × (-, ), v 1 (t, x) = 0, (t, x) ∈ (0, T ) × ∂ω × (-, ), v 1 (t, x , ± ) = ∂ x d ũ(t, x , ± ) = ±∂ x d u(t, x , ), (t, x ) ∈ (0, T ) × ω, v 1 (0, x) = 0, x ∈ ω × (-, ), (3.3)
and v 2 solves

                   ∂ α t v 2 + Av 2 = -qv(t, x) + H(t, x , x d ) + f (t, x )∂ x d R(t, x), (t, x , x n ) ∈ (0, T ) × ω × (-, ), v 2 (t, x) = 0, (t, x) ∈ (0, T ) × ∂ω × (-, ), v 2 (t, x , ± ) = ∂ x d ũ(t, x , ± ) = 0, (t, x ) ∈ (0, T ) × ω, v 2 (0, x) = 0, x ∈ ω × (-, ). (3.4) Since u ∈ C 1 ([0, T ]; H 2γ ( Ω)) ∩ C([0, T ]; H 2(1+γ) ( Ω)), we have (t, x ) → ∂ x d u(t, x , ) ∈ C 1 ([0, T ]; H 2γ-3 2 (ω)) ∩ C([0, T ]; H 1 2 +2γ (ω))
and ∂ x d u(0, •, ) = 0. therefore, using a classical lifting argument (e.g. [24, Theorem 8.3, Chapter 1]),

we can find G ∈ C 1 ([0, T ]; H 2γ-1 (Ω)) ∩ C([0, T ]; H 2γ+1 (Ω)) such that G(t, x , ) = ∂ x d u(t, x , ), (t, x ) ∈ (0, T ) × ω, G(0, x , x d ) = 0, (x , x d ) ∈ ω × (-, ), G W 1,∞ (0,T ;L 2 (ω×(-, ))) + G L ∞ (0,T ;H 2 (ω×(-, ))) C( ∂ x d u(•, •, ) L ∞ (0,T ;H 3 2 (ω)) + ∂ x d u(•, •, ) W 1,∞ (0,T ;H δ (ω)) ). (3.5)
Therefore, we can decompose

v 1 into v 1 = G + w 1 with w 1 solving                    ∂ α t w 1 + Aw 1 = G 1 , in (0, T ) × ω × (-, ), w 1 (t, x) = 0, (t, x) ∈ (0, T ) × ∂ω × (-, ), w 1 (t, x , ± ) = 0, (t, x ) ∈ (0, T ) × ω, w 1 (0, x) = 0, x ∈ ω × (-, ),
where

G 1 = -∂ α t G -AG ∈ C([0, T ]; L 2 (Ω))
. Thus, we have

w 1 (t, •) = t 0 S(s)G 1 (t -s)ds,
where S(t) corresponds to the operator valued function defined in Theorem 1.2. Therefore, applying (H00) and formula (1.148) of [29, Theorem 1.6], we deduce that w 1 ∈ C([0, T ]; H 2γ (Ω)) with

w 1 L ∞ (0,T ;H 2γ (Ω)) C T 0 s α(1-γ)-1 ds G 1 L ∞ (0,T ;L 2 (Ω)) . It follows that v 1 L ∞ (0,T ;H 2γ (Ω)) C( G 1 L ∞ (0,T ;L 2 (Ω)) + G L ∞ (0,T ;H 2 (Ω)) ) C( G W 1,∞ (0,T ;L 2 (ω×(-, ))) + G L ∞ (0,T ;H 2 (ω×(-, ))) )
and combining this with (3.5), we get

v 1 L ∞ (0,T ;H 2γ (ω×(-, )) C( ∂ x d u(•, •, ) L ∞ (0,T ;H 3 2 (ω)) + ∂ x d u(•, •, ) W 1,∞ (0,T ;H δ (ω)) ). (3.6)
Moreover, the estimate

∂ x d v 1 (•, •, ) L ∞ (0,T ;L 2 (ω)) C v 1 L ∞ (0,T ;H 2γ (ω×(-, )) , implies ∂ x d v 1 (•, •, ) L ∞ (0,T ;L 2 (ω)) C( ∂ x d u(•, •, ) L ∞ (0,T ;H 3 2 (ω)) + ∂ x d u(•, •, ) W 1,∞ (0,T ;H δ (ω)) ). (3.7)
On the other hand, we have

v 2 (t, •) = t 0 S(t -s)[-qv(s, •) + H(s, •) + ∂ x d F (s, •)]ds.
Thus, applying (H00) and repeating the arguments of Lemma 1.3, we get

v 2 (t, •) H 2γ (ω×(-, )) C t 0 (t-s) α(1-γ)-1 [ f (s, •) L 2 (ω) + v(s, •) L 2 (ω×(-, )) + H(s, •) L 2 (ω×(-, )) ]ds
which, for all t ∈ (0, T ], implies that

∂ x d v 2 (t, •, ) L 2 (ω) C t 0 (t-s) α(1-γ)-1 [ f (s, •) L 2 (ω) + v(s, •) L 2 (ω×(-, )) + H(s, •) L 2 (ω×(-, )) ]ds.
In light of (1. Proof of Corollary 1.5. Let u = v 1 -v 2 and notice that u solves (1.2) with q = q 1 , f = (q 2 -q 1 ) and R = v 2 . Then, using the fact that v 2 , ∂ x d v 2 ∈ C([0, T ] × Ω) and the fact that, thanks to (1.13),

(1.20), we are in position to apply Theorem 1.4 from which we deduce (1.20).

Appendix

We consider the following result which can be deduced from other known results (see e.g [30,

Theorem 2.2] and [START_REF] Li | Initial-boundary value problem for distributed order time-fractional diffusion equations[END_REF]Theorem 1.3]) considered for α ∈ (0, 1) that we extend to α ∈ (0, 2).

∂ α t u L 2 (Q) + Au L 2 (Q) + u L 2 (Q) C F L 2 (Q) . (3.8)
Proof. We prove this result for sake of completeness. We fix A the operator A acting on L 2 (Ω) with the boundary condition (1.3)- (1.4). We fix also the non-decreasing sequence of non-negative eigenvalues with G k (t) = t α-1 E α,α (-λ k t α ). Then, we find

T 0 ∞ k=1 λ 2 k |u k (t)| 2 dt = ∞ k=1 λ 2 k (G k 1 (0,T ) ) * (F k 1 (0,T ) ) 2 L 2 (0,T )
and an application of the Young inequality yields

T 0 ∞ k=1 λ 2 k |u k (t)| 2 dt = ∞ k=1 λ 2 k G k 2 L 1 (0,T ) F k 2 L 2 (0,T ) .
On the other hand, we have In the same way, we find

∂ α t u = -Au + F ∈ L 2 (Q)
which implies at the same time that ∂ α t u ∈ L 2 (Q) and (3.8). This proves (3.8).

1 . Introduction 1 . 1 .

 111 Statement. Let d 2, Ω = ω × (-, ) and Ω = ω × (0, ), and ω ⊂ R d-1 be a bounded domain with C 2 boundary. We set Q = (0, T ) × Ω and Q = (0, T ) × Ω. Let ν = ν(x) be the outward unit normal vector to ∂Ω or ∂ Ω at x. In what follows, we define A by the differential operator Au(x) = -d i,j=1

1. 7 .

 7 Comments about our results. To our best knowledge Theorem 1.2 and 1.4 are the first results of recovery of a source term depending on both time and space variables for fractional diffusion equations of the form (1.2) when α = 1. For α = 1, we refer to[START_REF] Isakov | Inverse source problems[END_REF] Section 6.3] addressing this inverse problem

3 2

 3 (Ω) ), t ∈ (0, T ].In view of the equation satisfied by v and according to the above arguments as well as the arguments used in Proposition 1.3, for all t ∈ (0, T ], we getv(t, •) H 1 (Ω) + u(t, •) H 1 (Ω) C( f L ∞ (0,t;L 2 (ω)) + ∂ x d u(•, •, ) L ∞ (0,T ;H (ω)) + ∂ x d u(•, •, ) W 1,∞ (0,T ;H δ (ω)) ).Thus, we findH(t, •) L 2 (Ω) C( f L ∞ (0,t;L 2 (ω)) + ∂ x d u(•, •, ) L ∞ (0,T ;Hand repeating the arguments of the preceding section, we can prove that v ∈ C([0, T ]; H 2(1+γ) ( Ω)) and for d 3, the Sobolev embedding theorem implies that v, ∂ x d v ∈ C [0, T ] × Ω . Applying the previous results about recovery of source terms we can complete the proof of Corollary 1.5.

Lemma 3 . 1 .

 31 Let F ∈ L 2 (Q), α ∈ (0, 2) and m, n = 0, 1. Then problem (1.2)-(1.4) admits a unique solution u ∈ L 2 (0, T ; H 1 (Ω)), satisfying Au, ∂ α t u ∈ L 2 (Q)and the following estimate holds true

(λ k ) k 1

 1 and associated eigenfunctions (ϕ k ) k 1 of A. Then, we considerk 0 := min{k ∈ N : λ k > 0}.Since u solves (1.2)-(1.4), we haveu k (t) := u(t, •), ϕ k L 2 (Ω) = t 0 (t -s) α-1 E α,α (-λ k (t -s) α )F k (s)ds, where F k (t) := F (t, •), ϕ k L 2 (Ω) . Therefore, we have u k (t) = (G k 1 (0,T ) ) * (F k 1 (0,T ) )

G k L 1 sλ 2 k

 12 (0,T ) = T 0 s α-1 E α,α (-λ k s α )ds = λ -1 α-1 E α,α (-s α )ds, k k 0 .Combining this with formula (1.148) of [29, Theorem 1.6], we getG k L 1 (0,T ) Cλ -1 )|u k (t)| 2 dt = C u L 2 (0,T ;D(A)) C F L 2 (Q) .

  [START_REF] Jiang | Weak unique continuation property and a related inverse source problem for time-fractional diffusion-advection equations[END_REF], we getH = ∂ x d a dd ∂ 2 x d ũ + 2 ∂ x d a dd ∂ x d v + 2 •) L 2 (Ω) C( v(t, •) H 1 (Ω) + u(t, •) H 1

	d-1						d-1
								∂ 2 x d a jd ∂ xj	ũ
							j=1
	= d-1	∂ x d a jd ∂ xj v +	 	d	∂ xj ∂ x d a jd	  v +	d-1	∂ 2 x d a jd ∂ xj	ũ
	j=1			j=1			j=1
	and it follows that
	H(t,

j=1 ∂ x d a jd ∂ xj ∂ x d ũ + d j=1 ∂ xj ∂ x d a jd ∂ x d ũ +

(ω)) + ∂ x d u(•, •, ) W 1,∞ (0,T ;H δ (ω)) )

Acknowledgements.

The work of the first author is partially supported by the French National Research Agency ANR (project MultiOnde) grant ANR-17-CE40-0029. The second author is supported by Grant-in-Aid for Scientific Research (S) 15H05740 of Japan Society for the Promotion of Science and by The National Natural Science Foundation of China (no. 11771270, 91730303), and the "RUDN University Program 5-100".

Combining this with the fact that F ∈ L p (0, T ; L 2 ( Ω)), we deduce that

Here we have used the fact that F |t=0 = 0. Repeating the above arguments and using the fact that

Now let us consider the case q = 0. We introduce the map

Then, using a classical fixed point argument combined with the preceding analysis, we deduce that G admits a unique fixed point lying in C([0, T ]; H 2(1+γ) ( Ω)) ∩ C 1 ([0, T ]; H 2γ ( Ω)) which will be the solution of (1.11). This completes the proof of the lemma.

From now on and in all the remaining part of this section, we assume that α ∈ (0, 1] and that the conditions of Proposition 1.3 are fulfilled. We consider ũ defined on (0, T ) × ω × (-, ) by

and it follows that

Combining this with (3.2) and (3.7), we find

Then, [7, Lemma 3] (see also [33, Theorem A]), implies that

from which we deduce (1.17).

3.1. Application to the recovery of coefficients. Consider v the solution of the problem

with w 0 , h k , k = 0, 1, given by (1.19). Then, we can write v = H + y where y solves

y(t, x) = 0, (t, x) ∈ (0, T ) × ∂ω × (0, ), y(t, x , ) = 0 (t, x ) ∈ (0, T ) × ω, ∂ x d y(t, x , 0) = 0 (t, x ) ∈ (0, T ) × ω, y(0, x) = 0, x ∈ Ω,

where H 1 = -(∂ α t + A + q)H ∈ C([0, T ]; H 2γ ( Ω)) ∩ W 1,p (0, T ; L 2 ( Ω)). Therefore, using the fact that, thanks to (1.19), we have H 1 (0, x) = 0, x ∈ Ω


X Bednarek

S Martin

A Ndiaye

V Peres

O Bonnefoy

Extrapolation of DEM simulations to large time scale. Application to the mixing of powder in a conical screw mixer

Keywords: Mixing powders, DEM, Large time scale simulation, Extrapolation

The paper proposes an original algorithm which allows a long time scale extrapolation of DEM results at a very low computational cost. This algorithm can be adapted to any periodic processes. In this study, it is applied to the mixing process of powders within a conical screw mixer. The results are then compared with long time DEM simulations. It appears that this method is able to predict the DEM results with a very good accuracy.

Introduction

Over the last decades, the Discrete Element Method (DEM) [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] has become a major tool to simulate the behavior of granular media. Since it is the only method able to consider the discrete aspect of materials, it is well designed to represent the mixing of powders. Thus, sev-ent adjustable parameters, like diffusion coefficient and segregation velocity, are fitted from DEM simulations. For more information on this approach, see the following references: [START_REF] Fan | Modelling size segregation of granular materials: the roles of segregation, advection and diffusion[END_REF]; [START_REF] Schlick | Modeling segregation of bidisperse granular materials using physical control parameters in the quasi-2D bounded heap[END_REF][START_REF] Schlick | A continuum approach for predicting segregation in flowing polydisperse granular materials[END_REF]. If this approach seems very promising, it has only been applied to very simple system so far.

In this context, the simulation of monodisperse powders remains very interesting. It does not exhibit the full complexity of industrial processes but allows the implementation of simple numerical approaches which can give very relevant results. For example, [START_REF] Doucet | Modeling of the mixing of monodisperse particles using a stationary DEM-based Markov process[END_REF] proposed to simulate the mixing of powders with a DEM-based stationary Markov process, which permits to extrapolate the DEM results to large time scale. This paper an original algorithm which allows a long time scale extrapolation of DEM results at a very low computational cost. Contrary to the works of [START_REF] Doucet | Modeling of the mixing of monodisperse particles using a stationary DEM-based Markov process[END_REF], the present method is fully deterministic and meshless. This original algorithm can be adapted to any periodical process with monodisperse powder. It can also be implemented in coarse grain simulations to scale up the simulation time. In this paper, the method is applied to the simulation of a Conical Screw Mixer with one million particles.

Material and method

The method is very simple and can be generalized to many different systems. However, its implementation requires the following conditions:

• The targeted process is periodic,

• The particles are identical,

• The granular medium is dense -i.e. the distance between the centers of neighboring particles is roughly equal to the diameter of a single particle.

Note that these two last conditions are soft constraints because they do not prevent the algorithm to be executed. Not fulfilling the second condition would lead to less accurate results.The method has indeed be designed assuming that all the particles have a similar behavior. Moreover, for different sized particles, segregation would not be taken into account. Then, not fulfilling the last condition would lead to minor errors. In the case of our screw mixer, the medium is very dense and this point is not an issue. However, the error could become significant for applications like spouted beds for which the bulk density is more heterogeneous.

The general idea of the method is first to run a short DEM simulation over a time-lapse which is greater than the period T characteristic of the process. It should be noticed that T is very different from the time step of the DEM simulation. T is indeed only a property of the process, it is independent of the simulation parameters. It is generally several order of magnitude greater than the DEM resolution time step. For a mixing process with a rotating impeller, this time-lapse T should be the rotation period of the screw impeller, while for a rotating drum this could be the interval between two avalanches. We note α≥1 the number of period T which is simulated:

t 1 -t 0 = αT (1)
where t 0 is the initial time fixed by the user and t 1 is the final time of the DEM simulation. Then, the purpose of the extrapolation is simply to repeat the track followed by particles between t 0 and t 1 for predicting t 2 and then t 3 , t 4 ...

The repetition of the trajectories must be done on a single structure: that of the state t 1 . One calls "structure" the set of the positions of particles for a given time. Since the structure at t 0 is not the same as that at t 1 , it is necessary to bring the first one closer to the second. Thus, the crucial step of the algorithm consists in pairing the particles: each particle j among the N particles of the simulation at the starting time t 0 , at position -→ r j (t 0 ), is associated with a particle i of the ending time t 1 , at position -→ r i (t 1 ), by finding the twinning function f , defined by Equation 2, which is bijective.

f : 1, N -→ 1, N i -→ j = f (i) (2)
This twinning function has to be defined in such a way that -→ r f (i) (t 0 ) is sufficiently close to -→ r i (t 1 ). The meaning of "sufficiently" is discussed later in the text. A graphical representation of the pairingie twinning -step is given in Figure 1. Note that the goal is to match particles at two distinct times.

The fact that the trajectories are repeated implies the hypothesis that the paired particles must have the same behavior and therefore be mechanically identical. So in the ideal case they differ at most only in their color.

Finally, the position of the particle i at time t 2 = t 1 + αT can be easily extrapolated:

- → r i (t 2 ) = - → r f (i) (t 1
). In other words, particle i will move, at time t 2 , to the position of the particle

f (i) at time t 1 . Let f n = f • . . . • f
be the composition of n functions f , and n a strictly positive integer. The position of particle i at the time t n+1 = t 1 + n αT is given by Equation 3.

- → r i (t n+1 ) = - → r f n (i) (t 1 ) (3)
Thus, knowing f , the extrapolation step can be reproduced many times at a very low computational cost. The general position vector and the displacement associated with the position during the time interval αT are indeed fixed over time.

Example

In order to make the method clearer, an example of the extrapolation approach is presented bellow. Lets consider a nondescript 9 particles system.

1. First of all, a DEM simulation is run. The initial steady state time t 0 is fixed and represented in the left picture of Figure 1. Then, the simulation is run up to t 1 = t 0 + αT . The position of particle at time t 1 is shown in the middle picture of Figure 1.

2. Then, the twining function has to be determined. Looking at the positions of particles at time t 0 and t 1 , it can be observed that particle 3 at time t 1 is closed to the position that was occupied by particle 1 at t 0 .

Thus f (3) = 1. In the same way one has f (1) = 9, f (2) = 7, f (4) = 8,

f (5) = 6, f (6) = 2, f (7) = 3, f (8) = 5 and f (9) = 4.
3. Finally, the positions of every particle can be extrapolated for the following times. Applying Equation 3, the position of particle 1 at time t 2 is predicted:

- → r 1 (t 2 ) = - → r f (1) (t 1
) that is -→ r 9 (t 1 ). At time t 3 , as one has f (f (1)) = 4, the position of particle 1 is predicted at -→ r 4 (t 1 ). Finally, since f (f (f (1))) = 8, the position of particle 1 at time t 4 is predicted at -→ r 8 (t 1 ). Figure 2 shows these three extrapolated steps, particle 1 being in green, and follows the paths presented in Figure 1, moving the green color on the structure of state t 1 . The same procedure is applied to all other particles.

The most time consuming part of the method is the pairing algorithm.

It is described in the following section.

Pairing algorithm

The principle of the pairing algorithm is to couple particles of time t 0 with particles of time t 1 in order to get a set of the closest pairs. The aim is then to find f that minimizes the objective function F obj defined as follows:

F obj = N i=1 - → r i (t 1 ) -- → r f (i) (t 0 ) (4)
The optimum pairing solution f corresponds to the minimum of the sum of the N center to center distances d pair i = -→ r i (t 1 ) --→ r f (i) (t 0 ) between pairs of particles, N being the number of particles in the simulation.

Since each particle has to be associated with another single particle, the full direct resolution of this problem would require the evaluation of the N ! possible combinations. Hopefully, such a problem can be viewed as an Euclidean assignment problem [START_REF] Rendl | On the Euclidean assignment problem[END_REF] and can be solved numerically by means of the Hungarian method introduced by [START_REF] Kuhn | The Hungarian method for the assignment problem[END_REF]. The input of this algorithm is the table of distances between each pair of particles candidate to coupling, whose size is N 2 . This leads to a prohibitive computational memory and time calculation for large scale simulations with N > 10 5

particles. Therefore, the number of combinations have to be substantially decreased. This aim is reached by introducing two simplifications:

1. The method can be improved by avoiding taking into account pairs of particles, the distance between which is too high. Thus, the first simplification step consists in allowing a maximum number N c of pairing possibilities for each particle, excluding the furthest ones. This can be seen as a cutoff distance, similar to those implemented in the SPH or Molecular Dynamics methods, except that, in the present case, the cutoff criterion is firstly based on a maximum number of neighoring particles. Each particle i at step t 1 has to get N c candidates j at step t 0 and vice versa. Hence, the number of candidates is at least equal to N c . Further explanations are given in Figure 3.

2. Simulation domain is divided into N x × N y × N z subdomains where N x , N y , and N z are respectively the number of slices in the x, y, and z directions. The Hungarian algorithm is then run separately on every subdomain. This is kind of a domain decomposition which simplifies the resolution and also allows its parallel implementation.

In this study, with N = 10 6 particles, the pairing algorithm has been run

with N c = 300 and

N x = N y = N z = 10.
It should be noticed that both the limitation to N c pairing candidates for each particle and the division of the global domain might not allow the convergence to the minimum for the objective function. Thus, the quality of the pairing solution should be evaluated before any extrapolation of the DEM results.

DEM model

A classical DEM model is used in this study. The particles are spherical and the particle size distribution (PSD) is uniform with a span (

d max p -d min p )
equal to 5% of the mean diameter d p . Since the purpose is to give an example of the extrapolation procedure, the contact law used here -the non-linear Hertz-Mindlin model -is not detailed but can be found elsewhere [START_REF] Zhu | Discrete particle simulation of particulate systems: Theoretical developments[END_REF]. The simulation has been run with LIGGGHTS 3.1.0 [START_REF] Kloss | Models, algorithms and validation for opensource DEM and CFD-DEM[END_REF]. The parameters used in the DEM simulation are given in Table 1.

Description of the conical screw mixer

This method has been applied to the DEM simulation of powder mixing within a conical screw mixer. The later is a conic blender with a screw

impeller. An accurate description of the mixer is given in Figure 4. The screw impeller rotates around the axis of the cone -→ e y , going along the wall with an angular velocity ω a . This rotation will be labeled as primary rotation.

The screw impeller also rotates around its own axis with a velocity ω s . This will be labeled as secondary rotation. In this study, ω s = 30 ω a .

Filling procedure

The conical vessel is filled with particles using a classical DEM procedure: initial positions of particles are randomly chosen in the upper part of the cone and particles are then allowed to fall by gravity. The filling step is considered to be accomplished when the mean translational velocity is bellow 10 -4 m.s -1 .

Friction coefficients and rolling friction coefficients are set to zero to ensure that the upper free surface is almost flat at the beginning of the mixing process.

I L = σ 2 0 -σ 2 m σ 2 0 -σ 2 r (5)
where σ 2 m is the variance of the composition of the component of interest over the sampling, σ 2 0 is the variance of the fully segregated state: σ 2 0 = p(1p), with p the ratio of particles of the component of interest in the whole mixture, σ 2 r the variance of the random system: σ 2 r = p(1-p) ns with n s the number of particles -regardless the nature of the particle -in each sample.

This index tends rapidly to 1 with respect to the mixing quality. Hence, we chose to also calculate the Ashton index I A [START_REF] Ashton | Mixing of powders and particles in industrial mixers[END_REF] which allows a better discrimination between mixing states when the system is almost homogeneous:

I A = log σ 2 0 σ 2 m log σ 2 0 σ 2 r (6)
This index also tends to 1 but slower to the Lacey index. It should be noticed that both indices depend strongly on the sampling protocol. Therefore, in this paper, the sampling is based on a 3D Voronoï diagram: 1000 particles are chosen randomly and 1000 samples are built around this particles. Each of the particles of the mixture is attached to the sample corresponding to the particle -among the 1000 initially selected -which is closest to it. This procedure is repeated 10 times and the indices presented are the average of the 10 mixing values thus obtained.

Dispersion

The dispersion is a mixing criterion that takes into account both convective and diffusive mixing mechanisms. In order to quantify it, N clusters are considered. Each cluster i is composed of a central particle i and neighbouring particles j. The neighbourhood limited to particles j that are separated from the central particle by a surface-to-surface distance of less than a mean radius i.e. R = 8.3 10 -1 mm. This translates into:

- → r j -- → r i < R i + R j + R (7)
Each particle is often part of several clusters. The distance of dispersion d i disp of one single cluster i is defined as the average distance between the central particle i and its M i neighbouring particles j at a time t minus the same mean distance but ∆t disp earlier.

d i disp (t; ∆t disp ) = M i j=1 - → r j (t) -- → r i (t) -- → r j (t -∆t disp ) -- → r i (t -∆t disp ) (8)
Most of the time, such a distance is positive but it can also be negative if any compression of the cluster occurs. [START_REF] Golshan | Granular mixing in nauta blenders[END_REF] used a similar criteria to define their tracers. Here, every particle is considered as the center of a cluster.

Finally, the dispersion d disp is interpolated to a tetrahedral mesh. This allows us to build a continuous dispersion field, for a given duration ∆t disp , over the volume of the conical screw mixer.

Results and discussion

The extrapolation method decribed above has been applied to the mixing process of powder.

The standard DEM simulation has been run for a full primary revolution of the screw impeller in the conical mixer, which corresponds to 30 s of physical time and 30 rotations of the screw impeller around its own axis (with

ω s = 60 rpm).
During post-processing, four initial configurations were created by assigning different colors to the particles regarding their initial positions. The different configurations are labeled according to the number of slices in the corresponding direction, respectively 2x, 2y, 2z and 8y. For example, 8y means that initialy there are eight layers of particles stacked on top of each other in the y direction. These configurations are represented in Figure 5.

Since only one type of particle is mixed, mixing indices cannot be calculated without any coloration procedure.

Finally, the extrapolation algorithm has been run to predict the two next primary revolutions of the screw impeller.

Determining the steady state

Extrapolation requires that the system be in a permanent state. It means that the velocity field has to reach its periodic state with respect to the screw rotation. After the filling step, the simulation starts with an acceleration step, which is due to the inertia of the system. Moreover, the shape of the upper free surface of the granular packing can evolve slightly between the static and dynamic state. This explains why applying the pairing algorithm before reaching the steady state would increase significantly the error of the pairing solution, particularly near the top of the mixer.

Therefore, the beginning steady state can be determined by considering the average velocity of the particles and the torque acting on the screw impeller, along its axis of rotation. It is shown in Figure 6 that the steady state is reached after ten secondary revolution of the impeller i.e. 10 s of mixing. In this study, the initial state was set at t 0 = 21 s for the extrapolation algorithm.

One natural time period of the process is T = 1 s, the period of rotation of the screw. We take α = 9 in Equation 1 and the final state of the simulation is at t 1 = 30 s. The extrapolation is so conducted with a αT = 9 s step which corresponds to 9 full revolutions of the screw impeller.

Accuracy of the pairing algorithm

If the pairing algorithm is executed directly on steps t 0 and t 1 , this will lead to poor results because of the non alignment of the empty volumes created by the screw impellers that are not in the same position for the two time steps. In that case, it is a critical point to superimpose the screws before the execution of the algorithm. Thus the whole state at time t 0 is rotated of an angle of αT ω a = -108 • around -→ e yi.e. all particles and the meshes. The extrapolated time step t n+1 is then obtained from extrapolated data after rotation by an angle -n × 108 • around -→ e y .

The efficiency of the pairing algorithm can be evaluated by looking at the distances between paired particles. The most accurate results correspond to the smallest distances. The pairing distance distribution is given in Figure 7 and Table 3. It appears that the distance d pair between paired particles is less than the mean diameter d p (resp. 1.2 d p ) for 90% (resp. 95%) of the particles. This result is considered as rather good. However, it should be noticed that for a very few particles, the pairing distance increases up to 27.7 times the mean diameter d p . This can be explained by the slight evolution of the free surface shape between t 0 and t 1 and the domain decomposition into N x × N y × N z = 1000 domains in order to accelerate the execution of the pairing algorithm. This can be seen in Figure 8. It is also interesting to notice that this algorithm works well in our dense packing simulation, but would lead to a higher error for looser packings. The pairing algorithm was run on a single computer with a 4 × 3.5 GHz intel i5 processor and took about 10 h for 10 6 particles.

Extrapolation of the DEM results

The extrapolation algorithm has been run to predict the two following primary revolutions of the screw impeller. To validate the extrapolation results, a DEM simulation has also been run for two extra full primary revolutions.

The time gains over the full simulations are summed up in Table 2. The later presents typical CPU times for full DEM simulation compared with the extrapolation protocole, and gives the ratio of CPU time between DEM simulation and extrapolation with respect to the number of impeller primary convolution. Hence, the last column represents the gain factor which is expected. Since the extrapolation process requires a short initial DEM simulation and then the resolution of the pairing algorithm, this method becomes very interesting for long time extrapolation. Here, it has been limited

to two extra revolutions because it was not possible to run longer DEM simulation for comparison. Thus, once the pairing algorithm has been solved, the extrapolation step for a full primary revolution takes about 1 min while it takes ≈ 5. × 10 3 hr of CPU time for DEM. This emphasizes the major gain reached by the extrapolation approach.

However, attention should be drawn to the fact that those results give only orders of magnitude. They could evolve slightly according to the computer configuration and optimization of the codes. In our study, DEM simulation were run with a 2.6 GHz cluster with 40 parallel threads while both pairing and extrapolation algorithms were run on a single 4 × 3.5 GHz intel i5 CPU.

For comparison, the cross sectional views of DEM and extrapolated results are given in Figure 9 and Figure 10 for the different configurations.

All the results for extrapolation and DEM look very similar. For example,

Table 2: Comparison between typical CPU times (hr) between DEM simulation and extrapolation method with respect to the number of impeller primary convolution. On the last line, +1 represent an extra primary convolution once the pairing algorithm has been resolved. For example, the CPU time for 3 full primary convolutions, is equal to the CPU time for 2 full convolutions plus the time for 1 extra primary convolution.

Primary DEM (hr) Extrapolation method (hr) Ratio convolution Initial DEM Pairing Extrapolation total

1 5 × 10 3 5 × 10 3 0 0 5 × 10 3 1 2 1 × 10 4 5 × 10 3 4 × 10 1 1.5 × 10 -2 5.04 × 10 3 2 +1 5 × 10 3 1.5 × 10 -2 3 × 10 5
the shape of the red layer of the 8y configuration is well reproduced at the different stages of the simulation.

In order to get a more accurate comparison between extrapolation and DEM simulations results, Lacey and Ashton indices have been calculated for the four different initial configurations.

The results are presented repectively in Figure 11 and Figure 12 for Lacey indices and in Figure 13 and Figure 14 for Ashton indices. Those make it clear that the mixing kinetic is correctly predicted by the extrapolation scheme. To go further and evaluate the accuracy of the method, the relative difference, between extrapolated and DEM results for the Lacey and Ashton indices is given in Figure 15 and Figure 16 respectively. It appears that the error remains lower than 0.7% for the Lacey index and lower than 1.5% for the Ashton index, for all the configurations. This shows that the extrapolation algorithm reproduces very accurately the different mixing mechanisms.

Finally, the dispersion field was calculated for ∆t disp = αT = 9 s, start-good agreement between extrapolated and DEM results. Thus, this method is able to predict accurately long time results at a very low computational cost compare to a full DEM simulation. This method has been applied to the simulation of mixing in a conical screw mixer but it could easily be applied to any other periodic mixers like rotating drums, four bladed mixers, Turbula R , V-blenders... The method could also be extended to simulate more complex processes, like spouted beds, the particles of which have a cyclic behavior.

2. 5 .

 5 Evaluation of the mixing quality 2.5.1. Mixing indices Multiple indices are available in the literature to evaluate the mixing quality for at least two components: the component "of interest" and the other one that can regroup several component excluding of course the one "of interest". The Lacey index (Lacey, 1954) I L is the most common. It is defined as follows:

Figure 1 :

 1 Figure 1: Representation of the pairing algorithm: particles position at time t 0 (left), t 1 (middle) and pairing step (right) where pairing distances are represented by doubleheaded arrows. Arrows (middle) arrows show the trajectories that have to be repeated by the algorithm.

Figure 2 :

 2 Figure 2: Representation of extrapolation algorithm, with the particle position at time t 2 (left), t 3 (middle) ans t 4 (right). Structure of state t 1 is conserved.

Figure 3 :Figure 4 :

 34 Figure 3: Example of pairing step (left) with N c = 2: searching N c candidates for particles of t 0 lead to the first table (middle) but particle i = 1 has only one candidate; adding one more candidate to i = 1 lead to the second table (right) with all particles having at least N c particles.

FigureFigure 6 :

 6 Figure 5: Representation of the four initial configurations

Figure 7 :

 7 Figure 7: Representation of the pairing distances with cumulative distribution function (CDF) for the N pair with an y-axis slightly dilated in order to get better indication on the highest distances.

Figure 8 :Figure 9 :

 89 Figure 8: Pairing algorithm efficiency, considering pre-classification into 10 × 10 × 10 domains (left), and regarding d pair repartition in 3D (right).

Figure 10 :Figure 11 :

 1011 Figure 10: Cross sectional views of the mixer at different stages for the 2z and 8y initial configurations. Comparison between DEM and extrapolated results.

Figure 12 :Figure 13 :Figure 14 :Figure 15 :Figure 16 :Figure 17 :Figure 18 :

 12131415161718 Figure 12: Comparison between extrapolated and DEM Lacey index I L . Results for configuration 8y.

Table 1 :

 1 DEM parameters

Table 3 :

 3 Pairing distances

	Frequency[%] Relative distance d pair / (2R)
	10	0.29
	50	0.55
	90	0.94
	95	1.18
	99	4.42
	99.9	9.85
	99.99	18.52
	100	27.70

Table 4 :

 4 Decile absolute relative dispersion error |(d extra

	disp -d DEM disp )/d DEM disp |

ing from t 1 , for DEM results and extrapolation results. This parameter is very interesting for two reasons: (1) the pairing algorithm adds an artificial dispersion by pairing the particles with particles that are not exactly at the same location and (2) it describe in details the diffusion mechanism. Thus, the comparison between the dispersion fields for DEM and extrapolated results is given in Figure 17. The relative error is presented in Figure 18 and Table 4. We can see that the relative differences can be high, greater than 1000%. Two reasons explain this observation. First, when DEM dispersion is near zero, the dispersion added by the extrapolation algorithm is comparatively high even it is if low in absolute value. Secondly, the highest pairing distances near the free surface lead to a poorer rendering of the dispersion magnitude as shown in Figure 17.

In order to investigate the influence of α, another extrapolation has been made with α = 1ie with an extrapolation step of 1 s. In that case, the relative errors are larger, about 5% percent for the Lacey index and 8% percent for the Ashton index. This is simply explained by the fact that the prediction of the same time t, needs 9 times more steps with α = 1 than with α = 9. This leads to greater errors comming from the dispersion error accumulated at each extrapolation step. Thus a larger α decreases the cumulative error but leads to less points to describe the evolution of the extrapolated mixture.

Conclusion

An original method was presented for the long time extrapolation of DEM simulation by using a pairing algorithm. The results show an extremely