Peter Lewintan

Patrizio Neff

L p -trace-free version of the generalized Korn inequality for incompatible tensor fields in arbitrary dimensions

Keywords: AMS 2010 subject classification: Primary: 35A23; Secondary: 35B45, 35Q74, 46E35 W 1, p (Curl)-Korn's inequality, Lions lemma, Nečas estimate, generalized Curl, incompatibility

For n ≥ 3 and 1 < p < ∞ we prove an L p -version of the generalized trace-free Korn-type inequality for incompatible, p-integrable tensor fields P : Ω → R n×n having p-integrable generalized Curln and generalized vanishing tangential trace P τ l = 0 on ∂Ω, denoting by {τ l } l=1,...,n-1 a moving tangent frame on ∂Ω. More precisely, there exists a constant c = c(n, p, Ω) such that

where the generalized Curln is given by (Curln P ) ijk := ∂iP kj -∂jP ki and devn X := X -1 n tr(X) • 1n denotes the deviatoric (trace-free) part of the square matrix X. The improvement towards the threedimensional case comes from a novel matrix representation of the generalized cross product.

Introduction

The estimate

∃ c > 0 ∀ u ∈ W 1, p 0 (Ω, R n ) : Du L p (Ω,R n×n ) ≤ c dev n sym Du L p (Ω,R n×n ) , (1.1)
for n ≥ 2 and p ∈ (1, ∞) where dev n X := X -1 n tr(X) • 1 denotes the deviatoric (trace-free) part of the square matrix X and its (compatible) generalizations of (1.1) are well known, cf. [START_REF] Dain | Generalized Korn's inequality and conformal Killing vectors[END_REF][START_REF] Fuchs | An application of a new coercive inequality to variational problems studied in general relativity and in Cosserat elasticity giving the smoothness of minimizers[END_REF][START_REF] Reshetnyak | Estimates for certain differential operators with finite-dimensional kernel[END_REF][START_REF] Schirra | New Korn-type inequalities and regularity of solutions to linear elliptic systems and anisotropic variational problems involving the trace-free part of the symmetric gradient[END_REF][START_REF] Smith | Formulas to represent functions by their derivatives[END_REF]. In [START_REF] Lewintan | L p -trace-free generalized Korn inequalities for incompatible tensor fields in three space dimensions[END_REF] another generalization to the (incompatible) case ∃ c > 0 ∀ P ∈ W 1, p 0 (Curl; Ω, R 3×3 ) :

P L p (Ω,R 3×3 ) ≤ c dev 3 sym P L p (Ω,R 3×3 ) + Curl P L p (Ω,R 3×3 ) (1.2)
has been given. The main objective of the present paper is to extend (1.2) to the trace-free case for n ≥ 3 dimensions. Such a result was expected, cf. [START_REF] Lewintan | L p -trace-free generalized Korn inequalities for incompatible tensor fields in three space dimensions[END_REF]Rem. 3.11] and was already proven to hold true for p = 2, cf. [START_REF] Bauer | Dev-Div-and DevSym-DevCurl-inequalities for incompatible square tensor fields with mixed boundary conditions[END_REF]. However, the latter used a Helmholtz decomposition and a Maxwell estimate and is not amenable to the L p -case. On the contrary, the argumentation scheme using the Lions lemma resp. Nečas estimate, known from classical Korn inequalities, turned out to be also fruitful in the case of Korn inequalities for incompatible tensor fields, cf. [START_REF] Lewintan | Korn inequalities for incompatible tensor fields in three space dimensions with conformally invariant dislocation energy[END_REF][START_REF] Lewintan | L p -trace-free generalized Korn inequalities for incompatible tensor fields in three space dimensions[END_REF][START_REF] Lewintan | L p -versions of generalized Korn inequalities for incompatible tensor fields in arbitrary dimensions with p-integrable exterior derivative[END_REF][START_REF] Lewintan | Nečas-Lions lemma revisited: An L p -version of the generalized Korn inequality for incompatible tensor fields[END_REF] and also [START_REF] Gmeineder | On Korn-Maxwell-Sobolev Inequalities[END_REF]. The secret of success is then to determine a linear combination of certain partial derivatives. One such expression in [START_REF] Lewintan | L p -trace-free generalized Korn inequalities for incompatible tensor fields in three space dimensions[END_REF] was D 2 (A + ζ • 1 3 ) = L(DCurl(A + ζ • 1 3 )) denoting by L a constant coefficients linear operator, for a skew-symmetric matrix field A and scalar field ζ. Here we catch up with a corresponding linear expression in all dimensions n ≥ 3. For that purpose a careful investigation of the generalized cross product, especially a corresponding matrix representation will be given. Indeed, it is this matrix representation which allows us to obtain suitable relations which are not easily visible in index notations. Korn's inequalities in higher dimensions for matrix-valued fields whose incompatibility is a bounded measure and corresponding rigidity estimates were obtained in the recent papers [START_REF] Conti | Sharp rigidity estimates for incompatible fields as consequence of the Bourgain Brezis div-curl result[END_REF][START_REF] Lauteri | Geometric rigidity estimates for incompatible fields in dimension ≥ 3[END_REF], however, without boundary conditions. More precisely, Conti and Garroni [START_REF] Conti | Sharp rigidity estimates for incompatible fields as consequence of the Bourgain Brezis div-curl result[END_REF] obtained as a consequence of a Hodge decomposition with critical integrability due to Bourgain and Brezis for

P ∈ L 1 (Ω, R n×n ) with Curl n P ∈ L 1 (Ω, R n× n(n-1)
2

) the sharp geometric rigidity estimate inf

R∈SO(n) P -R L 1 * (Ω) ≤ c dist(P, SO(n)) L 1 * (Ω) + Curl n P L 1 (Ω) (1.3)
with a constant c = c(n, Ω), the Sobolev-conjugate exponent 1 * := n n-1 and where the generalized Curl n is seen without a matrix representation as (Curl n P ) ijk := ∂ i P kj -∂ j P ki . Replacing the geometric rigidity by Korn's inequality they deduced from (1.3) furthermore inf A∈so(n)

P -A L 1 * (Ω) ≤ c sym P L 1 * (Ω) + Curl n P L 1 (Ω) .
(1.4)

This estimates remain true for Curl n P being a Radon measure. In that case the L 1 -norm of Curl n P has to be substituted by the total variation of the measure Curl n P , cf. [START_REF] Conti | Sharp rigidity estimates for incompatible fields as consequence of the Bourgain Brezis div-curl result[END_REF]. Lauteri and Luckhaus [START_REF] Lauteri | Geometric rigidity estimates for incompatible fields in dimension ≥ 3[END_REF] obtained the rigidity estimate (1.3) in the Lorentz space L 1 * ,∞ . In [START_REF] Lewintan | L p -versions of generalized Korn inequalities for incompatible tensor fields in arbitrary dimensions with p-integrable exterior derivative[END_REF] we have already established the corresponding results in the L p -setting. Here we focus on the trace-free case showing that the symmetric part can even be replaced by the symmetric deviatoric part.

Preliminaries and auxilliary results

By . ⊗ . we denote the dyadic product and by ., . the scalar product, so(n) := {A ∈ R n×n | A T = -A} is the Lie-Algebra of skew-symmetric matrices and Sym(n) := {X ∈ R n×n | X T = X}. Recall that usually the higher dimensional generalization of the curl is an operation curl n :

D (Ω, R n ) → D (Ω, R n(n-1) 2
) given by

curl n a :=           ∂ 1 a 2 -∂ 2 a 1 ∂ 1 a 3 -∂ 3 a 1 ∂ 2 a 3 -∂ 3 a 2 ∂ 1 a 4 -∂ 4 a 1 ∂ 2 a 4 -∂ 4 a 2 ∂ 3 a 4 -∂ 4 a 3 . . .           . (2.1)
Thus, in order to express this operation using the Hamiltonian formalism as a generalized cross product with ∇, we focus on the following bijection a n : so

(n) → R n(n-1) 2
given by

a n (A) := (A 12 , A 13 , A 23 , . . . , A 1n , . . . , A (n-1)n ) T for A ∈ so(n) (2.2a)
as well as its inverse A n : R

n(n-1) 2
→ so(n), so that

A n (a n (A)) = A for all A ∈ so(n) and a n (A n (a)) = a for all a ∈ R n(n-1) 2
(2.2b) and in coordinates it looks like

A n (a) =           0 a 1 a 2 a 4 . . .
-a 1 0 a 3 a 5 . . .

-a 2 -a 3 0 a 6 . . .

-a 4 -a 5 -a 6 0 . . .

• • • • • • • • • • • • 0           . (2.3)
Moreover, we will make use of the following notations

b = (b, b n ) T ∈ R n with b ∈ R n-1 (2.4a)
and

A =     0     -A e n T A e n A ∈ so(n) with A ∈ so(n -1).
(2.4b)

A generalized cross product

Regarding our goal to express curl n by the Hamiltonian formalism, we apply the following generalization of the cross product for n ≥ 2 acting as

× n : R n × R n → R n(n-1) 2 by a × n b := a n (a ⊗ b -b ⊗ a) for a, b ∈ R n . (2.5a)
Since for a fixed a ∈ R n the operation a × n . is linear in the second component there exists a unique matrix denoted by a ×n ∈ R n(n-1) 2

×n such that

a × n b =: a ×n b for all b ∈ R n .
(2.5b)

The matrices . ×n can be characterized inductively, and for a = (a, a n ) T the matrix a ×n has the form

a ×n =             a -a n • 1 n-1 0 a × n-1 where a 1 a 2 ×2 = -a 2 , a 1 , (2.6) so,
  a 1 a 2 a 3   ×3 =     -a 2 a 1 0 -a 3 0 a 1 0 -a 3 a 2     and     a 1 a 2 a 3 a 4     ×4 =           -a 2 a 1 0 0 -a 3 0 a 1 0 0 -a 3 a 2 0 -a 4 0 0 a 1 0 -a 4 0 a 2 0 0 -a 4 a 3           etc.
(2.7)

Remark 2.1. There are many possible identifications of skew-symmetric matrices with vectors. However, it is this matrix representation . ×n of the generalized cross product × n which allows us to establish the main identities needed for Lemma 2.9. Indeed, they were not easily visible to us before in index notations. Moreover, with this matrix representation in hand, the discussion of the boundary condition (see Observation 2.7) as well as the partial integration formula (2.56) are more transparent.

Remark 2.2. The entries of the generalized cross product a × 3 b, with a, b ∈ R 3 , are permutations (with a sign) of the entries of the classical cross product a × b. Recall, that also the operation a × . can be identified with a multiplication with the skew-symmetric matrix

Anti(a) =   0 -a 3 a 2 a 3 0 -a 1 -a 2 a 1 0   (2.8)
which differs from the expression a ×3 and from A 3 (a) which reads

A 3 (a) =   0 a 1 a 2 -a 1 0 a 3 -a 2 -a 3 0   .
(2.9)

Moreover, we have

a × n b = -b × n a = -b ×n a = a T -b ×n T T (2.10)
this allows us to define a generalized cross product of a vector b ∈ R n and a matrix P ∈ R n×m from the left and with a matrix B ∈ R m×n from the right via

b × n P := b ×n P ∈ R n(n-1) 2 ×m
to be seen as column-wise cross product, (

and

B × n b := B -b ×n T ∈ R m× n(n-1)
2 to be seen as row-wise cross product.

(

In such a way we obtain for all b ∈ R n :

1 n × n b = -b ×n T and b × n 1 n = b ×n . (2.12) Furthermore, for a, b ∈ R n it holds a ×n × n b = a ×n -b ×n T =             a -a n • 1 n-1 0 a × n-1         -b T 0 b n • 1 n-1 -b × n-1 T =     a ×n-1 -b ×n-1 T b n a ×n-1 a n b ×n-1 T -a ⊗ b -a n b n • 1 n-1     ∈ R n(n-1) 2 × n(n-1) 2 (2.13) end especially for a = b: b ×n × n b =     b ×n-1 -b ×n-1 T b n b ×n-1 b n b ×n-1 T -b ⊗ b -b n b n • 1 n-1     ∈ Sym n(n -1) 2 (2.14) Hence, for all a, b ∈ R n : tr( a ×n × n b) = a ×n , -b ×n = -(n -1) • a, b (2.15)
by induction over n, and, especially for a = b:

tr( b ×n × n b) = -(n -1) • b 2 . (2.16)
The entries of b ×n × n b are, by definition, linear combinations of b i b j , the entries of b ⊗ b. Interestingly, for n ≥ 3 also the converse holds true, i.e., the entries of b ⊗ b are linear combinations of the entries of b ×n × n b which will be assertion of the subsequent lemma. Moreover, we will use this as a key observation to achieve the existence of linear combinations in D

2 (A + ζ • 1 n ) = L(DCurl n (A + ζ • 1 n ))
for n ≥ 3, so that we can follow the argumentation scheme presented in n = 3 dimensions, cf. [START_REF] Lewintan | L p -trace-free generalized Korn inequalities for incompatible tensor fields in three space dimensions[END_REF], also in the higher dimensional case.

Lemma 2.3. For all b ∈ R n with n ≥ 3 we have:

b ⊗ b = L( b ×n × n b), (2.17)
denoting by L a corresponding constant coefficients linear operator.

Remark 2.4. There are no linear combinations (2.17) in n = 2 dimensions. Indeed, we only have

b 1 b 2 ×2 × 2 b 1 b 2 = -b 2 b 1 b 2 -b 1 = -b 1 2 -b 2 2 , (2.18) so that there are no linear expressions of b 1 2 , b 2 2 nor of b 1 b 2 in terms of the sole entry of b 1 b 2 ×2 × 2 b 1 b 2 .
Proof of Lemma 2.3 by induction over n ≥ 3. For the base case n = 3 we have

  b 1 b 2 b 3   ×3 × 3   b 1 b 2 b 3   (2.14) =     -b 1 2 -b 2 2 b 3 • -b 2 b 1 b 3 • -b 2 b 1 - b 1 b 2 ⊗ b 1 b 2 -b 3 2 • 1 2     (2.19) and tr     b 1 b 2 b 3   ×3 × 3   b 1 b 2 b 3     = -2 • (b 1 2 + b 2 2 + b 3 2 ). (2.20) Thus, for all b ∈ R 3 b 3 2 = e 1 , ( b ×3 × 3 b) e 1 - 1 2 tr( b ×3 × 3 b) = L( b ×3 × 3 b) (2.21)
and consequently from the expression (2.19) we conclude

b ⊗ b = L( b ×3 × 3 b) ∀ b ∈ R 3 . (2.22)
Now, assume for the inductive step that for all b ∈ R n-1 with n ≥ 4 we have

b ⊗ b = L( b ×n-1 × n-1 b). (2.23) For b ∈ R n we have b ⊗ b =   b ⊗ b b n • b b n • b T b n 2   . (2.24) Surely, b n • b = L(b n • b ×n-1 ) = L( b ×n × n b) by the expression (2.

14). The induction hypothesis gives

b ⊗ b = L( b ×n-1 × n-1 b) (2.14) = L( b ×n × n b), (2.25) hence, also b 2 = L( b ×n × n b), so that finally b 2 n (2.16) = - 1 n -1 tr( b ×n × n b) -b 2 = L( b ×n × n b).
By definition (2.11b), the entries of B × n b are linear combinations of the entries B ij b k , i.e., of the entries of the matrix B multiplied with the entries of the vector b. For skew-symmetric matrices also the converse holds true. This is the assertion of the next lemma.

Lemma 2.5. For all A ∈ so(n) and b ∈ R n with n ≥ 2 we have

(i) A × n b = L(a n (A) ⊗ b) (ii) a n (A) ⊗ b = L(A × n b)
denoting by L a corresponding constant coefficients linear operator which can differ in both cases.

Proof. For A ∈ so(n) and b ∈ R n we have

a n (A) ⊗ b = a n-1 (A) A e n ⊗ b b n =   a n-1 (A) ⊗ b b n • a n-1 (A) A e n ⊗ b b n • A e n   (2.26)
and on the other hand

A × n b =     0     -A e n T A e n A         -b T 0 b n • 1 n-1 -b × n-1 T (2.27) =     A -b ×n-1 T b n • A -A e n ⊗ b A e n T b ×n-1 T -b n • A e n T     =   A × n-1 b b n • A -A e n ⊗ b -A e n × n-1 b T -b n • A e n T   .
Thus, the conclusions of the lemma follow by induction over the dimension n. Indeed, for the base case n = 2 we have 0 α

-α 0 × 2 b 1 b 2 = 0 α -α 0 -b 2 b 1 = -α b 1 -α b 2 = -(α ⊗ b) T (2.28)
which establishes (i) and (ii) of the Lemma for n = 2.

For the inductive step let us assume that the statement of the Lemma holds for all A ∈ so(n -1) and all b ∈ R n-1 , i.e., it holds:

(i) A × n-1 b = L(a n-1 (A) ⊗ b) and (ii) a n-1 (A) ⊗ b = L(A × n-1 b).
Thus, returning to A ∈ so(n) and b ∈ R n we have by the expressions (2.26) and (2.27), respectively,

(i) A × n-1 b = L(a n (A) ⊗ b) and (ii) a n-1 (A) ⊗ b = L(A × n b)
and the conclusion of part (i) of the Lemma follows then from the definition of the generalized cross product, indeed,

-A e n × n-1 b T = a n-1 (b ⊗ A e n -A e n ⊗ b) T (2.26) = L(a n (A) ⊗ b). (2.29)
On the other hand we have This finishes the proof of (ii) since we have shown that all the entries of a n (A) ⊗ b can be written as linear combinations of the entries of A × n b.

b n • A -A e n ⊗ b -b n • A -A e n ⊗ b T + A n-1 A e n × n-1 b = 2 b n • A, ( 2
Remark 2.6. The identity (2.30) is not a new result, and usually it is expressed using coordinates:

(A × n b) kij -(A × n b) kji + (A × n b) jik = 2 A ij b k ∀ i, j, k = 1, . . . , n.
(2.31) However, we included the statement as well as the proof not only for the sake of completeness, but also since the use of the matrix representation of the generalized cross product allows us to give a coordinate-free proof and provides a deeper insight in the algebra needed in the present paper.

For a square matrix P ∈ R n×n we can take the generalized cross product with a vector b ∈ R n both left and right, and simultaneously we obtain for the identity matrix

b × n 1 n × n b = b ×n 1 n ( -b ×n ) T = b ×n ( -b ×n ) T = b ×n × n b (2.14) ∈ Sym n(n -1) 2 (2.32)
and for a general matrix P ∈ R n×n

(b × n P × n b) T = b ×n P ( -b ×n ) T T = -b ×n P T ( b ×n ) T = b × n P T × n b. (2.33)
Thus, especially for a symmetric matrix S ∈ Sym(n) and a skew-symmetric matrix A ∈ so(n) we obtain

b × n S × n b ∈ Sym n(n -1) 2 and b × n A × n b ∈ so n(n -1) 2 .
(2.34)

Observation 2.7. Let A ∈ so(n) and α ∈ R. Then (A + α • 1 n ) × n b = 0 for b ∈ R n \{0} implies A = 0 and α = 0.
Proof. Taking the generalized cross product from the left on both sides of 0

= A × n b + α • 1 n × n b gives 0 = b × n A × n b skew-symmetric + α • b × n 1 n × n b symmetric , (2.35)
so that taking the trace of the symmetric part on both sides we obtain

0 = α • tr(b × n 1 n × n b) (2.32) = α • tr( b ×n × n b) (2.16) = -α • (n -1) • b 2 (2.36)
which implies α = 0. Consequently we moreover have A × n b = 0 which by Lemma 2.5 (ii) yields a n (A) ⊗ b = 0, and thus a n (A) = 0 and A = 0.

Considerations from vector calculus

Subsequently we make use of the algebraic behavior of the vector differential operator ∇ as a vector for formal calculations. So, the derivative and the divergence of a vector field a ∈ D (Ω, R n ) can be seen as Da = a ⊗ ∇ and div a = a, ∇ = Da, 1 n = tr(Da).

(2.37)

In a similar way the generalized curl is related to the generalized cross product

curl n a := a × n (-∇) = ∇ × n a = ∇ ×n a = a T ( ∇ ×n ) T T . ( 2

.38)

The latter expression gives a generalized row-wise matrix Curl n operator for B ∈ D (Ω, R m×n ) via

Curl n B := B × n (-∇) = B ∇ ×n T . (2.39)
This differential operator kills the Jacobian matrix of a vector field (a compatible field), indeed

Curl n Da = Da ∇ ×n T = (a ⊗ ∇) ∇ ×n T = a ∇ ×n ∇ T = a (∇ × n ∇) T ≡ 0, (2.40) since b × n b = 0 for all b ∈ R n . Furthermore, for a scalar field ζ ∈ D (Ω, R) we obtain Curl n (ζ • 1 n ) = 1 n ∇ζ ×n T = ∇ζ ×n T .
(2.41)

For P ∈ D (Ω, R n×m ) we consider also the column-wise differential operator of first order coming from a cross product form the left, namely ∇ × n P = ∇ ×n P, (2.42) which kills the transposed Jacobian (Da) T . It is clear, that Curl n B = L(DB), i.e., that the entries of Curl n B are linear combinations of the entries of DB for all B ∈ D (Ω, R m×n ). However, for skew-symmetric matrix fields also the converse holds true:

Corollary 2.8. For all A ∈ D (Ω, so(n)) with n ≥ 2 it holds: DA = L(Curl n A).
It is a well known result and follows from the linear expression (ii) in Lemma 2.5 replacing b by -∇ as well as from its analogous statement written out in coordinates (2.31).

We now catch up with the crucial linear relation needed in our argumentation scheme.

Lemma 2.9.

Let n ≥ 3, A ∈ D (Ω, so(n)) and ζ ∈ D (Ω, R). Then the entries of D 2 (A + ζ • 1 n ) are linear combinations of the entries of DCurl n (A + ζ • 1 n ).
Remark 2.10. The statement is false in n = 2 dimensions. Indeed, with α, ζ ∈ D (Ω, R) we have

Curl 2 ζ α -α ζ = ζ α -α ζ -∂ 2 ∂ 1 = ∂ 1 α -∂ 2 ζ ∂ 2 α + ∂ 1 ζ (2.43) so that DCurl 2 ζ α -α ζ = ∂ 1 ∂ 1 α -∂ 1 ∂ 2 ζ ∂ 1 ∂ 2 α -∂ 2 ∂ 2 ζ ∂ 1 ∂ 2 α + ∂ 1 ∂ 1 ζ ∂ 2 ∂ 2 α + ∂ 1 ∂ 2 ζ (2.44)
and we cannot extract ∂ 1 ∂ 1 α from the components of (2.44).

Proof of Lemma 2.9. The proof is divided into the two observations

1. D 2 ζ = L(DCurl n (A + ζ • 1 n )), 2. D 2 A = L(DCurl n (A + ζ • 1 n ))
denoting by L a corresponding constant coefficients linear operator which can differ in both cases. To show that the entries of the Hessian D 2 ζ can be written as linear combinations of the entries of DCurl n (A + ζ • 1 n ) we introduce the following second order derivative operator for square matrix fields P ∈ D (Ω, R n×n ):

inc n P := ∇ × n P × n ∇ = -∇ × n (Curl n P ) (2.45)
in the style of the incompatibility operator known from the three-dimensional case. In regard of (2.32) we see

inc n (ζ • 1 n ) = ∇ × n (ζ • 1 n ) × n ∇ = ∇ ×n × n ∇ ζ ∈ Sym n(n -1) 2 , (2.46)
so that substituting b by -∇ in Lemma 2.3 we obtain

D 2 ζ = (∇ ⊗ ∇) ζ = L ∇ ×n × n ∇ ζ = L(inc n (ζ • 1 n )). (2.47)
Moreover, with regard to (2.34) 2 we have for a skew-symmetric matrix field A ∈ D (Ω, so(n)):

inc n A = ∇ × n A × n ∇ ∈ so n(n -1) 2 , (2.48)
concluding for the 1. part that

D 2 ζ = L(inc n (ζ • 1 n )) = L(sym(inc n (A + ζ • 1 n ))) = L(inc n (A + ζ • 1 n )) = L(DCurl n (A + ζ • 1 n )) (2.49)
where in the last step we have used that the entries of inc n P = -∇ × n (Curl n P ) are, of course, linear combinations of the entries from DCurl n P .

To establishes part 2 recall that the entries of DA for a skew-symmetric matrix field are linear combinations of the entries of Curl n A, giving

DA = L(Curl n A) = L(Curl n (A + ζ • 1 n )) -L(Curl n (ζ • 1 n )) (2.41) = L(Curl n (A + ζ • 1 n )) -L ∇ζ ×n T . (2.50)
The conclusion follows by taking the ∂ j -derivative of (2.50) together with the observation from the 1. part:

∂ j DA = L(∂ j Curl n (A + ζ • 1 n )) -L ∂ j ∇ζ ×n T (2.49) = L(DCurl n (A + ζ • 1 n )).
In the last result of this section we focus on the kernel of dev n sym and Curl n :

Lemma 2.11. Let n ≥ 3, A ∈ L p (Ω, so(n)) and ζ ∈ L p (Ω, R). Then Curl n (A + ζ • 1 n ) ≡ 0 in the distributional sense if and only if there exist constant b ∈ R n , d ∈ R n(n-1) 2 , β ∈ R such that A + ζ • 1 n = A n -b ×n x + d + ( b, x + β) • 1 n almost everywhere in Ω.
Remark 2.12. The "only if"-part is false in n = 2 dimensions. To see this, take in (2.43) α and ζ to be the real an imaginary part of a holomorphic function.

Proof of Lemma 2.11. For the "if"-part we have

Curl n (A n -b ×n x + b, x • 1 n ) = Curl n (A n -b ×n x ) + ∇ b, x ×n
b ×n T = Curl n (ζ • 1 n ) = -Curl n A = -B T -C T (2.53)
implying that C ≡ 0 and B = b ×n almost everywhere in Ω.

Remark 2.13. The expression of the kernel follows also from the consideration for the classical trace-free Korn inequalities. Indeed, on simply connected domains, Curl n P ≡ 0 implies that P = Du for a vector field u ∈ W 1, p (Ω, R n ). Thus, the condition dev n sym P ≡ 0 reads dev n sym Du ≡ 0, whose solutions are well known as infinitesimal conformal mappings, cf. [START_REF] Bauer | Dev-Div-and DevSym-DevCurl-inequalities for incompatible square tensor fields with mixed boundary conditions[END_REF][START_REF] Dain | Generalized Korn's inequality and conformal Killing vectors[END_REF][START_REF] Jeong | Existence, uniqueness and stability in linear Cosserat elasticity for weakest curvature conditions[END_REF][START_REF] Neff | Subgrid interaction and micro-randomness -novel invariance requirements in infinitesimal gradient elasticity[END_REF][START_REF] Reshetnyak | Estimates for certain differential operators with finite-dimensional kernel[END_REF][START_REF] Reshetnyak | Stability Theorems in Geometry and Analysis[END_REF][START_REF] Schirra | New Korn-type inequalities and regularity of solutions to linear elliptic systems and anisotropic variational problems involving the trace-free part of the symmetric gradient[END_REF].

Function spaces

Having above relations at hand we can now catch up the arguments from [START_REF] Lewintan | L p -trace-free generalized Korn inequalities for incompatible tensor fields in three space dimensions[END_REF]. Let Ω ⊆ R n , we start by defining the space

W 1, p (Curl n ; Ω, R n×n ) := {P ∈ L p (Ω, R n×n ) | Curl n P ∈ L p (Ω, R n× n(n-1)

2

)} (2.54a) equipped with the norm

P W 1, p (Curln;Ω,R n×n ) := P p L p (Ω,R n×n ) + Curl n P p L p (Ω,R n× n(n-1)
2

)
1 p

(2.54b) and its subspace W 1, p 0 (Curl n ; Ω, R n×n ) as the completion of C ∞ 0 (Ω, R n×n ) in the W 1, p (Curl n ; Ω, R n×n )norm.

In our proofs we shall use an important equivalence of norms due to Nečas [13, Théorème 1] valid on bounded Lipschitz domains, cf. also discussion in [START_REF] Lewintan | L p -trace-free generalized Korn inequalities for incompatible tensor fields in three space dimensions[END_REF][START_REF] Lewintan | Nečas-Lions lemma revisited: An L p -version of the generalized Korn inequality for incompatible tensor fields[END_REF] and the references cited therein.

Thus, in what follows Ω ⊂ R n will be a bounded domain with Lipschitz boundary and we are allowed to define boundary conditions in the distributional sense, so that

W 1, p 0 (Curl n ; Ω, R n×n ) = {P ∈ W 1, p (Curl n ; Ω, R n×n ) | P × n ν = 0 on ∂Ω} = {P ∈ W 1, p (Curl n ; Ω, R n×n ) | P τ l = 0 on ∂Ω for all l = 1, . . . , n -1}, (2.55)
where ν stands for the outward unit normal vector field and {τ l } l=1,...,n-1 denotes a moving tangent frame on ∂Ω, cf. [START_REF] Lewintan | L p -versions of generalized Korn inequalities for incompatible tensor fields in arbitrary dimensions with p-integrable exterior derivative[END_REF]. Here, the generalized tangential trace P × n ν is understood in the sense of W -1 p , p (∂Ω, R n× n(n-1)

2

) which is justified by partial integration, so that its trace is defined by

∀ Q ∈ W 1-1 p , p (∂Ω, R n× n(n-1)

2

) :

(2.56)

P × n (-ν), Q ∂Ω = Ω Curl n P, Q R n× n(n-1) 2 + P, Q ∇ ×n R n×n dx having denoted by Q ∈ W 1, p (Ω, R n× n(n-1) 2
) any extension of Q in Ω, where, ., . ∂Ω indicates the duality pairing between W -1 p , p (∂Ω, R n× n(n-1)

2

) and W

1-1 p , p (∂Ω, R n× n(n-1) 2
). Indeed, for smooth P and Q on Ω we have

∂Ω P × n (-ν), Q R n× n(n-1) 2 dS (2.11b) = ∂Ω P ( ν ×n ) T , Q R n× n(n-1) 2 dS = ∂Ω ( ν ×n ) T , P T Q R n× n(n-1) 2 dS ( * ) = Ω ( ∇ ×n ) T , P T Q R n× n(n-1) 2 dx (2.57) = Ω P ( ∇ ×n ) T , Q R n× n(n-1) 2 + P, Q ∇ ×n R n×n dx = Ω Curl n P, Q R n× n(n-1) 2 + P, Q ∇ ×n R n×n dx,
where in ( * ) we have used the fact that we only deal with linear combinations of partial derivatives and from the classical divergence theorem it holds

∂Ω ν i ζ dS = Ω ∂ i ζ dx (2.58)
for a smooth scalar function ζ on Ω. Further, following [START_REF] Lewintan | Nečas-Lions lemma revisited: An L p -version of the generalized Korn inequality for incompatible tensor fields[END_REF] we introduce also the space W 1, p Γ,0 (Curl n ; Ω, R n×n ) of functions with vanishing tangential trace only on a relatively open (non-empty) subset Γ ⊆ ∂Ω of the boundary by completion of C ∞ Γ,0 (Ω, R n×n ) with respect to the W 1, p (Curl n ; Ω, R n×n )-norm.

Remark 3.4. For compatible displacement gradients P = Du we get back from (3.4) the quantitative version of the classical trace-free Korn's inequality, cf. [START_REF] Dain | Generalized Korn's inequality and conformal Killing vectors[END_REF][START_REF] Reshetnyak | Stability Theorems in Geometry and Analysis[END_REF][START_REF] Schirra | New Korn-type inequalities and regularity of solutions to linear elliptic systems and anisotropic variational problems involving the trace-free part of the symmetric gradient[END_REF].

Finally, we examine the effect of tangential boundary conditions P × n ν ≡ 0.

Theorem 3.5. Let n ≥ 3, Ω ⊂ R n be a bounded Lipschitz domain and 1 < p < ∞. There exists a constant c = c(n, p, Ω) > 0, such that for all P ∈ W 1, p 0 (Curl n ; Ω, R n×n ) we have P L p (Ω,R n×n ) ≤ c dev n sym P L p (Ω,R n×n ) + Curl n P L p (Ω,R n× n(n-1)

2

)

.

(3.6)

Proof. We follow the same argumentation scheme as in the proof of [START_REF] Lewintan | Nečas-Lions lemma revisited: An L p -version of the generalized Korn inequality for incompatible tensor fields[END_REF]Theorem 3.5] and consider a sequence {P k } k∈N ⊂ W 1, p 0 (Curl n ; Ω, R n×n ) converging weakly in L p (Ω, R n×n ) to some P * so that dev n sym P * = 0 almost everywhere and Curl n P * = 0 in the distributional sense, i.e. P * ∈ K dS,Cn . By (2.56) we obtain that P * × n (-ν), Q ∂Ω = 0 for all Q ∈ W 1, p (Ω, R n× n(n-1)

2

). However, the boundary condition P * × n ν = 0 is also valid in the classical sense, since P * ∈ K dS,Cn has an explicit representation. Using the explicit representation of P * = A n -b ×n x + d + ( b, x + β) • 1 n , we conclude using Observation 2.7 that, in fact, P * ≡ 0: Remark 3.6. Estimate (3.6) should persist also in n = 2 dimensions. So, the case p = 2 is already contained in [START_REF] Bauer | Dev-Div-and DevSym-DevCurl-inequalities for incompatible square tensor fields with mixed boundary conditions[END_REF]. However, for the general case p ∈ (1, ∞) we need a different approach and it will be the subject of a forthcoming note.

  .30) so that, b n •A = L(A× n b) and also b n •a n-1 (A) = L(A× n b) which by (2.27) implies that A e n ⊗b = L(A× n b).

T= 2 ,

 2 -b ×n T + b ×n T = 0. Conversely, Curl n (A + ζ • 1 n ) ≡ 0 in the distributional sense implies by (2.49): D 2 ζ ≡ 0 and D 2 A ≡ 0 in the distributional sense, (2.51) thus, ζ(x) = b, x + β and A = A n (B x + d) a.e. (2.52) for some b ∈ R n , β ∈ R, B ∈ R n(n-1) 2 ×n and d ∈ R n(n-1) and we have Curl n (ζ • 1 n ) = ∇ζ ×n T = b ×n T as well asCurl n A = Curl n (A n (B x)) = B T + C T , where C ∈ R n(n-1)2×n has only possibly non-zero entries at those positions at which the matrix b ×n has zeros. Hence, the condition, Curl n (A + ζ • 1 n ) ≡ 0 gives:

[

  A n -b ×n x + d + ( b, x + β) • 1 n ] × n ν = 0 Obs. 2.7 ⇒ -b ×n x + d = 0 and b, x + β = 0 for all x ∈ ∂Ω ⇒ β = 0, b = 0, d = 0 .

Remark 3 . 7 .Remark 3 . 8 .

 3738 For compatible P = Du we recover from (3.6) a tangential trace-free Korn inequality. For n ≥ 3, the previous results also hold true for tensor fields with vanishing tangential trace only on a relatively open (non-empty) subset Γ ⊆ ∂Ω of the boundary, cf. discussion in[START_REF] Lewintan | Nečas-Lions lemma revisited: An L p -version of the generalized Korn inequality for incompatible tensor fields[END_REF]. But, this is not the case in n = 2 dimensions. Indeed, already the trace-free version of Korn's first inequality (1.1) with only partial boundary condition is false in the n = 2 case, cf. e.g. the counterexample contained in [1, section 6.6].

Chair for Nonlinear Analysis and Modeling, Faculty of Mathematics, University of Duisburg-Essen, Thea-Leymann-Str. 9, 45127 Essen, Germany * Corresponding author, email: peter.lewintan@uni-due.de.

Acknowledgment: The authors thank the referees for their valuable comments. This work was initiated in the framework of the Priority Programme SPP 2256 'Variational Methods for Predicting Complex Phenomena in Engineering Structures and Materials' funded by the Deutsche Forschungsgemeinschaft (DFG, German research foundation), Project-ID 422730790. The second author was supported within the project 'A variational scale-dependent transition scheme -from Cauchy elasticity to the relaxed micromorphic continuum' (Project-ID 440935806). Moreover, both authors were supported in the Project-ID 415894848 by the Deutsche Forschungsgemeinschaft.

Trace-free Korn inequalities for incompatible tensors in higher dimensions

With the auxilliary results in hand we can now catch up with the higher dimensional versions of the results presented in [START_REF] Lewintan | L p -trace-free generalized Korn inequalities for incompatible tensor fields in three space dimensions[END_REF].

Lemma 3.1. Let n ≥ 3, Ω ⊂ R n be a bounded Lipschitz domain and 1 < p < ∞. Then P ∈ D (Ω, R n×n ), dev n sym P ∈ L p (Ω, R n×n ) and Curl n P ∈ W -1, p (Ω, R n× n(n-1)

2

) imply P ∈ L p (Ω, R n×n ). Moreover, we have the estimate

with a constant c = c(n, p, Ω) > 0.

Proof. We have to show that skew P + 1 n tr P • 1 n ∈ L p (Ω, R n×n ) follows from the assumptions of the lemma. By the linearity of differential operator DCurl n and the orthogonal decomposition

Thus, by the assumed regularity of the right hand side, it follows that the left hand side belongs to

×n ). Furthermore, we have

By Lemma 2.9 we obtain D 2 (skew P + 1 n tr P • 1 n ) ∈ W -2, p and an application of the Lions lemma resp. Nečas estimate [9, Thm. 2.7 and Cor. 2.8] to skew P + 1 n tr P • 1 n yield the conclusions. Eliminating the first term on the right hand side of (3.1) gives: Theorem 3.2. Let n ≥ 3, Ω ⊂ R n be a bounded Lipschitz domain and 1 < p < ∞. There exists a constant c = c(n, p, Ω) > 0, such that for all P ∈ L p (Ω, R n×n ) we have inf

where the kernel is given by

, β ∈ R} .


Charlène Lambert
email: charlene.lambert@agroparistech.fr

Julien Cartailler
email: julien.cartailler@agroparistech.fr

Sandy Rouchouse
email: s.rouchouse@tecaliman.com

Giana Almeida Perré
email: giana.perre@agroparistech.fr

Francis F Courtois
email: francis.courtois@agroparistech.fr

Giana Almeida

Characterization and modeling of cooling and drying of pellets for animal feed

Keywords: cooling, drying, pellets, animal feed, characterization, model

Introduction

The process by which pellets are manufactured consists of several steps: grinding, weighing raw materials, mixing, pelleting to agglomerate small particles into larger ones, cooling, and conditioning. Despite its economic importance, this process has been largely ignored in the literature, with only a few exceptions. The limited research conducted to date has mainly focused on the effect of pellet formulation on the manufacturing process [START_REF] Wood | The functional properties of feed raw materials and their effect on the production and quality of feed pellets[END_REF][START_REF] Thomas | Physical quality of pelleted animal feed 1: Criteria for pellet quality[END_REF][START_REF] Thomas | Physical quality of pelleted animal feed 3. contribution of feedstuff components[END_REF][START_REF] Hemmingsen | Water adsorption in feed ingredients for animal pellets at different temperatures, particle size, and ingredient combinations[END_REF], or the conditioning [START_REF] Thomas | Physical quality of pelleted animal feed. 2. contribution of processes and its conditions[END_REF] or pelleting [START_REF] Thomas | Physical quality of pelleted animal feed. 2. contribution of processes and its conditions[END_REF][START_REF] De Blank | On-line control and optimisation of the pelleting process of animal feed[END_REF] steps of the process.

Even if cooling is not one of the most energy-consuming steps of this process, industrial investment has recently focused on the use of adjustable-speed drives to reduce fan energy consumption. This technology conserves energy by adapting air flow to the actual needs of the pellets at any given time. However, for most of these coolers, there is no real-time feedback between measurements of pellet characteristics (e.g., humidity, temperature) and processing parameters. Ideally, measurements of air humidity taken from a cooler's exhaust could be used to optimize the drying-bed height or fan velocity. Moreover, the lack of automatic control of the cooling process may lead to overdrying of pellets, resulting in higher energy consumption and cracks in the pellet surface which affect their durability [START_REF] Thomas | Physical quality of pelleted animal feed. 2. contribution of processes and its conditions[END_REF]. Hence, precise control of this step of the process is of the utmost importance for the animal-feed industry in order to optimize energy consumption and ensure product quality. For this reason, the French Environment and Energy Management Agency (ADEME) funded the OPSERA project ("Optimizing Dryer-Cooler Performance for Animal Feed"), which aimed to optimize the drying-cooling process in real time using an advanced control algorithm. The control scheme centers on the estimation of pellet moisture content by a smart sensor, which utilizes observation/estimation techniques based on a dynamic simulation model of the drying-cooling process. This paper describes the characterization and modeling of the drying and cooling of pellets as required for the real-time estimation of their moisture content.

In some respects, the lack of scientific attention given to the cooling process (with the exception of [START_REF] Kelley | Evaluation of chilled cooling of feed pellets[END_REF]Maier & Bakker-Arkema, 1992;Maier et al., 1992)) is surprising, given that the cooling operation is also a drying one: the cold air blown is used to decrease both pellet moisture content and temperature. Drying/cooling thus represents a crucial final step in pellet production. However, the strong coupling between cooling and drying phenomena makes the operation difficult to optimize in the absence of accurate simulations. Over two decades ago, a drying-cooling model was developed by (Maier & Bakker-Arkema, 1992) and its simulations were validated at a pilot scale. However, this model was only able to predict the trend of moisture and temperature effects on cooler performance [START_REF] Kelley | Evaluation of chilled cooling of feed pellets[END_REF]Maier & Bakker-Arkema, 1992;Maier et al., 1992) and did not accurately predict moisture loss and temperature in pellets (up to 60% deviation was observed between experimental and simulated data). Hence, what is needed today is a more accurate model that can control an industrial cooler and predict the effects of changes in pellet formulation on operating conditions. Ideally, such a model would also be able to take into account changes in pellet composition, as pellet recipes are changed continuously in response to variations in the cost of ingredients (soya, corn, etc.). These changes are likely to affect drying-cooling performance and represent a challenge to model development.

The purpose of this work was to develop a convective drying-cooling model that could be implemented in two versions: 1) an unsteady-state version, in order to control process parameters, and 2) a steady-state version, in computer-aided engineering (CAE) software. In a recent work discussing classical modeling approaches of convective drying processes (Lambert et al., 2015a), it was argued that the main characteristics of a successful model should be a low computational requirement and a wide range of validity. Simulation of the drying or drying-cooling operation should not last longer than a few seconds, and the model should also have multi-product compatibility and extrapolation capability. It should be able to simulate all classical drying and drying-cooling conditions. In addition, it should allow the late addition of new products or at least new formulations without the need to modify the equations. Thus, parameters of the model should have a physical meaning and should be calculated with correlations, experimentally measured, or identified by reverse methods (Lambert et al., 2015a).

The aim of our model was to predict the final moisture content of pellets and the energy consumption of the drying-cooling process. However, obtaining reliable drying-cooling kinetics is difficult, especially in the thin-layer configuration. Hence, we decided to first focus on a drying-only model. We then attempted to validate this model using deep-bed cooling-drying kinetics.

Materials and methods

There is almost an infinite number of formulations of animal feed pellets. This study focused on the experimental characterization of two formulations, with various recipes and pelleting conditions. Water activity and effective moisture diffusivity were identified using DVS data and drying kinetics, respectively.

Origin of the pellets

The two formulations of animal feed pellets studied were pellets for chickens (P4C) and for rabbits (P4R). In the latter case, we studied two different geometries of the same composition (same raw materials and recipe). These were designated P4R -2.5 mm and P4R -4 mm. All formulations had a cylindrical shape. Raw seed composition and dry matter composition of both pellet formulations are given in tables 1 and 2, respectively. For the pelleting process, ingredients were mixed with steam in a conditioner then pelleted in a press. The temperatures of the pellets at the press outlet were 84°C (P4C), 87°C (P4R -2.5 mm), and 83°C (P4R -4 mm). The pellet moisture content at the press outlet was 0. For each formulation and geometry, 2.5-kg samples of pellets were vacuum-packed in plastic boxes and later stored in a cold room.

Moisture content measurement

Pellet moisture content was measured according to the French AFNOR norm NF V 18-109, before thinlayer and deep-bed drying and after thin-layer drying. The pellets were first ground using a water-cooled miller (IKA WERKE KM 20). The measurement method assumes a weight loss of 5.00 g (± 0.05 g) of pellet powder during a 4-hour stay in a 103°C oven. Samples were then placed in a desiccator for at least 1 hour before being weighed. The moisture content was averaged over three samples.

Geometry measurement

Apparent density of pellet dry matter was measured using Fontainebleau sand. Real density of pellet dry matter was determined using an air pycnometer (ACCUPYC 1330). Pellet dimensions were measured using two methods: a manual approach using a caliper and an automated approach based on image analysis (using a flatbed scanner and ImageJ software). Porosity of the bed was determined according to equation 1.

ε = ρ real dm -ρ ap dm ρ r dm ( 1
)
where ε is the bed porosity (decimal), ρ real dm the real density of the dry matter (kg.m -3 ) and ρ ap dm the apparent density of the dry matter (kg.m -3 ).

Sorption isotherm measurement

Isothermal analyses were performed using a DVS Intrinsic apparatus (Surface Measurement Systems, London, UK). An ultra-sensitive micro-balance enabled the measurement of mass variations as small as 1 part per 10 million. First, samples were ground using a water-cooled miller (IKA WERKE KM 20).

Approximately 10 mg of each sample was inserted in the pan. In the DVS device, a constant flow of nitrogen gas was mixed with a preset amount of water vapor. Nitrogen then passed through the chamber to maintain the desired relative humidity level, ± 0.5%. All sorption cycles started with 0% relative humidity (RH). were based on actual measured RH values (which occasionally differed from target values). The instrument maintained a constant RH until the change in sample moisture content was less than 0.002% per minute over a 10-min period. This threshold value was obtained from [START_REF] Hill | The water vapor sorption behavior of natural fibers[END_REF], in order to obtain a final sample moisture content within 0.1% of the equilibrium value over an extended period of time.

Equilibrium moisture content (X eq ) and relative humidity (RH) / water activity (a w ) data were fitted using a modified Henderson equation (equation 2) [START_REF] Henderson | A basic concept of equilibrium moisture[END_REF][START_REF] Thompson | Mathematical simulation of corn drying: A new model[END_REF] without temperature dependency and the Chung-Pfost correlation (equation 3) [START_REF] Chung | Adsorption and desorption of water vapor by cereal grains and their products. part i: Heat and free energy changes of adsorption and desorption[END_REF]. The Henderson equation is one of the five equations recommended by ASAE Standard D245.6 for fitting sorption data from several agri-food materials [START_REF] Asae | D245.6. moisture relationships of plant-based agricultural products[END_REF]. The Chung-Pfost correlation was used by (Maier & Bakker-Arkema, 1992) and [START_REF] Brook | Handbook of transportation and marketing in agriculture. chapter Drying, cleaning and condtionning[END_REF] but, as will be discussed later, this equation (3) cannot guarantee that when X eq = 0 then a w = 0, which is a major flaw for simulators.

a w = 1 -exp (-p 1 (X eq ) p2 )
(2)

X eq = p 1 -p 2 ln (-p 3 (Θ + p 4 ) ln (a w )) (3)
where p 1 -p 4 are constants that need to be identified using experimental data, X eq is the equilibrium moisture content (% d.b.), and Θ is the product surface temperature (°C).

Drying experiment

Pellets were dried using two different protocols: a conventional hot-air drying approach (i.e. several thin layer kinetics and one deep bed kinetic) and a drying-cooling method (deep-bed kinetics only). Measurements of thin-layer kinetics were performed in triplicate. For both methods, experimental data were recorded using an automated pilot dryer located in an AgroParisTech/INRA drying laboratory (Massy, France) under constant conditions (see below).

Description of the drying pilot

A schematic of the drying pilot is given in figure 1. In this pilot, the air passes, from the top to the bottom, through a cylindrical drying chamber (25 cm in diameter and 28 cm in height). The following variables were recorded every minute (with an OPTO22 SNAP-PAC-EB2 data logger): For measurements of relative humidity and temperature in the pilot dryer, the following uncertainties (using ± 3 σ) were determined:

Inlet air temperature: ± 1.5°C (σ = 0.5°C);

Outlet air temperature: ± 2.5°C (σ = 0.8°C);

Inlet and outlet relative humidities: ± 5% (σ = 1.7%)

However, actual deviations in measurements of outlet air temperature and relative humidity may have been greater than those estimated here, for three reasons. First, due to the intermittent nature of the mass measurement system, minor air leaks and locally reduced thermal insulation below the sieve were unavoidable. Second, it is possible that having the sensor close to the sieve introduced bias due to thermal radiation. Finally, due to the depth of the beds, air flow below the sieve was less homogeneous. Because of this, the outlet air temperature may have been slightly underestimated. In addition, the response time of the relative humidity sensor (45 s) was not negligible compared to the sudden variations in relative humidity in the first minutes of the drying process. As a consequence, a comparison of experimental and simulated data during this initial period was undertaken very carefully (see discussion).

Drying and drying-cooling protocols

Drying method. Thin-layer and deep-bed drying kinetics were recorded using P4Cs and P4Rs (P4R -4 mm and P4R -2.5 mm). About 200 g (thin layers) and between 2.5 to 5 kg (deep beds approximately 0.10-to 0.20-m high) of pellets were introduced into the dryer. Operating conditions of the dryer are summarized in table 3.

Drying-cooling method. As stated in the introduction, it is quite difficult to precisely control the initial conditions of the pellets entering the drying-cooling experiment. It takes considerable effort to warm up the pellets without drying them and then transfer them to the dryer without the premature loss of water or heat.

Even when this can be done, it is extremely challenging to then obtain accurate measurements of the actual initial product temperature and moisture content at the beginning of the drying-cooling step. In this work, no reliable data could be obtained for thin-layer drying-cooling kinetics, hence, the only drying-cooling data available were for the deep-bed drying kinetics of P4Cs. About 5 kg of pellets were vacuum-packed using a vacuum heat sealer (MULTIVALUED GASTRONOMIC) and pouches composed of low-density polyethylene (LEXINGTON). Packing pellets in a vacuum-sealed bag limited pellets' contact with air and minimized mass transfers during preheating. The pellets and the sieve were then pre-heated as follows:

The sieve of the dryer was preheated to about 100°C using a heat gun (METABO H 1600) in order to be at nearly 80°C at the beginning of the experiment;

Pellets, in their vacuum-sealed bags, were pre-heated in a 100°C oven for one hour. The core temperature of the pellets was about 60 -70°C when taken out of the oven.

Pellets were then introduced into the drying zone after stabilization of the dryer's operating conditions (tested conditions are listed in table 3). A sub-sample was withdrawn for moisture content measurement.

The initial temperature of the pellets was initially estimated through the use of an infrared thermometer, but these measurements were not used because of suspected bias. A better estimate was obtained from the recorded outlet air temperature after two minutes of drying-cooling.

Identification of D ef f

This study continues previous work that discussed three different approaches for the identification of effective moisture diffusivity D ef f . [START_REF] Lambert | Reverse methodology to identify moisture diffusivity during air-drying of foodstuffs[END_REF]. Here, these three strategies were used for the identification of the effective moisture diffusivity of pellets as summarized in table 4.

As a reminder, the classical strategy used in food drying consists of four steps:

1. Measure in the laboratory the accessible properties of the product (water activity, heat capacity, geometry, and porosity);

2. Collect experimental measurements of thin-layer drying kinetics under constant air drying conditions over a wide range of air temperatures;

3. Identify the effective moisture diffusivity coefficient by fitting a model to these thin-layer drying kinetics (using variations of recorded product moisture content);

4. Validate the model using some thin-layer (using product moisture content) and/or deep-bed drying kinetics (using product moisture content and possibly output air moisture and temperature).

This strategy was used to identify the effective moisture diffusivity of P4Rs -2.5 mm. The validation set contained both thin-layer and deep-bed drying kinetics. As explained in a previous study [START_REF] Lambert | Reverse methodology to identify moisture diffusivity during air-drying of foodstuffs[END_REF], in this strategy, the mass variation of the product is the only experimental data used to identify the effective moisture diffusivity. However, using only mass-transfer-related information is problematic since, in drying, this transfer is strongly coupled to energy change. In a deep-bed experiment, though, inlet and outlet air temperatures and relative humidities may be recorded in order to provide energy-related data [START_REF] Lambert | Reverse methodology to identify moisture diffusivity during air-drying of foodstuffs[END_REF]. By providing data on both heat and mass transfers, an approach based on deepbed drying kinetics brings far more information to the identification procedure compared to those based on several measurements of thin-layer kinetics. With this in mind, the two other developed strategies were derived from the classical one, but with a major difference: effective moisture diffusivity was identified using deep-bed drying kinetics as opposed to thin-layer kinetics [START_REF] Lambert | Reverse methodology to identify moisture diffusivity during air-drying of foodstuffs[END_REF]. In the first approach, the identification set contained drying kinetics from a single deep-bed setup, obtained at a high temperature (90°C). The validation set contained both thin-layer drying kinetics and a set of deep-bed drying-cooling measurements of P4Cs. With the second approach, two sets of data were obtained for deep-bed kinetics, one at 40°C and one at 80°C, to identify effective moisture diffusivity of P4Rs -4 mm. The validation set was obtained from a thin-layer configuration.

For all identification strategies, operating conditions of the identification set are summarized in table 4 and operating conditions of the validation set are in table 3. The Abud correlation [START_REF] Abud-Archila | A compartmental model of thin layer drying kinetics of rough rice[END_REF] was used to find the effective moisture diffusivity (equation 4):

D ef f = exp (-d 1 + d 2 X T p) (4)
where d 1 and d 2 are constants that are optimized to force the simulation to best fit the experimental data.

Modeling of drying-cooling

In order to develop a versatile food-drying simulation, we selected the following compromise among the different approaches, published in a previous work (Lambert et al., 2015a), and considered some major assumptions. A convective drying model was created for both deep-bed and thin-layer configurations. One usually treats the deep bed as a series of thin layers for which the only transfer occurs between air and grain, and a thin layer as an equivalent single (average) particle. Moreover, deep-bed drying involves energy and mass transfers in the air, after each thin layer crossing [START_REF] Courtois | Dynamic Modelling of Drying to Improve Processing Quality of Corn[END_REF].

Heat and mass balances

For the particles

The following assumptions were considered for the particles:

Ideal geometry (infinite cylinder or spherical shape) with negligible shrinkage or swelling;

Isotropic behavior (heat and mass transfers only in the radial direction);

Negligible heat conduction between particles;

The heat and mass Biot values (Bi h and Bi m ) were, respectively, between 0.1 and 100 and much higher than 100. Therefore, heat conduction and liquid water diffusion within the particles were considered, along with water evaporation, heat convection, and mass convection at their surfaces [START_REF] Bird | Transport phenomena[END_REF]. The heat and mass Peclet values (P e h and P e m ) were much higher than 1, so a plug-type airflow and negligible pressure variation were assumed [START_REF] Bird | Transport phenomena[END_REF]. Heat and mass balances for the particles, liquid water diffusion (equation 5) and thermal conduction (equation 6) phenomena were based on Fick's and Fourier's law, respectively [START_REF] Bird | Transport phenomena[END_REF].

dX r, z t dt = 1 r n ∂ ∂r r n D ef f. ∂X r, z t ∂r (5) d (ρ dm (Cp dm + Cp w X r, z t ) T p r, z t ) dt = 1 r n ∂ ∂r r n λ ef f ∂T p r, z t ∂r (6)
where X is the moisture content (d.b.), r the radial position in the particle (m), z the particle position in the height of the bed (m), t the elapsed time (s), n the shape factor (1 or 2 for cylindrical or spherical symmetry, respectively), D ef f. the effective water diffusivity (m 2 .s -1 ), T p the product temperature (K), λ ef f. the effective thermal conductivity (W.m -1 .K -1 ), ρ dm the density of the dry matter (kg.m -3 ) and

Cp dm and Cp w the heat capacities at constant pressure of dry matter and water, respectively (J.kg -1 .K -1 ),

for the considered particle.

Temperature and moisture content within the particles were assumed to be initially uniform. In addition, boundary conditions are given in equations 7 to 10.

∂X 0, z t ∂r = 0 (7) ∂T p 0, z t ∂r = 0 (8) -ρ dm ∂ (D ap X rmax, z t ) ∂r = k m M w R a w P v sat. T p rmax,z t - HR P v sat T a z t = ϕ m (9) λ ap ∂T p rmax, z t ∂r = h glob (T a z t -T p rmax, z t ) -ϕ m ∆H v = ϕ Q (10)
where k m is the convective mass transfer coefficient (m.s -1 ), M w the molar mass of water (kg.mol -1 ), R the perfect gas constant (J.kg -1 .K -1 ), a w the water activity of the particle (P a/P a), P v sat the pressure of saturated vapor (P a), RH the relative humidity (decimal), h glob the external heat transfer coefficient including the heat transfer effectiveness by convection (W.m -2 .K -1 ), T a the air temperature (K), H v the enthalpy of vaporization of pure water (J.kg -1 ), ϕ m the mass flux density (kg.s -1 .m -2 ), and ϕ Q the heat flux density (J.s -1 .m -2 ).

For the air

A microscopic balance was created in the air after each thin-layer crossing (equations 11 and 12). Air moisture content and temperature were fixed at the inlet of the first thin layer but were not necessarily constant during a drying experiment.

d (ρ da Y z t ) dt = - ∂ (ρ da v a Y z t ) ∂z + ϕ m a (1 -ε) ε (11) d(ρ da (Cp da + Y z t Cp v ) T a z t ) dt = -∂ ∂z (ρ da (Cp da + Y z t Cp v ) T a z t ) - a (1-ε) ϕ Q ε (12)
where ρ da is the density (kg.m -3 ) of dry air, a the particle specific surface area (surface area per particle volume) (m -1 ), ε the bed porosity (decimal), Y the moisture content (d.b.), T a the temperature (K), and

Cp da and Cp v the specific heat capacities at a constant pressure of dry air and water vapor, respectively (J.kg -1 .K -1 ), for the considered air.

Simulation and optimization

The dynamic system consisted of a set of six Partial Differential Equations (PDEs 5 to 10). These equations were discretized in space into a set of Ordinary Differential Equations (ODEs) using an explicit finite volume scheme to be numerically integrated with respect to both time and space by the explicit embedded Runge-Kutta Cash-Karp (4, 5) method [START_REF] Galassi | Gnu scientific library reference manual[END_REF] (for stiff systems).

To quantify the quality of the model predictions, the maxima of relative errors between experimental and simulation data were calculated according to equations 13 and 14 for thin-layer and deep-bed drying kinetics, respectively, valid for any time t. The aim was to minimize the maximum relative error (MINMAX algorithm) obtained on the measurable variables using the Nelder and Mead method (also known as Downhill Simplex method) [START_REF] Nelder | A simplex method for function minimization[END_REF].

M RE T L = max   X t sim -X t exp X t exp 2   (13) M RE DB = max      X t sim -X t exp X t exp   2 , T a t out, sim -T t out, exp T a t out, exp 2    (14)
where M RE T L and M RE DB are the respective maximum relative errors for thin-layer and deep-bed kinetics, X sim1 and X sim2 the simulated mean moisture content (d.b.), X exp and X exp the experimental mean moisture content (d.b.) for thin-layer and deep-bed kinetics, respectively, and T a out, sim and T a out, exp the simulated and experimental output air temperature (°C).

The model simulation required 13 parameters, each with a physical meaning. Some were measured in the laboratory (pellet dimensions and density, porosity of the bed, initial temperature and moisture content of the pellets), some were found in the literature (heat capacity, thermal conductivity, latent heat of vaporization, coefficients of heat and mass transfers by convection), and others identified by reverse methods (effective moisture diffusivity and water activity). All parameters required for the simulation (parameters of the model and numerical parameters of the solver) are summarized in table 5.

Results and discussion

Sorption isotherms

The sorption isotherms of two pellet formulations (P4C or P4R -4 mm) are presented in figure 2. There was no visible effect of temperature within the [25°C -40°C] range (see figure 2.A). In addition, the standard deviation for the desorption and adsorption isotherm data of P4Rs at 25°C was less than 0.5% between the two repetitions, demonstrating the reproducibility of our measurements. By comparing the sorption isotherms for the two formulations, one can observe that (see figure 2.B):

when a w < 0.7, pellets for chickens were more hygroscopic than pellets for rabbits, meaning that the latter pellets had fewer sorption sites available at low water activity values; when a w > 0.7, pellets for rabbits were more hygroscopic than pellets for chickens, which is probably related to differences in the pellets' micro-porosity and composition. The water activity of mixtures is dependent on the fraction size of each ingredient [START_REF] Hemmingsen | Water adsorption in feed ingredients for animal pellets at different temperatures, particle size, and ingredient combinations[END_REF]. In particular, sugar beet pulp has a higher water-holding capacity than other ingredients at temperatures below 80°C [START_REF] Hemmingsen | Water adsorption in feed ingredients for animal pellets at different temperatures, particle size, and ingredient combinations[END_REF][START_REF] Serena | Chemical and physicochemical characterisation of co-products from the vegetable food and agro industries[END_REF]. Sugar beet pulp was present only in P4Rs, and was probably responsible for their more-hygroscopic behavior at high water activity values.

However, it should be noted that variations in water activity had little influence on the accuracy of the drying model. We therefore disregarded discrepancies between sorption isotherm values for water activities below 0.8 in the framework of the drying or drying-cooling model simulation. As stated previously, two different correlations were used to fit the desorption isotherm data: the Chung-Pfost correlation [START_REF] Chung | Adsorption and desorption of water vapor by cereal grains and their products. part i: Heat and free energy changes of adsorption and desorption[END_REF]) and a modified Henderson correlation [START_REF] Henderson | A basic concept of equilibrium moisture[END_REF][START_REF] Thompson | Mathematical simulation of corn drying: A new model[END_REF]. The former correlation was used by Brook [START_REF] Brook | Handbook of transportation and marketing in agriculture. chapter Drying, cleaning and condtionning[END_REF] to fit equilibrium moisture content and relative humidity data from soybeans. Maier et al. then showed that the identified parameters (p 1 , p 2 , p 3 and p 4 provided in equation 15) could be applied to different formulations of pellets for animal feed for relative humidity / water activity data between 20 and 90% (Maier & Bakker-Arkema, 1992).

X eq = 0.375 -0.0668 ln (-1.98 (Θ + 24.6) ln (a w ))

The identified parameters (p 1 , p 2 , p 3 and p 4 ) allowed for a good agreement with experimental data for water activities between 0.4 and 0.8, as shown in figure 3. However this correlation gives non-zero (negative) values of a w when X eq = 0, which is physically not acceptable and can cause numerical divergence during the simulations. water activity values given by the modified Henderson correlation were positive or zero when the equilibrium moisture content was set to zero, which was expected. Therefore this latter correlation was preferred for drying and drying-cooling simulations.

Identification of D ef f

The identified values of effective moisture diffusivity D ef f are given in table 6 for each formulation and geometry. In the case of pellets for rabbits (P4R -2.5 mm and P4R -4 mm) with a shared initial recipe, it was possible to consider a single set of coefficients (d 1 , d 2 ) for the effective moisture diffusivity of both geometries. However, generally speaking, the operating conditions of the pelleting stage drastically influenced the intrinsic thermophysical properties of the final product. Therefore it was preferable to identify the effective moisture diffusivity for each new formulation or set of pelleting operating conditions, which explains the slight differences observed between correlations of coefficients of effective moisture diffusivity (see table 6).

Final validation using drying and drying-cooling kinetics

Mean absolute errors observed on mean moisture content and outlet air temperature are given in table 7 using the kinetics of the validation set of each pellet formulation. Mean errors of prediction for product moisture content were close to the actual uncertainties of the experimental data. Unsurprisingly, the largest errors were obtained from drying-cooling experiments in which the initial state of the pellets was known with less certainty.

The drying-cooling process was well-simulated in the case of pellets for chickens (figure 4). More experiments would probably be needed for full validation, but these results confirm that the overall approach based on drying kinetics gives acceptable results in predicting drying-cooling experiments.

The simulations also yielded acceptable results for the drying of the P4R -4 mm pellets for rabbits (see figure 5) over the range of experimental conditions considered. The case of the P4R -2.5 mm pellets for rabbits was distinct from the others in that an identification set of thin-layer, instead of deep-bed, kinetics was used. At higher temperatures (in the 60 -80°C range; see figure 7), the model simulated with an acceptable degree of accuracy the deep-bed drying process for various air conditions (inlet air temperature in the range of 60 -90°C, inlet air velocity in the range of 0.7 -2 m.s -1 ) and bed heights (10 -18 cm) (see figure 6). However at lower temperatures (e.g. 40°C), an overall underestimation of pellet mean moisture content was observed (see left sub-figure in figure 7). The mean absolute error was below 0.032 d.b. Using thin-layer kinetics, the product mass variation was the only experimental data available. During the drying process, though, heat and mass transfers are strongly coupled. Hence, identification using one or two datasets of deep-bed kinetics brings far more information (pellet mass variation, inlet and outlet air temperatures and relative humidities)

and the identified parameters may be valid on a wider range of drying operating conditions. In addition, the identified D ef f parameters were also validated on classical drying-cooling operating conditions. Given the difficulties of recording drying-cooling kinetics, these results justify the decision to characterize formulations of animal feed pellets using only drying experiments.

Conclusion

Several formulations of pellets for animal feed were characterized experimentally and using correlations found in the literature. The sorption isotherms for all the different pellets yielded mostly similar results, close to those given previously by (Maier & Bakker-Arkema, 1992). The modified Henderson equation was preferred for numerical reasons and fitted onto the DVS data. To estimate the diffusion coefficient as a function of pellet moisture and temperature, several strategies were tested in the course of a larger project.

These were all based on inverse methods, i.e. to identify D ef f such that the simulated drying kinetics were as close as possible to experimental ones (identification set). With all parameters, the simulated drying kinetics, both in thin-layer and deep-bed arrangements, were in acceptable agreement with experimental ones (validation set), being close to the experimental uncertainties of the latter. Moreover, the model was proven to correctly simulate drying-cooling experiments. This approach is a possible workaround when dealing with the practical difficulties of performing reliable experiments on drying-cooling methods, especially at the thin-layer level. The same approach may be applied to different formulas of pellets with minimal laboratory work: measuring pellet size, adjusting the sorption isotherms (or neglecting their impact on the drying-

  204 d.b. (P4C), 0.211 d.b. (P4R -2.5 mm), and 0.202 d.b. (P4R -4 mm).

  inlet and outlet air temperatures (with T-type thermocouples and platinium probes) and relative humidities (with VAISALA 135Y and ROTRONIC HF532-HC2 probes); inlet air velocity (using an ANNUBAR ANR-73 Pitot tube and a FURNESS FCO-53 differential pressure sensor); mass of the product (with a METTLER PE 16 discontinuous weighing system with air bypassed).The precision of measurements of mean moisture content was estimated for each formulation of pellets by repeating measurements of thin-layer kinetics. The following uncertainties were determined: P4C: ± 0.015 d.b. (σ = 0.005 d.b.); P4R -2.5 mm and P4R -4 mm: ± 0.018 d.b. (σ = 0.06 d.b.).

X

  eq = exp ln [ -ln(1.0-aw) /20.802] identified parameters for the modified Henderson correlation (provided in equation 16) gave acceptable results, as shown in figure 3. More specifically, good agreement with experimental data was observed for water activity values below 0.8. Contrary to what was found with the Chung-Pfost correlation,

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: Sorption isotherms of pellets. (A) Effect of temperature on equilibrium moisture content, P4C pellets for chickens. (B) Effect of pellet formulation on equilibrium moisture content at 25°C. (→ adsorption curves), (← desorption curves).

Figure 5 :Figure 6 :

 56 Figure 5: Comparison of experimental (dots with error bars) and simulated (-) data of thin-layer drying kinetics of the P4R -4 mm pellets for rabbits (validation set)

X

  X Initial mean moisture content of the pellets respectively in the case of thin layer and deep bed kinetics T a Initial inlet air temperature va Inlet air velocity T p Initial temperature of the bed of pellets Table 7: Mean absolute errors observed on product mean moisture content and outlet air temperature of the kinetics of the validation set Mean moisture content of the particles (case of thin layer kinetics) X Mean moisture content of the particles within the bed (case of deep bed kinetics) T aout Outlet air temperature T a in Inlet air temperature va Inlet air velocity T p 0 Initial temperature of the bed of pellets

Table 2 :

 2 Dry matter composition of pellets for chickens and rabbits (P4R -2.5 mm and P4R -4 mm)

		Composition (% d.b.)
		Pellets for chickens	Pellets for rabbits
	Proteins	20.9	15.8
	Lipids	19.1	3.9
	Total sugar	54.5	70.5
	Ashes	5.5	9.8

Table 3 :

 3 Experimental conditions of the drying and the drying-cooling of pellets for chickens and rabbits

X is the mean moisture content of each layer and is calculated as the volume-weighted average moisture content.

X, the mean moisture content within the bed, is calculated as the mean of the mean moisture content of each thin layer.

Acknowledgments

The authors wish to thank several undergraduate, Master's, and Ph.D. students (Perrine Leclercq, Analia Aparecida Vanzo, Abir Hosni, Azziz Zemmouri, Monica Pinto, and Benoit Hareng) for their valuable help in acquiring the experimental data that supported the modeling and simulation work. The authors wish to thank the French Agence Nationale de la Recherche (ANR) for the funding of the Ph.D. of C. Lambert, and the partners of the COOPERE-2 project, led by Veolia Research and Innovation, in collaboration with ProSim and ENSIACET. The authors also wish to thank the French Agency ADEME (Agence Nationale de l'Environnement et de la Maîtrise de l'Energie) for funding the experimental characterization and modeling of the drying of pellets.


L Amini

R Sameni
email: reza.sameni@gipsa-lab.inpg.fr

C Jutten

G A Hossein-Zadeh

H Soltanian-Zadeh

MR ARTIFACT REDUCTION IN THE SIMULTANEOUS ACQUISITION OF EEG AND FMRI OF EPILEPTIC PATIENTS

Integrating high spatial resolution of functional magnetic resonance imaging (fMRI) and high temporal resolution of electroencephalogram (EEG) is promising in simultaneous EEG and fMRI analysis, especially for epileptic patients. The EEG recorded inside an MR scanner is interfered with MR artifacts. In this article, we propose new artifact reduction approaches and compare them with the conventional artifact reduction methods. Our proposed approaches are based on generalized eigenvalue decomposition (GEVD) and median filtering. The proposed methods are applied on experimental simultaneous EEG and fMRI recordings of an epileptic patient. The results show significant improvement over conventional MR artifact reduction methods.

INTRODUCTION

Although simultaneous acquisition of EEG and fMRI may be essential for exploring the dynamics and localization of neural activity, several artifacts interfere the resulting EEG. These artifacts include: 1) MR artifacts, 2) RF pulse artifacts, and 3) Balistocardiogram (BCG) artifacts. The MR artifact is caused by the switching of the magnetic filed gradients used in the image acquisition, which induces unwanted voltages on the EEG recordings. The RF pulse artifact is due to the time varying electromagnetic field pulses (RF pulses) used for excitation in MRI recordings [START_REF] Negishi | Origin of the Radio Frequency Pulse Artifact in Simultaneous EEG-fMRI Recording: Rectification at the Carbon-Metal Interface[END_REF]. Although, the RF pulse has a very high frequency (in the order of 10 to 100 megahertz); but it is nonlinearly rectified to low frequencies (below 100Hz), within the EEG bandwidth. BCG artifact is caused by the micro motions of head EEG leads and wires within the static magnetic field. These motions are related to the pulsatile blood-flow in the head.

Different methods have been proposed for removing MR artifacts. Sparse component decomposition on the wavelet and discrete cosine basis [START_REF] Xu | A New Method Based on Sparse Component Decomposition to Remove MRI Artifacts in the Continuous EEG Recordings[END_REF], blind source extraction (BSE) followed by averaging-and-subtraction [START_REF] Jing | Scanner Artifact Removal in Simultaneous EEG-fMRI for Epileptic Seizure Prediction[END_REF], and adaptive finite impulse response (FIR) filtering [START_REF] Wan | Artifact Reduction for Simultaneous EEG-fMRI Recording-Adaptive FIR Reduction of Imaging Artifacts[END_REF], are among these methods. Grouiller et al. [START_REF] Grouiller | A Comparative Study of Different Artefact Removal Algorithms for EEG Signals Acquired During Functional MRI[END_REF] have compared some common methods such as image artifact reduction (IAR) [START_REF] Allen | A Method for Removing Imaging Artifact from Continuous EEG Recorded During Functional MRI[END_REF], independent component analysis (ICA), fMRI artifact slice template removal (FASTR) [START_REF] Niazy | Removal of fMRI environment artifacts from EEG data using optimal basis sets[END_REF], and filtering in the frequency domain using the Fourier transform on simulated and real data. IAR [START_REF] Allen | A Method for Removing Imaging Artifact from Continuous EEG Recorded During Functional MRI[END_REF] is a conventional method based on subtraction of the averaging artifact waveform followed by adaptive noise cancellation (ANC). FASTR is the combination of principal component analysis (PCA) and IAR. Among the compared methods of [START_REF] Grouiller | A Comparative Study of Different Artefact Removal Algorithms for EEG Signals Acquired During Functional MRI[END_REF], IAR was shown to be the most effective method. This method is widely used for MR artifact removal in the literature.

In this article, we propose new approaches for MR artifact reduction. The proposed methods are based on generalized eigenvalue decomposition (GEVD) [START_REF] Sameni | Multichannel Electrocardiogram Decomposition using Periodic Component Analysis[END_REF], [START_REF] Parra | Blind Source Separation via Generalized Eigenvalue Decomposition[END_REF], [START_REF] Tong | Indeterminacy and Identifiability of Blind Identification[END_REF], [START_REF] Strang | Linear Algebra and Its Applications[END_REF], and median filtering. In fact, EEG signals having MR artifacts have both stationary and non-stationary properties at different time scales. MR artifacts include temporal structures that can be considered stationary. On the other hand, due to the switching of the magnetic field gradient, EEG contaminated with MR artifact is non-stationary in the time domain. In this work, two representations of MR artifact using this prior information have been utilized in an iterative GEVD scheme, together with a median and low-pass filter to extract the MR artifact. The results of the proposed method are compared with the IAR method over real simultaneous EEG and fMRI of epileptic patients.

The organization of the paper is as follows. Section 2, explains the GEVD and its application in source separation. The proposed MR artifact reduction approaches are presented in Section 3. Section 4 is devoted to the experimental results by reporting the results and quantitative comparison of the conventional IAR method and the proposed methods. Concluding remarks are presented in Section 5.

REVIEW OF GENERALIZED EIGENVALUE DECOMPOSITION

We assume zero-mean N-dimensional non-stationary observations x(t) that are mutually dependent in different N channels 1 . We are interested in linear mixtures of the form y(t) = w T x(t) that satisfy some measure of signal separability. For this, we define the following cost function:

J(w) . = E θ {y(θ )y(θ + τ θ )} E t {y 2 (t)} (1)
where E t {•} represents averaging over t, and θ and τ θ are respectively the time intervals and the time-varying lags that are found from our a priori knowledge of the signal's structure and non-stationarity. By maximizing this cost function, we are looking for the y(t) with bounded energy and with a maximal lagged correlation over specific time intervals (θ ) and time-lags (τ θ ).

Equation ( 1) can be rewritten as follows:

J(w) = w T B x w w T C x w (2)
where

C x . = E t {x(t)x(t) T } (3)
and

B x . = E θ {x(θ )x(θ + τ θ ) T } (4)
The matrix C x is the covariance matrix of x(t), which is known to be symmetric and positive definite. On the other hand, the matrix B x is not generally symmetric and we need to make it symmetric (B x ← (B x + B T x )/2) for further use 2 . Equation ( 2) is in the form of the Rayleight quotient [START_REF] Strang | Linear Algebra and Its Applications[END_REF], and its maximum value is achieved through the joint diagonalization of the matrix pair (B x ,C x ):

W T B x W = Λ W T C x W = I ( 5
)
where Λ is a diagonal matrix containing real generalized eigenvalues on its diagonal (in descending order), and W = [w 1 , ...w N ] is the matrix containing the corresponding generalized eigenvectors as its columns. The input signals may next be decomposed as follows:

y(t) = W T x(t) (6)
where y(t) = [y 1 (t), ..., y N (t)] T has decorrelated components, with the first component maximizing defined quotient in (2).

PROPOSED MR ARTIFACT REDUCTION APPROACHES

In Fig. 1, a typical segment of EEG contaminated with MR artifacts of two channels is depicted. MR artifacts are MR scanner dependent. MR artifact of the MR scanner used for our data are repeated with a period of three seconds, with two seconds of activation in each period (while magnetic field gradient was on). A closer look at the MR artifacts of each period, shows sharp regular peaks with dominant high amplitudes that are repeated every 62-63 samples (indicated with circles in Fig. 1). Each peak is also surrounded by several smaller peaks.

In the previous section, we noticed that the components extracted by GEVD are ranked according to an order that depends on the statistics carried by the matrix B x . The spatial whiteness of the extracted components is guaranteed by the diagonalization of C x . Here, our objective is to use the above mentioned properties of the MR artifact structure to form the matrix B x , such that it contains the statistical properties of these artifacts. Therefore, since the y i (t) calculated from [START_REF] Allen | A Method for Removing Imaging Artifact from Continuous EEG Recorded During Functional MRI[END_REF] are sorted in descending order of their corresponding eigenvalues, the MR artifacts are expected to be most concentrated in the first few components. In the following, two models of MR artifact (B x ) using this prior information are presented. a. First Statistical Measure

In the first approach, we first detect the dominant MR peaks from an arbitrary channel. Due to the stationarity of these peaks, a robust method for detecting the peaks is a matched filter that uses a typical MR artifact segment as its template filter. A window is considered around each sharp peak of MR artifact. The resulted windowed segment is assumed to be stationary. Hence the average of the correlation matrices of windowed segments of two proceeding peaks can be a good model of statistical properties of MR artifact. In other words since the MR artifacts of different channels are rather regular, the samples that are equally distanced from each peak are highly correlated with one another. We can therefore, calculate the correlation between the windows surrounding each peak as a statistical property of the MR artifacts. Carrying out this procedure over the different channels, i.e., by considering each time sample of the different channels as an N-dimensional vector, we can redefine the matrix B x defined in (4) as follows:

B x1 = E θ {x(θ )x(θ + ∆p i ) T } ( 7
)
where ∆p i = p i+1p i and θ ∈ {MR artifact time interval}. p i is the location of the i-th sharp peak of the MR artifact in the original data.

b. Second Statistical Measure

In the first model, the covariance matrix B x1 was calculated by direct vector-by-vector multiplication of the samples of short windows of x(θ ). According to Fig. 1, each of these windows, consists of a high-amplitude MR artifact peak and a few smaller peaks around it. Therefore, in the statistical measure B x1 , each sample will have a contribution as strong as its energy. In other words, the smaller peaks around the dominant peaks of MR artifacts will have less influence on the statistics of B x1 and the GEVD procedure will make less effort to remove such peaks from the input signal. However, the strong and weak MR artifacts are somehow equally destructive in EEG analysis. We can therefore attempt in equalizing the peaks before the calculation of our statistical measure. This idea leads us to the following definition of the MR artifact covariance matrix: of the windowed segments. Fig. 2 shows from top to bottom one row of F(θ ), a typical windowed segment and the corresponding weighted windowed segment, respectively. B x2 is the average of weighted windowed segments energies.

B x2 = E θ {[F(θ ) ⊗ x(θ )][F(θ ) ⊗ x(θ )] T } (8

Nonlinear Median Filtering

The GEVD procedure provides a means of transforming the input data, into components that are ranked according to their similarity with the MR artifacts. However, due to the linearity of the transformation ( 6), the performance of this primary stage is limited and it can not fully separate the MR artifacts from other components (including the EEG) except if the contaminated EEG signal perfectly satisfies a linear mixture (x = As + n), with a moderately low noise and a sufficient number of observation channels. Here the model of contaminated EEG signal is rather complex as proved by our obtained results of applying linear ICA. Here we propose the GEVD followed by a nonlinear denoising procedure that is applied to the first D (D < N) components of the GEVD outputs while keeping the rest of (N -D) components unchanged. This nonlinear denoising may indeed be applied to the original channels too (without the GEVD procedure) that is our third proposed approach (Fig. 3.c); but in this case we would not benefit from the mutual information of the different channels. In fact, since the GEVD procedure concentrates the energy of the MR artifacts in the first few components, the EEG (non-MR artifact) components are less influenced by the nonlinear filtering when we apply this filter at the output of the GEVD procedure.

In this study, a two-step moving window median (MWM) filter was used for nonlinear filtering. A theoretical study of the MWM filter and its extensions may be found in [START_REF] Arce | Nonlinear Signal Processing: A Statistical Approach[END_REF]. The first median filter highly attenuates the sharp peaks, while the second one (with a wider window length) reduces the remaining peaks and smoothes the results. As a rule of thumb, the median filter window lengths should be wider than the width of the MR artifacts. Therefore, with proper choice of median filter width, peaks such as the MR artifacts are less likely to pass the MWM filter and by subtracting the output of the second MWM filter from the input signal, we can assure that the MR artifact peaks are effectively removed, while the EEG components are least influenced. The MWM filter is further followed by a low-pass filter, to suppress the possible out of band components introduced by the nonlinear MWM filter and higher frequencies that are not interested in EEG analysis. After applying the filtering procedure denoted above (Median-LPF) to the first D components of GEVD, the D filtered and (N -D) unchanged components are back-projected to the input space using the inverse of the decomposing matrix W T (Fig. 3.a). The preprocessed data, denoted by X in Fig. 3, is the matrix including preprocessed EEG channels in each row. Preprocessed data is contaminated EEG with MR artifact without dc and linear trend. To detrend, the best straight-line fit is removed from each processed channel segment.

Iterative Denoising

Up to this point, we have improved the input signal quality by applying GEVD (using either B x1 or B x2 ), the nonlinear filtering procedure, and the back-projection. This procedure may be repeated in several iterations, each time over the output of the previous run. This leads to a deflation procedure that removes one (or more) dimensions of the MR artifacts in each iteration and is repeated until the residual signals are 'sufficiently clean'. For this, we need some measure of signal cleanness that can be used as the stopping criterion of the iterative procedure. Due to the switching magnetic field gradient, we propose to use the ratio between the signal energy of the MR artifact time intervals (θ ) over the rest of data, i.e. when the magnetic field gradient is off ( θ )):

C = E ch E θ { x2 ch (θ )} E θ { x2 ch ( θ )} (9)
where x(t) is the cleaned data of each channel in each iteration and E ch {•} is averaging over channels. The iterative procedure is stopped when the normalized difference of cleaning criteria between two iterations is less than a predefined threshold.

We hereby refer to this algorithm as Iterative GEVD. The overall flowchart of this algorithm is depicted in Fig. 3.a. Following, the two statistical measures defined in ( 7) and ( 8), the iterative procedure can either be based on B x1 , which we call Iterative GEVD-Similarities (IGS), or on B x2 , which we call Iterative GEVD-Weighted (IGW). The flowcharts of IGS and IGW are depicted in Fig. 3.b.

EXPERIMENTAL RESULTS

Data

The recordings were made in the 3T scanner (3T Bruker BioSpin, Bruker Medizintechnik GmbH, Ettlingen, Germany) utilizing an MR compatible EEG amplifier (SD32, Micromed, Treviso, Italy) with 17 c-shaped electrodes positioned according to the 10/20 system (O1 and O2 were not used for subjects comfort). The reference electrode was at Oz. A Gradient-Echo Echo Planar Imaging (GE-EPI) sequence used for MR image acquisition [START_REF] Grouiller | A Comparative Study of Different Artefact Removal Algorithms for EEG Signals Acquired During Functional MRI[END_REF]. The sampling rate of EEG acquisition was 1024Hz. An anti-aliasing hardware low-pass filter with cut-off frequency of 286.8Hz was applied. EEG signals were calibrated with a square wave of 100µV utilizing an external calibrator plugged on all inputs. The epileptic patient was required to keep eyes closed and relax during the experiment.

Results

Initially, the data were preprocessed for DC and trend removal. The three proposed approaches were applied on the preprocessed data. For the first two approaches (IGS and IGW), as explained before, the sharp peaks of the MR artifacts were detected using a matched filter (Fig. 1). Next, C x , B x1 and B x2 were calculated according to the procedures explained in previous sections. In the IGS method B x1 and C x were jointly diagonalized, while in the IGW approach B x2 and C x were jointly diagonalized. In each iteration of the proposed procedure, the first D = 3 of the N = 15 components extracted by GEVD, were denoised by the nonlinear MWM filter. The choice of D was based on empirical study of the output of the GEVD stage; the final results are not very sensitive to the value of D > 3. In the nonlinear filtering step, the first and second MWM window lengths used in the denoising were 20 and 40 samples (about 20ms and 40ms), respectively, since the typical width of sharp MR artifact peaks are about 20ms. The low-pass filter following the MWM filter was a first-order filter with a cut-off frequency of 45Hz, which is close to the effective bandwidth of the EEG. Next, the filtered and unchanged components were back-projected using the inverse of the decomposing matrix. This procedure was repeated until the normalized difference of the cleaning criteria defined in (9) between two consecutive iterations became less than 0.001. This limit was typically reached around 16 iterations. In Fig. 4, a typical segment of contaminated EEG and the results of IGS and IGW are shown. For comparison, the results of the Median-LPF and IAR [START_REF] Allen | A Method for Removing Imaging Artifact from Continuous EEG Recorded During Functional MRI[END_REF][START_REF] Grouiller | A Comparative Study of Different Artefact Removal Algorithms for EEG Signals Acquired During Functional MRI[END_REF] are added in Fig. 4.

Evaluation Methods

Fig. 4, compares the proposed methods with the IAR method presented in [START_REF] Allen | A Method for Removing Imaging Artifact from Continuous EEG Recorded During Functional MRI[END_REF][START_REF] Grouiller | A Comparative Study of Different Artefact Removal Algorithms for EEG Signals Acquired During Functional MRI[END_REF]. It is necessary to define an objective criterion for comparing the methods quantitatively. Since there is no ground truth, measuring the performance of each method is challenging. A first idea is to compare the signal with a reference which could be the signal acquired out of scanner or the time segment without MR artifact. The EEG signal inside and outside of MR scanner of a patient even with the same electrode locations are not necessarily similar. The EEG itself is very non-stationary in the time domain, especially for epileptic patients, since the IEDs (Interictal Epileptiform Discharges) appear spontaneously and the shape of IEDs are different. Due to this problem, the EEG contaminated with MR artifact and cleaned signal were compared. Comparable methods in the time domain were evaluated. The cleaned data were segmented around each peak of the MR artifact. The similarity between each two proceeding segments or all the segments was calculated. In other words, the periodicity of the cleaned data around each peak was measured. The cleaned data preserves the periodicity under each peak if the artifact removal method has not reduced MR artifact perfectly. In another method, the similarity between cleaned and original data under each sharp peak was measured. The cleaned data has some overshoots under each peak when the artifact removal method has not worked efficiently. However, since these methods of evaluation were not consistent for different MR artifact removal approaches, we compared spectral indexes of the EEG that are known to be more consistent. The following three criteria were defined for this purpose:

F 1 = 10log 10 E ch { E i { P( f p i )} E i {P( f p i )} } F 2 = 10log 10 E ch {E i { P( f p i ) P( f p i ) }} F 3 = 10log 10 E ch { E i=17,34 { P( f p i )} E i=17,34 {P( f p i )} } (10)
where E ch {•} and E i {•} are averages over channels and original data PSD peaks, respectively. P and P are the power spectral density (PSD) of the original and cleaned data, respectively. f p i is an interval of 4Hz around the peak of the original data PSD in frequencies higher than 14Hz. This range is chosen since the interested spikes or IEDs appear below 14Hz. Since the sampling rate is not a multiple of the MR artifact frequencies, the energy of MR artifact frequencies can be distributed among a few frequency bins of the PSD. Therefore, an interval of 4Hz is considered around each PSD peak.

F 1 and F 2 calculate the PSD attenuation (in dB) of cleaned data around each peak of the original data PSD (first row of Fig. 5). The method that most attenuates the PSD peaks of the MR artifact is preferred. Since sharp peaks are repeated in MR artifact every 62 or 63 samples and the sampling rate is 1024Hz, peaks around 16-17 Hz and their harmonics in the PSD of the original data is related to MR artifact. Therefore, F 3 measures the attenuation around peaks of 16Hz and 32Hz that are shown in Fig. 5. The peaks in higher harmonics are not considered in F 3 , because frequencies higher than 45Hz are not of interest in EEG analysis. PSD of contaminated EEG, EEG signal out of MR scanner, and the processed signal by IGS, IGW, Median-LPF, and IAR are depicted respectively from up to bottom in Fig. 5. To compute the PSD, the Discrete Fourier Transforms (DFT) of 0.5s windowed segments are calculated. To compare the results in high frequency resolution, the same criteria (F 1 , F 2 , and F 3 ) have been utilized, but the PSD is computed by the DFT of 9s windowed segments. The window length is chosen 9s since the gradient is switched on every 3 seconds. To measure how the cleaned data is following this period, several periods (e.g., 3 periods) are needed. F 1 , F 2 and F 3 in low and high frequency resolutions were calculated for IAR and the three proposed methods. It is important to mention that the low-pass filter applied on the results of the proposed methods utilized on the results of IAR method before F 1 -F 3 calculation for comparing all the methods equally. The calculations of F 1 -F 3 were carried out over 50 different 30s EEG segments for an epileptic patient. The average of 50 calculations of F 1 -F 3 are reported in Table 1. The method that better reduces the MR artifact has greater absolute values of F 1 -F 3 that is a necessary condition but not sufficient. Since the outer original signal has no MR artifact (no PSD peaks after 14Hz), F 1 -F 3 for the outer scanner signal has the maximum absolute values that is the minimum resemblance to contaminated EEG with MR artifact. Hence the F 1 -F 3 of the outer signal shows the minimum value in each column of Table 1.

By comparing the results of Table 1, it can be noticed that the proposed methods have outperformed the conventional IAR method. Among, the proposed methods first IGS, second IGW and third Median-LPF have better reduced the MR artifacts (using the defined criteria). In the IGS method, components of GEVD are well amplified in the first few components. But in the IGW the MR methods, MR artifact is not extracted as good as the IGS method. The Median-LPF method is very simple, but it does not consider the information between channels. It also does not use the prior information as much as the first two methods. The PSD attenuation in MR artifact peaks of the IGS method is comparable with the EEG acquired outside of the scanner. This method has effectively mitigated MR artifacts from contaminated EEG.

CONCLUSION

In simultaneous EEG and fMRI recording, MR artifact reduction is an important issue and the development of an efficient artifact removal algorithm has significant effect on further processing. We proposed three artifact removal methods using two models of MR artifact. The performance of our proposed artifact removal methods was evaluated on real epileptic patient data and appeared to outperform the conventional IAR method. The IGS method has shown better performance in MR artifact reduction compared to the other proposed methods. It is important to mention that as MR artifacts are MR scanner dependent, our method has been cus-tomized to the MR scanner used for our data. The presented algorithms were based on a linear model. Convolutive models will be studied in the future work for this purpose.

Figure 1 :

 1 Figure 1: Peak detection and windowing on MR artifact of EEG signal. Peaks are shown with circles.

Figure 2 :

 2 Figure 2: Weighting effect on windowed segments. From top to bottom, the window function, a typical windowed segment before and after weighting.

Figure 3 :

 3 Figure 3: Flowchart of proposed MR artifact reduction algorithms. (a) Iterative GEVD; (b) IGS and IGW, and (c) Median-low-pass filtering. D is the number of first few components in which the MR artifact is most amplified.

Figure 4 :

 4 Figure 4: Comparison of methods in reducing MR artifact from EEG signal. From up to bottom, contaminated EEG signal with MR artifact, the processed signal by IGS, IGW, Med-LPF, and IAR.

Figure 5 :

 5 Figure 5: Comparison of the PSDs. From up to bottom, PSD of contaminated EEG signal, outer scanner EEG signal, and the processed signal by IGS, IGW, Med-LPF, and IAR.

Table 1 :

 1 Frequency domain evaluation results (in dB) at low and high frequency resolutions. Orig-Out, IGS, IGW, Med-LPF, and IAR stand for original signal outside the MR scanner, Iterative GEVD-similarities, Iterative GEVD-weighted, Median-low-pass filter and Image Artifact Reduction[START_REF] Grouiller | A Comparative Study of Different Artefact Removal Algorithms for EEG Signals Acquired During Functional MRI[END_REF],[START_REF] Allen | A Method for Removing Imaging Artifact from Continuous EEG Recorded During Functional MRI[END_REF]

	Frequency	Low Resolution	High Resolution
	Method(dB)	F 1	F 2	F 3	F 1	F 2	F 3
	Orig-Out	-34 -33 -23 -35 -30 -25
	IGS	-28 -27 -16 -28 -24 -17
	IGW	-27 -25 -14 -26 -23 -15
	Med-LPF	-26 -25 -13 -25 -22 -14
	IAR	-21 -23 -11 -20 -21 -12

Such signals can for example be the output of a linear process of the form x(t) = As(t) + n where A ∈ R N×M

Note that here we are not concerned by the non-positive definiteness of B x ; since it is verified experimentally that most of the eigenvalues have positive values and negative eigenvalues have significantly smaller absolute values than the positive ones.

ACKNOWLEDGMENTS

We gratefully acknowledge Olivier David and Frédéric Grouiller for acquiring the data and the assistance at the beginning of this study.