Anass Bouchriti
email: anass.bouchriti@edu.uca.ac.ma

Morgan Pierre
email: morgan.pierre@math.univ-poitiers.fr

Nour Eddine Alaa
email: n.alaa@uca.ac.ma

Gradient stability of high

Keywords: BDF methods, gradient flow, Allen-Cahn equation, Lojasiewicz-Simon inequality

come

Introduction

We consider k-step backward differentiation formulae (BDF) applied to gradient flows in finite and infinite dimension. We are interested in the asymptotic behaviour, as time goes to infinity, of the sequences generated by these methods. Our model problem in the infinite dimensional case is the Allen-Cahn equation, which was introduced to describe the process of phase transition [START_REF] Allen | A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsing[END_REF]. It reads

u t -∆u + g(u) = 0, t > 0, x ∈ Ω, ( 1
)
where Ω is a bounded subset of R d with smooth boundary and g ∈ C 1 (R, R) is the derivative of a potential G. A standard choice is the double-well potential G(s) = (s 2 -1) 2 /4. The unknown function u represents the order parameter. Throughout the paper, we will consider Dirichlet boundary conditions.

The Allen-Cahn equation ( 1) is a gradient flow of the energy

E(u) = Ω 1 2 |∇u| 2 + G(u) dx (2)
for the L 2 (Ω) scalar product. The theoretical picture for (1) is well-known and there is a huge litterature on this equation (for the well-posedness and some asymptotic properties, see e.g. [START_REF] Haraux | Systèmes dynamiques dissipatifs et applications[END_REF][START_REF] Temam | Infinite-dimensional dynamical systems in mechanics and physics[END_REF]). If g is real analytic, then Simon [START_REF] Simon | Asymptotics for a class of nonlinear evolution equations, with applications to geometric problems[END_REF] proved that every bounded solution converges to a steady state as time goes to infinity. It is not known if this result still holds for every choice of nonlinearity g ∈ C ∞ (R, R) (see [START_REF] Haraux | Systèmes dynamiques dissipatifs et applications[END_REF][START_REF] Poláčik | Nonconvergent bounded solutions of semilinear heat equations on arbitrary domains[END_REF] and eqreferences therein). The proof of Simon is based on a so-called Lojasiewicz-Simon inequality, which is an extension to the infinite dimension of the celebrated Lojasiewicz inequality for real analytic functions [START_REF] Lojasiewicz | Ensembles semi-analytiques[END_REF]. We eqrefer the reader to the recent review by Haraux and Jendoubi [START_REF] Haraux | The convergence problem for dissipative autonomous systems[END_REF]. Our aim in this paper is to prove that convergence to equilibrium also holds for the BDF methods of high-order (3, 4 or 5) applied to equation [START_REF] Absil | Convergence of the iterates of descent methods for analytic cost functions[END_REF] with an analytic nonlinearity.

BDF methods of low order (1 or 2) are widely used for the time discretization of gradientflows such as [START_REF] Absil | Convergence of the iterates of descent methods for analytic cost functions[END_REF]. A crucial feature is that the discrete-in-time dynamical systems defined by these methods have properties similar to the continuous-in-time problem, such as the existence of a Lyapunov functional [START_REF] Shen | Numerical approximations of Allen-Cahn and Cahn-Hilliard equations[END_REF][START_REF] Stuart | Dynamical systems and numerical analysis[END_REF]. We note that for the BDF method of order 1, the Lyapunov functional is the energy E (2), whereas for higher-order BDF schemes, a pseudo-energy (i.e. a modification of E) has to be used.

The one-step BDF method (BDF1) is simply the backward Euler scheme, also known as the proximal algorithm in the context of optimization. Convergence to equilibrium for this scheme has been shown by means of a Lojasiewicz-like inequality in finite and infinite dimension [START_REF] Attouch | On the convergence of the proximal algorithm for nonsmooth functions involving analytic features[END_REF][START_REF] De Carvalho Bento | Dual descent methods as tension reduction systems[END_REF][START_REF] Merlet | Convergence to equilibrium for the backward Euler scheme and applications[END_REF]. The finite dimensional case arises naturally when one considers a space semi-discretization of (1) which preserves the gradient-flow structure. This happens with standard finite element or finite difference methods [START_REF] Chen | Applications of semi-implicit Fourier-spectral method to phase field equations[END_REF][START_REF] Elliott | The global dynamics of discrete semilinear parabolic equations[END_REF].

Convergence to equilibrium for the two-step BDF scheme (BDF2) was proved in finite dimension in [START_REF] Alaa | Convergence to equilibrium for discretized gradient-like systems with analytic features[END_REF], and in infinite dimension in [START_REF] Antonietti | Convergence to equilibrium for a secondorder time semi-discretization of the Cahn-Hilliard equation[END_REF] for the Cahn-Hilliard equation, which is a gradient flow of E (2) for the H -1 scalar product. We note that a restriction on the time step (independent of the initial data) has to be introduced in order for the gradient stability to hold. For the BDF1 scheme, the restriction on the time step can be removed by considering the proximal algorithm [START_REF] Bolte | Characterizations of Lojasiewicz inequalities: subgradient flows, talweg, convexity[END_REF]. The analysis of the infinite dimensional case is interesting because it is independent of the choice of a particular space discretization.

Since the gradient stability of the three-step BDF method is known [START_REF] Stuart | Dynamical systems and numerical analysis[END_REF], convergence to equilibrium can be obtained in finite dimension by applying general results on descent methods, see e.g. [START_REF] Alaa | Convergence to equilibrium for discretized gradient-like systems with analytic features[END_REF][START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods[END_REF][START_REF] De Carvalho Bento | Dual descent methods as tension reduction systems[END_REF]. The infinite dimensional case is more involved, but we work it out here in a way similary to the BDF2 case, with a stability restriction on the time step (Theorem 5.4). In contrast, there exist several second-order schemes which are unconditionnally gradient stable, such as the secant scheme, but for which convergence to equilibrium in the infinite dimensional setting remains to investigate (see [START_REF] Elliott | The global dynamics of discrete semilinear parabolic equations[END_REF][START_REF] Feng | Nonlinear stability of the implicit-explicit methods for the Allen-Cahn equation[END_REF][START_REF] Gomez | Provably unconditionally stable, second-order time-accurate, mixed variational methods for phase-field models[END_REF][START_REF] Shen | The scalar auxiliary variable (SAV) approach for gradient flows[END_REF][START_REF] Tierra | Numerical methods for solving the Cahn-Hilliard equation and its applicability to related energy-based models[END_REF] and references therein). If we consider for instance the secant scheme, which is a gradient-stable Crank-Nicolson like scheme, one issue is to obtain the precompactness of trajectories.

In [27, p. 424], Stuart and Humphries conjectured that the BDF schemes of order 4, 5 and 6 are gradient stable. To the best of our knowledge, this assertion has not yet been proved. In this paper, by introducing a notion of quadratic-stability, we prove that this conjecture is true for the BDF4 and BDF5 schemes. For k ≥ 7, the k-step BDF schemes are not zero-stable [START_REF] Hairer | Solving ordinary differential equations. I[END_REF], so they are not expected to be gradient stable. We make this clear in Remark 2.6. We note that the gradient-stability of the BDF6 scheme remains to be investigated (Remark 2.10).

The convergence result of Simon has been generalized to a great variety of gradient-like equations, in particular to damped wave equations such as the sine-Gordone equation or the modified Allen-Cahn equation (cf. [START_REF] Haraux | The convergence problem for dissipative autonomous systems[END_REF] and references therein). Convergence to equilibrium for the backward Euler equation applied to the modified Allen-Cahn equation has been proved in [START_REF] Grasselli | Convergence to equilibrium of solutions of the backward Euler scheme for asymptotically autonomous second-order gradient-like systems[END_REF][START_REF] Pierre | Convergence to equilibrium for a time semi-discrete damped wave equation[END_REF]. We note however that BDF schemes of order greater than one do not seem to preserve the gradient-like structure in this situation, so that other energy stable schemes should be considered.

The paper is organized as follows. In Section 2, we introduce the BDF methods applied to a gradient flow in finite dimension, and we define the notion of quadratic-stability. In Section 3, we prove that the BDF4 and BDF5 methods are gradient stable. Convergence to equilibrium for the BDF3 scheme applied to the Allen-Cahn equation is proved in Sections 4 and 5. In the last section, we perform some numerical simulations for the BDF3 scheme applied to [START_REF] Absil | Convergence of the iterates of descent methods for analytic cost functions[END_REF]. They nicely illustrate the theoretical results.

Quadratic-stability of BDF methods

The continuous problem

We denote •, • the scalar product in R M and • the euclidean norm. We consider the gradient flow in

R M U ′ (t) = -∇F (U (t)), t ≥ 0, (3)
where

U : [0, +∞) → R M and F ∈ C 1 (R M , R).
Throughout the paper, we assume that F satisfies the following two conditions,

∇F (V ) -∇F (W ), V -W ≥ -c F V -W 2 , ∀V, W ∈ R M , (4)
for some constant c F ≥ 0, and lim

V →∞ F (V ) = +∞. (5)
Condition ( 4) is known as a one-sided Lipschitz condition [START_REF] Stuart | Dynamical systems and numerical analysis[END_REF]. This condition is equivalent to the fact that the functional

V → ∇F (V ) + c F V is a maximal monotone operator on R M , or that the function V → F (V ) + (c F /2) V 2 is convex. Using the Cauchy-Peano theorem, this guarantees, for every U 0 ∈ R M , the existence of a unique solution U ∈ C 1 ([0, +∞), R M ) of (3) such U (0) = U 0 [8].
Condition ( 5) is called a coercivity condition. On multiplying equation ( 3) by U ′ (t), we obtain that

d dt F (U (t)) = -U ′ (t) 2 . (6)
In particular, F (U (t)) is nonincreasing, and the coercivity condition [START_REF] Attouch | On the convergence of the proximal algorithm for nonsmooth functions involving analytic features[END_REF] implies that the solution U is bounded. Using Lasalle's invariance principle [START_REF] Haraux | Systèmes dynamiques dissipatifs et applications[END_REF], we obtain that the ω-limit set of U (0), defined by

ω(U (0)) := {U ⋆ ∈ R M : ∃t n ↑ +∞ s. t. U (t n ) → U ⋆ },
is a compact and connected subset of R M , which is included in the set of critical points of F ,

S = {V ∈ R M : ∇F (V ) = 0}. (7)
If the cricital points of F are isolated, this implies that U (t) → U ⋆ as t → +∞. If the critical points of F are not isolated, then convergence to equilibrium may fail if M ≥ 2, even if F is C ∞ (see e.g. [START_REF] Absil | Convergence of the iterates of descent methods for analytic cost functions[END_REF]). Lojasiewicz [START_REF] Lojasiewicz | Ensembles semi-analytiques[END_REF] proved that if F : R M → R is real analytic, then we always have U (t) → U ⋆ as t → +∞.

General k-step BDF methods

Let ∆t > 0 be the time step. The general k-step backward differentiation formula (BDF) for ( 3) is defined by

k j=1 1 j ∂ j U n+k = -∆t∇F (U n+k ), n ≥ 0, (8)
where, for a sequence (U n ) n≥0 , the (backward) difference operator ∂ j is defined recursively by

∂ j U n = ∂ j-1 (U n -U n-1 ) (j ≥ 2, n ≥ j). When j = 1, we have ∂U n = U n -U n-1 .
The one-step BDF method is the backward Euler scheme:

U n+1 -U n = -∆t∇F (U n+1 ), n ≥ 0. (9)
The two-step BDF method reads

3 2 U n+2 -2U n+1 + 1 2 U n = -∆t∇F (U n+2 ), n ≥ 0. ( 10
)
The k-step BDF method (8) is obtained by Newton's divided difference formula which interpolates U by a polynomial of degree k (see Section 4.2 for the case k = 3). It is known to have a consistency error of order k. Moreover, it is zero-stable if and only if k ∈ {1, . . . , 6}, so that BDF methods are not used for k ≥ 7 [START_REF] Hairer | Solving ordinary differential equations. I[END_REF]. This also implies that if

F ∈ C k+2 (R M , R) (k ≤ 6)
and if the initial conditions are well chosen, the error between the solution U of (3) and its approximation given by the BDF scheme (8) is of order O(∆t k ) on finite time intervals [START_REF] Stuart | Dynamical systems and numerical analysis[END_REF]Theorem 3.5.7].

In [START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF], the initial conditions U 0 , . . . U k-1 are given in R M ; at every step n ≥ 0, the vector U n+k is computed from U n , U n+1 , . . . , U n+k-1 by solving the (generally nonlinear) equation [START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]. Thanks to the coercivity condition [START_REF] Attouch | On the convergence of the proximal algorithm for nonsmooth functions involving analytic features[END_REF], there exists at least one such U n+k . Indeed, [START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]

reads ∆t∇F (U n+k ) + α k U n+k + L n = 0, (11)
where

α k = k j=1 1/j > 0 ( 12
)
and L n ∈ R M is a linear combination of the vectors U n , . . . , U n+k-1 . The function

V → ∆tF (V ) + α k 2 V 2 + L n , V
is continous on R M and coercive. Thus, it has at least one minimizer U n+k in R M , which solves the Euler-Lagrange equation [START_REF] Elliott | The global dynamics of discrete semilinear parabolic equations[END_REF]. By using the one-sided Lipschitz condition (4), one can guarantee that the solution U n+k of ( 11) is unique for ∆t small enough (for a given L n ). Indeed, if Ũ n+k is another solution, then the difference

V n = U n+k -Ũ n+k satisfies ∆t[∇F (U n+k ) -∇F ( Ũ n+k ] + α k V n = 0.
On multiplying by V n and using (4), we see that

(α k -c F ∆t) V n 2 ≤ 0. ( 13
)
Thus, if c F ∆t < α k , V n = 0 and the sequence (U n ) n is uniquely defined by its initial values U 0 , . . . U n+k-1 and the k-step formula (8).

Quadratic-stability of BDF methods

In [START_REF] Stuart | Dynamical systems and numerical analysis[END_REF], the proof of gradient stability for the BDF methods of order 1, 2 and 3 is essentially based on the following algebraic computations. The authors multiply the left-hand side of (8) by ∂U n+k , and they consider the quantity

Γ k = k j=1 1 j ∂ j U n+k , ∂U n+k . ( 14
)
They use that

Γ 1 = ∂U n+1 2 , (15)
Γ 2 = ∂U n+2 2 + 1 4 ∂U n+2 -∂U n+1 2 + 1 4 ∂U n+2 2 -∂U n+1 2 , (16)
and

Γ 3 = 5 6 ∂U n+3 2 + 1 4 ∂U n+3 -∂U n+2 2 + 1 6 ∂U n+3 -∂U n+2 + ∂U n+1 2 + 5 12 ∂U n+3 2 -∂U n+2 2 + 1 6 ∂U n+3 -∂U n+2 2 -∂U n+2 -∂U n+1 2 . ( 17
)
In order to extend similar computations to higher-order schemes, it seems natural to consider Γ k as a quadratic form depending on the variables ∂U n+k , . . . , ∂U n+1 . For this purpose, we consider first the case M = 1 (the general case will follow immediately).

We set x 1 = ∂U n+k , . . . , x k = ∂U n+1 , and we define the following quadratic form on R k :

γ k (x 1 , . . . , x k ) = k j=1 1 j ( ∂j-1 x 1 )x 1 , ( 18
)
where ∂j is the forward difference operator defined recursively by ∂j x 1 = ∂j-1 (x 1 -x 2 ). We note that

∂j x 1 = j i=0 (-1) i j i x 1+i ,
where j i is the binomial coefficient.

We recall that a quadratic form q on R k is uniquely associated to a symmetric matrix

A q = (a ij ) 1≤i,j≤k of size k through q(x 1 , . . . , x k ) = k i=1 k j=1 a ij x i x j .
We will denote q ∼ A q this univoque correspondence. The quadratic form q is positive definite if q(x 1 , . . . , x k ) > 0 for all (x 1 , . . . , x k ) ∈ R k \ {0}, or equivalently, if the matrix A q is positive definite. By Sylvester's criterion, a symmetric matrix A = (a ij ) 1≤i,j≤k is positive definite if and only if its leading principal minors are positive, i.e. ∀p ∈ {1, . . . , k},

∆ p := det[(a ij ) 1≤i,j≤p ] > 0. ( 19
)
Definition 2.1. We say that the k-step BDF method is quadratic-stable if

γ k (x 1 , . . . , x k ) = q k (x 1 , . . . , x k-1 ) -q k (x 2 , . . . , x k ) + r k (x 1 , . . . , x k ), (20)
for all (x 1 , . . . , x k ) ∈ R k , where q k is a positive definite quadratic form on R k-1 and r k is a positive definite quadratic form on R k .

Example 2.2. From ( 16), we see that γ 2 is quadratic-stable with q 2 (x 1 ) = x 2 1 /4 and

r 2 (x 1 , x 2 ) = x 2 1 + (x 1 -x 2 ) 2 /4
. By [START_REF] Haraux | Systèmes dynamiques dissipatifs et applications[END_REF], γ 3 is also quadratic-stable (and by convention, so is γ 1 ).

Remark 2.3. The quadratic forms q k and r k , if they exist, are not unique. Indeed, let A ε qk be a symmetric matrix of size k which is a small perturbation of A qk in R k 2 , where q k ∼ A qk . Then A ε rk , defined through [START_REF] Lojasiewicz | Ensembles semi-analytiques[END_REF] in terms of A γk and A ε qk , is also a perturbation of A rk . By Sylvester's criterion [START_REF] Kavian | Introduction à la théorie des points critiques et applications aux problèmes elliptiques[END_REF], if A ε qk is close enough to A qk , then both A ε qk and A ε rk are positive definite.

We first show that quadratic-stability implies asymptotic stability. For this purpose, we apply the BDF method to the linear ode

y ′ (t) = -λy(t), t ≥ 0, (21)
with λ > 0. Since y(t) → 0 as t → ∞, we expect a similar behaviour for the time discrete version, which reads

k j=1 1 j ∂ j y n+k = -λ∆ty n+k , n ≥ 0. ( 22
)
Proposition 2.4. If the k-step BDF method is quadratic-stable, then for all λ > 0 and for any initial conditions y 0 , . . . , y n+k-1 , the sequence (y n ) n≥0 generated by [START_REF] Pierre | Convergence to equilibrium for a time semi-discrete damped wave equation[END_REF] tends to 0.

Proof. We multiply ( 22) by ∂y n+k , and by definition ( 18) of γ k , we obtain γ k (∂y n+k , . . . , ∂y n+1 ) = -λ∆ty n+k (y n+k -y n+k-1 ), n ≥ 0.

Using [START_REF] Lojasiewicz | Ensembles semi-analytiques[END_REF] and the well-known identity a(a

-b) = [a 2 -b 2 + (a -b) 2 ]/2, we obtain E n+k -E n+k-1 + λ 2 (∂y n+k ) 2 + r k (∂y n+k , . . . , ∂y n+1 ) = 0, n ≥ 0,
where

E n+k := λ 2 (y n+k ) 2 + q k (∂y n+k , . . . , ∂y n+2 ).
By summing on n, we find that

∞ n=0 r k (∂y n+k , . . . , ∂y n+1 ) ≤ E k-1 < +∞.
In particular, r k (∂y n+k , . . . , ∂y n+1 ) → 0 and since r k is positive definite, this shows that (∂y n+k , . . . , ∂y n+1 ) → (0, . . . , 0) in R k .

By [START_REF] Pierre | Convergence to equilibrium for a time semi-discrete damped wave equation[END_REF], y n+k is a linear combination of ∂y n+k ,. . . ,∂y n+1 . Thus, y n+k → 0 and the claim is proved.

In the following result, we use that asymptotic stability implies zero-stability.

Corollary 2.5. If k ≥ 7, the k-step BDF method is not quadratic-stable.

Proof. We consider now the BDF method applied to the ode y ′ (t) = 0, that is

k j=1 1 j ∂ j y n+k = 0, n ≥ 0. ( 23
)
This is a linear difference equation of order k. Its characteristic polymonial ρ(ζ), which is obtained on replacing y n+k by ζ n+k in ( 23) and dividing by ζ n , is

ρ(ζ) = k j=1 1 j j i=0 (-1) i j i ζ k-i = k j=1 1 j ζ k-j (ζ -1) j . ( 24
)
By definition, the k-step BDF method is zero-stable if and only if all the roots of ρ lie in the unit disc and the ones on the unit circle are simple. This is true if and only if k ∈ {1, . . . , 6} [START_REF] Hairer | Solving ordinary differential equations. I[END_REF].

Let now k ≥ 7. It is shown in [START_REF] Hairer | Solving ordinary differential equations. I[END_REF] that ρ has at least one root λ k ∈ C outside the unit disc, i.e. |λ k | > 1. For ε > 0, we consider now the polynomial

ρ ε (ζ) = ρ(ζ) + εζ k .
This is the characteristic polynomial of the linear difference equation ( 22) with ε = λ∆t. As ε → 0 + , the roots of ρ ε converge to the roots of ρ, so that for ε > 0 small enough, ρ ε has at leat one root λ k,ε such that |λ k,ε | > 1. Therefore, we have a sequence (λ n k,ε ) n≥0 of complex numbers which complies with [START_REF] Pierre | Convergence to equilibrium for a time semi-discrete damped wave equation[END_REF], and which is unbounded. By Proposition 2.4, this k-step BDF method cannot be quadratic-stable. We note here that the result in Proposition 2.4 is proved for a sequence (y n ) n of real numbers, but it extends immediately to complex numbers, by considering the real and imaginary parts.

Remark 2.6. The proof above shows that for k ≥ 7, the k-step BDF schemes are not gradient stable in the sense of [START_REF] Stuart | Dynamical systems and numerical analysis[END_REF]Definition 1.8.9]. Indeed, consider the function F (y) = λy 2 /2 with λ > 0, so that ( 21) is the gradient flow of F . Assume by contradiction that a k-step BDF scheme is gradient stable for some k ≥ 7. Then every sequence (y n ) n defined by ( 22) would be bounded, by [START_REF] Stuart | Dynamical systems and numerical analysis[END_REF]Theorem 1.8.10]. But we have seen that for ∆t > 0 small enough, this k-step BDF scheme produces an unbounded sequence.

Proposition 2.7. The 4-step BDF method is quadratic-stable. Namely,

γ 4 (x 1 , x 2 , x 3 , x 4 ) = q 4 (x 1 , x 2 , x 3 ) -q 4 (x 2 , x 3 , x 4 ) + r 4 (x 1 , x 2 , x 3 , x 4 ), for all (x 1 , x 2 , x 3 , x 4 ) ∈ R 4 , where q 4 ∼ 1 144   225 -
    . (25)
Proof. We give a constructive proof. We seek A q4 in the form of a symmetric matrix

A q4 =   a d g d b f g f c   . The quadratic form q4 (x 1 , x 2 , x 3 , x 4 ) = q 4 (x 1 , x 2 , x 3 ) -q 4 (x 2 , x 3 , x 4 )
is given by

A q4 =     a d g 0 d b -a f -d -g g f -d c -b -f 0 -g -f -c     .
We note that the sum of the elements in each diagonal of A q4 amounts to 0. Since A r4 = A γ4 -A q4 , this shows that the sum of the elements in each diagonal of A r4 and A γ4 are equal.

By definition, we have

A γ4 =     25
γ 5 (x 1 , x 2 , x 3 , x 4 , x 5 ) = q 5 (x 1 , x 2 , x 3 , x 4 ) -q 5 (x 2 , x 3 , x 4 , x 5 ) + r 5 (x 1 , x 2 , x 3 , x 4 , x 5 ), (26)
for all (x 1 , x 2 , x 3 , x 4 , x 5 ) ∈ R 5 , where

q 5 ∼     1.
      . ( 28
)
Proof. By definition, we have

γ 5 ∼ 1 120       274 -163 137 -63 12 -163 0 0 0 0 137 0 0 0 0 -63 0 0 0 0 12 0 0 0 0      
, so it is easily seen that ( 26) holds, with q 5 and r 5 as above. Using Sylvester's criterion, we can check numerically that q 5 and r 5 are positive definite (the computation is exact since we work with rational numbers of reasonable size), and the proof is complete (see Remark 2.9 for a constructive approach).

Remark 2.9. The argument used in Proposition 2.7 is not sufficient in order to obtain q 5 and r 5 . We have actually obtained q 5 by numerical investigation with a Scilab program. More precisely, let

q 5 ∼ A q5 = (a i,j ) 1≤i,j≤4 ,
where the ten coefficients a i,j for 1 ≤ i ≤ j ≤ 4 are unknown (the other coefficients are obtained by symmetry of A q5 ). The idea is to check a large number of values of the coefficients a i,j in order to obtain the positivity of q 5 and r 5 . We have used a grid of stepsize 0.15 in R 10 for these coefficients. Remark 2.3 shows that for a stepsize small enough, this approach should give a result, if γ 5 is indeed quadratic stable.

In order to reduce the computational time, we used the following simple facts. The diagonal coefficients of A r5 being positive, we deduce from (26) that 274/120 > a 1,1 > a 2,2 > a 3,3 > a 4,4 > 0.

For each nondiagonal coefficient a i,j , we note that the minor a i,i a j,j -a 2 i,j corresponding to q 5 is positive, that is |a i,j | < √ a i,i a j,j . For every choice of A q5 in our program, we checked wheter the matrices A q5 and A r5 were positive definite or not by a numerical calculation of there eigenvalues.

Remark 2.10. We have investigated numerically the quadratic-stability of γ 6 in a manner similar to γ 5 (cf. Remark 2.9). We used a grid of stepsize 0.1 in R 15 for the fifteen unknown coefficients of A q6 . We have not been successful. It could be interesting to perform numerical simulations with a smaller stepsize, but we were limited by the computational time.

3. Gradient stability of the BDF4 and BDF5 methods

3.1. From R to R M
Now, we turn back to the general case M ≥ 1. If

q(x 1 , . . . , x k ) = k i=1 k j=1 a ij x i x j
is a quadratic function on R k , associated to the symmetric matrix A q = (a ij ) 1≤i,j≤k , we define the following quadratic form on (R M ) k :

Q(V 1 , . . . , V k ) = k i=1 k j=1 a ij V i , V j .
Then Q inherits the properties of q. In particular, if q is positive definite, then so is Q. Indeed, if A q is positive definite, there exists a matrix P ∈ O(k) such that

A q = P t      λ 1 0 • • • 0 0 λ 2 . . . . . . . . . . . . . . . 0 0 • • • 0 λ k      P, where 0 < λ 1 ≤ λ 2 ≤ • • • λ k are the eigenvalues of A q . This reads, for all X = (x 1 , . . . , x k ) t ∈ R k , q(x 1 , . . . , x k ) = X t A q X = k i=1 λ i ( k j=1 p ij x j ) 2 ,
where

P = (p ij ) 1≤i,j≤k . Thus, for (V 1 , . . . , V k ) ∈ (R M ) k , we have Q(V 1 , . . . , V k ) = k i=1 λ i k j=1 p ij V j 2 ≥ λ 1 k i=1 k j=1 p ij V j 2 ≥ λ 1 k i=1 k j=1 k m=1 p ij p im V j , V m . (29)
Since P t P = Id k , we have k i=1 p ij p im = δ jm , where δ jm is the Kronecker symbol. Thus, (29) reads

Q(V 1 , . . . , V k ) ≥ λ 1 k j=1 V j 2 ,
and the claim is proved.

Gradient stability of the BDF4 method

From Proposition 2.7, we deduce:

Lemma 3.1. The quadratic forms Γ 4 , Q 4 and R 4 associated to γ 4 , q 4 and r 4 satisfy

Γ 4 (V 1 , V 2 , V 3 , V 4 ) = Q 4 (V 1 , V 2 , V 3 ) -Q 4 (V 2 , V 3 , V 4 ) + R 4 (V 1 , V 2 , V 3 , V 4 ) ( 30
)
for all (V 1 , V 2 , V 3 , V 4 ) ∈ (R M ) 4 .
In particular, Q 4 and R 4 are positive definite. Moreover,

R 4 (V 1 , V 2 , V 3 , V 4 ) = 41 144 V 1 2 + R4 (V 1 , V 2 , V 3 , V 4 ), ( 31
)
where R4 is positive definite on (R M ) 4 .

Proof. It remains only to prove that R4 is positive definite. For this purpose, we consider the quadratic form on R 4 defined by

r4 (x 1 , x 2 , x 3 , x 4 ) = r 4 (x 1 , x 2 , x 3 , x 4 ) - 41 144 x 2 1 .
Using Sylvester's criterion, we easily check that A r4 is positive definite. The claim on R4 follows. We point out that the coefficient 41/144 was found by numerical investigation and that it is not optimal.

For V = (V 0 , V 1 , V 2 , V 3 ) ∈ (R M ) 4 , we define F4 ( V ) = F (V 0 ) + 1 ∆t Q 4 (V 1 , V 2 , V 3 ),
and for a sequence (U n ) n≥0 in R M , we define

Û n+4 = (U n+4 , ∂U n+4 , ∂U n+3 , ∂U n+2 ), so that F4 ( Û n+4 ) = F (U n+4 ) + 1 ∆t Q 4 (∂U n+4 , ∂U n+3 , ∂U n+2 ).
Theorem 3.2. If c F ∆t ≤ 41/72, the 4-step BDF method (8) is gradient-stable, i.e.

F4 ( Û n+4 ) + 1 ∆t R4 (∂U n+4 , ∂U n+3 , ∂U n+2 , ∂U n+1 ) ≤ F4 ( Û n+3 ), n ≥ 0. ( 32
)
Proof. We multiply (8) by ∂U n+4 in R M , and we find

Γ 4 (∂U n+4 , ∂U n+3 , ∂U n+2 , ∂U n+1 ) = -∆t ∇F (U n+4 ), U n+4 -U n+3 ,
where Γ 4 is defined by [START_REF] Gomez | Provably unconditionally stable, second-order time-accurate, mixed variational methods for phase-field models[END_REF]. Since V → F (V ) + (c F /2) V 2 is convex, we have

F (V ) -F (U ) + c F 2 V -U 2 ≥ ∇F (U ), V -U , ∀U, V ∈ R M . ( 33
)
We apply this relation with U = U n+4 and V = U n+3 and we use [START_REF] Tierra | Numerical methods for solving the Cahn-Hilliard equation and its applicability to related energy-based models[END_REF]. This yields

F4 ( Û n+4 ) + 1 ∆t R 4 (∂U n+4 , ∂U n+3 , ∂U n+2 , ∂U n+1 ) ≤ F4 ( Û n+3 ) + c F 2 ∂U n+4 2 ,
for all n ≥ 0. Finally, using (31) and 1 ∆t 41 144 -cF 2 ≥ 0, we find (32). Let

ω((U n ) n ) = {U ⋆ ∈ R M : ∃n p → +∞ s.t. U np → U ⋆ }
denote the ω-limit set of a sequence (U n ) n , and recall that S denotes the set of critical point of F (cf. ( 7)).

From the gradient stability of the BDF4 method, we infer:

Corollary 3.3.
Assume that c F ∆t ≤ 41/72, and let (U n ) be a sequence in R M defined by the 4-step BDF method. Then ω((U n ) n ) is a compact and connected subset of R M included in S.

Proof. Since F is continuous on R M and coercive (assumption ( 5)), F is bounded from below, i.e. there exists a real number m such that

F (V ) ≥ m, ∀V ∈ R M .
By (32), the sequence F4 ( Û n+3 ) is nonincreasing; it is also bounded from below (by m), so its converges to a value F ⋆ ∈ R. By induction, we obtain from (32) that for all N ≥ 0,

F4 ( Û N +4 ) + 1 ∆t N n=0 R4 (∂U n+4 , ∂U n+3 , ∂U n+2 , ∂U n+1 ) ≤ F4 ( Û 3 ).
Letting N tend to +∞, we find that

1 ∆t +∞ n=0 R4 (∂U n+4 , ∂U n+3 , ∂U n+2 , ∂U n+1 ) ≤ F4 ( Û 3 ) -m < +∞.
In particular, R4 (∂U n+4 , ∂U n+3 , ∂U n+2 , ∂U n+1 ) → 0. Since R4 is positive definite, we know that

λ ∂U n+4 2 + ∂U n+3 2 + ∂U n+2 2 + ∂U n+1 2 ≤ R4 (∂U n+4 , ∂U n+3 , ∂U n+2 , ∂U n+1 ),
for some constant λ > 0 independent of n, and so ∂U n+1 → 0. The boundedness of F4 ( Û n+3 ) also implies that (U n ) n is bounded. A standard argument shows then that ω

((U n ) n ) is a compact and connected subset of R M . If U ⋆ = lim p→+∞ U np+4 is an element of ω((U n ) n ),
then the right hand-side of (8) (with k = 4) tends to -∆t∇F (U ⋆ ), while the left hand-side tends to 0. Thus U ⋆ ∈ S, and the proof is complete.

Corollary 3.4. Assume that c F ∆t ≤ 41/72, and let (U n ) be a sequence in R M defined by the 4-step BDF method. Assume moreover that either one of the following assumptions hold: i) the space dimension is M = 1; i) the critical points of F are isolated; ii) the function F is real analytic on R M .

Then the whole sequence (U n ) n converges to a single point U ⋆ which is a critical point of F .

Proof. If i) holds, a proof similar to [15, Corollary 2] can be carried out. It uses that R is an ordered field. Details are left to the reader.

If assumption ii) holds, then the result is a direct consequence of Corollary 3.3.

If assumption iii) holds, the claim follows from general results on descent methods. We refer for instance to [START_REF] Alaa | Convergence to equilibrium for discretized gradient-like systems with analytic features[END_REF]Theorem 2.4]. In this case, the crucial point is the Lojasiewicz inequality for real analytic functions: it applies here to the Lyapunov function F4 .

Gradient stability of the BDF5 method

From Proposition 2.8, we deduce: Lemma 3.5. The quadratic forms Γ 5 , Q 5 and R 5 associated to γ 5 , q 5 and r 5 satisfy 5 . In particular, Q 5 and R 5 are positive definite. Moreover,

Γ 5 (V 1 , V 2 , V 3 , V 4 , V 5 ) = Q 5 (V 1 , V 2 , V 3 , V 4 ) -Q 5 (V 2 , V 3 , V 4 , V 5 ) + R 5 (V 1 , V 2 , V 3 , V 4 , V 5 ) (34) for all (V 1 , V 2 , V 3 , V 4 , V 5 ) ∈ (R M )
R 5 (V 1 , V 2 , V 3 , V 4 , V 5 ) = 1 20 V 1 2 + R5 (V 1 , V 2 , V 3 , V 4 , V 5 ), ( 35
)
where R5 is positive definite on (R M ) 5 .

Proof. It remains only to prove that R5 is positive definite. For this purpose, we consider the quadratic form on R 5 defined by

r5 (x 1 , x 2 , x 3 , x 4 , x 5 ) = r 5 (x 1 , x 2 , x 3 , x 4 , x 5 ) - 1 20 x 2 1 .
Using Sylvester's criterion, we easily check that A r5 is positive definite. The claim on R5 follows. As previously, the coefficient 1/20 was found by numerical investigation and it is not optimal.

For V = (V 0 , V 1 , V 2 , V 3 , V 4 ) ∈ (R M ) 5 , we define F5 ( V ) = F (V 0 ) + 1 ∆t Q 5 (V 1 , V 2 , V 3 , V 4 ),
and for a sequence (U n ) n≥0 in R M , we define

Û n+5 = (U n+5 , ∂U n+5 , ∂U n+4 , ∂U n+3 , ∂U n+2 ), so that F5 ( Û n+5 ) = F (U n+5 ) + 1 ∆t Q 5 (∂U n+5 , ∂U n+4 , ∂U n+3 , ∂U n+2 ).
By arguing as in the proof of Theorem 3.2 and using Lemma 3.5, we find: Theorem 3.6. If c F ∆t ≤ 1/10, the 5-step BDF method (8) is gradient-stable, i.e.

F5 ( Û n+5 ) + 1 ∆t R5 (∂U n+5 , ∂U n+4 , ∂U n+3 , ∂U n+2 , ∂U n+1 ) ≤ F5 ( Û n+4 ), n ≥ 0. ( 36
)
As a consequence of the gradient stability of BDF5 method, Corollaries 3.3 and 3.4 are also valid for a sequence (U n ) n generated by the BDF5 scheme (assuming the stability condition c F ∆t ≤ 1/10 for the time step).

In Table 1, we sum up the results obtained so far regarding uniqueness and gradient stability of k-step BDF schemes. Uniqueness is ensured if c F ∆t < α k (cf. ( 12)-( 13)) and gradient stability is ensured if c F ∆t < β k . The values of β k (1 ≤ k ≤ 3) obtained in [START_REF] Stuart | Dynamical systems and numerical analysis[END_REF] have been improved here by a factor two, thanks to (33). For k = 1, we can have β 1 = +∞ if the BDF1 scheme is the proximal algorithm [START_REF] Attouch | On the convergence of the proximal algorithm for nonsmooth functions involving analytic features[END_REF][START_REF] Merlet | Convergence to equilibrium for the backward Euler scheme and applications[END_REF]. For k = 1, 2, we have α k < β k : a particular choice of initial values may produce several gradient stable sequences. In contrast, for 3 ≤ k ≤ 5, we have α k > β k , so that provably gradient stable sequences are uniquely defined. We note that the values of β 4 and β 5 are not optimal (cf. Lemmas 3.1 and 3.5), and we do not know if β 3 is optimal. If the function F is convex, then c F = 0 and the BDF methods are unconditionnally stable for 1 ≤ k ≤ 5. 4. The BDF3 scheme applied to the Allen-Cahn equation

Notation and assumptions

Let H = L 2 (Ω) be equipped with the L 2 (Ω) norm |•| 0 and product (•, •). We denote V = H 1 0 (Ω) the standard Sobolev space, and

|•| 1 = |∇•| 0 the associated hilbertian norm. Recall that -∆ : V -→ V ′ is an isomorphism associated to the inner product on V through < -∆u, v > V ′ ×V = (∇u, ∇v) ∀u, v ∈ V.
We assume that the nonlinearity g : R → R in ( 1) is real analytic on R and if d ≥ 2, we assume that there exist c > 0 and p ≥ 0 such that

(d -2)p < 4 and |g ′ (s)| ≤ c(1 + |s| p ) ∀s ∈ R. ( 37
)
No growth assumption is required if d = 1. Moreover, there exists a constant c g ≥ 0 such that

g ′ (s) ≥ -c g ∀s ∈ R. ( 38
)
The last assumption is a coercivity condition which yields

lim |s|→∞ inf g(s) s > -λ 1 , (39)
where λ 1 > 0 is the first eigenvalue of -∆, that is

λ 1 = inf v∈V \{0} |v| 2 1 |v| 2 0 . (40)
The Poincaré inequality implies that λ 1 > 0. We note that H is continuously imbedded in V ′ and for all v ∈ H, we have

v V ′ := sup |w|1≤1 v, w V ′ ×V = sup |w|1≤1 (v, w) ≤ sup |w|1≤1 |v| 0 |w| 0 ≤ 1 √ λ 1 |v| 0 . ( 41
)
By the mean value theorem, we have

(g(b) -g(a))(a -b) ≥ -c f (a -b) 2 , ∀a, b ∈ R. ( 42
)
Conversely, this relation implies (38), which can therefore be considered as a one-sided Lipschitz condition on g (cf. ( 4)).

We define G as the functional

G(s) = s 0 g(r)dr. Assumption (38) implies that G is semi- convex, i.e. s → G(s) + (c g /2)s 2 is convex.
We note that the standard choice g(s) = s 3 -s satisfies (37)-(39) in space dimension 1 ≤ d ≤ 3.

Derivation of the BDF3 scheme

Consider a uniform grid on [0, +∞) with a positive stepsize ∆t and times t n = n∆t (n = 0, 1, . . . ). We shall compute u(t) at these points, and we denote (momentarily) u n = u(t n ). Consider the third degree polynomial P u interpolating u n , u n-1 , u n-2 and u n-3 at times t n , t n-1 , t n-2 and t n-3 respectively,

P u (t) = 3 k=0 3 j=0 j =k t -t n-j t n-k -t n-j u n-k .
Since the grid is uniform, this interpolation polynomial P u can also be written in Newton form (divided differences),

P u (t) = u n + 3 k=1 1 ∆t k k! k-1 j=0 (t -t n-j ) ∂ k u n , where ∂ k u n is recursively defined by ∂ k u n = ∂ k-1 (∂u n ) = ∂ k-1 (u n -u n-1 ), for k ≥ 1.
If the last equation is differentiated with respect to t, one gets

d dt P u (t) = 3 k=1 1 ∆t k k!   d dt k-1 j=0 (t -t n-j )   ∂ k u n = 3 k=1 1 ∆t k k!    k-1 p=0 k-1 j=0 j =p (t -t n-j )    ∂ k u n ,
that is

d dt P u (t) = 3 k=1 1 ∆t k k!       k-1 p=1 k-1 j=0 j =p (t -t n-j ) + k-1 j=1 (t -t n-j )       ∂ k u n . P (t)
Since t n is a root of the polynomial P , by setting t = t n , we obtain

d dt P u (t n ) = 3 k=1 1 ∆t k k!   k-1 j=1 (t n -t n-j )   ∂ k u n = 3 k=1 1 ∆t k k!   k-1 j=1 j∆t   ∂ k u n .
Thus, we have

d dt P u (t n ) = 1 ∆t 3 k=1 1 k ∂ k u n . (43)
Substituting ( 43) into (1) yields the 3-step backward differentiation formula applied to the Allen-Cahn equation,

1 ∆t 3 k=1 1 k ∂ k u n -∆u n + g(u n ) = 0, ( 44
)
that is,

1 ∆t 11 6 u n -3u n-1 + 3 2 u n-2 - 1 3 u n-3 -∆u n + g(u n ) = 0.
We note that (44) holds only up to order (∆t) 3 if u n = u(t n ), so that the sequence (u n ) n generated by (44) only defines an approximation of the values (u(t n )) n . The third order accuracy of the BDF3 scheme is well-known in the finite dimensional case [START_REF] Stuart | Dynamical systems and numerical analysis[END_REF]. It could also be carried out in our situation by using standard techniques and by assuming enough regularity on the solution of (1). We refer to [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF] for an error analysis of a fully discrete version of (1) based on the BDF1 (i.e. backward Euler) scheme for the time discretization; see also [START_REF] Shen | Numerical approximations of Allen-Cahn and Cahn-Hilliard equations[END_REF] for an error analysis in the case of a second-order scheme. We shall study the scheme in the following form: let u 0 , u 1 , u 2 ∈ V , and for n = 2, 3, . . . , let u n+1 ∈ V solve

u n+1 - 18 11 u n + 9 11 u n-1 - 2 11 u n-2 - 6 11 ∆t ∆u n+1 + 6 11 ∆t g(u n+1 ) = 0 (45) in V ′ .

Existence, uniqueness and stability

Proposition 4.1 (Existence). For any (u 0 , u 1 , u 2 ) ∈ V 3 , there exists at least one sequence (u n ) n≥0 in V which complies with (45).

Proof. Let u n-2 , u n-1 , u n be fixed in V . Then u n+1 can be obtained as a minimizer of the functional

G n : V -→ R u -→ 1 2 |u| 2 0 + 1 11 (u, -18u n + 9u n-1 -2u n-2 ) + 3 11 ∆t |u| 2 1 + 6 11 ∆t(G(u), 1).
The growth restriction (37) and Sobolev imbeddings imply that the functional u → (G(u), 1) is well-defined and of classe C 1 on V [START_REF] Kavian | Introduction à la théorie des points critiques et applications aux problèmes elliptiques[END_REF]. The coercivity condition (39) implies the existence of constants k 1 < λ 1 and k 2 ≥ 0 such that, for all s in R,

G(s) ≥ - k 1 2 s 2 -k 2 .
By integration on Ω,

(G(u), 1) ≥ - k 1 2 |u| 2 0 -k 2 |Ω| ≥ - k 1 2λ 1 |u| 2 1 -k 2 |Ω|, (46)
where in the last inequality, we use the Poincaré inequality (40). A direct consequence is that

G n (u) ≥ 1 4 |u| 2 0 + 1 - k 1 λ 1 3 11 ∆t |u| 2 1 -c n ( 47
)
for some constant c n ≥ 0 independent of u ∈ V . In particular, G n is bounded from below. Denote then m = inf v∈V G n (v) and consider a minimizing sequence (v k ) k≥0 ⊂ V . By (47), (v k ) k is bounded in V and we know by Rellich's theorem that V is compactly imbedded in H. Hence we can extract a subsequence of (v k ), also denoted (v k ), such that (v k ) converges weakly in V , strongly in H, and a.e. in Ω to some function u ⋆ ∈ V . Using Fatou's lemma and the lower semi-continuity of the V -norm, we obtain

G n (u * ) ≤ lim inf k G n (v k ) = m = inf v∈V G n (v).
As a consequence, u ⋆ is a minimizer of G n : it satisfies the associated Euler-Lagrange equation, i.e ( d du G n (u * ), w) = 0, ∀w ∈ V , which reads

(u * , w) + 1 11 (-18u n + 9u n-1 -2u n-2 , w) + 6 11 ∆t [(∇u * , ∇w) + (g(u * ), w)] = 0.
A possible value of u n+1 is then u * .

Proposition 4.2 (Uniqueness). Assume that λ 1 + 11 6∆t > c g . Then for every choice of (u n-2 , u n-1 , u n ) ∈ V 3 , there exists at most one u n+1 ∈ V which satifies (45).

Proof. Let u n-2 , u n-1 and u n be fixed in V . Suppose that u n+1 and v n+1 are two solutions of the scheme. We have

     u n+1 - 18 11 u n + 9 11 u n-1 - 2 11 u n-2 - 6 11 ∆t ∆u n+1 + 6 11 ∆t g(u n+1 ) = 0, v n+1 - 18 11 u n + 9 11 u n-1 - 2 11 u n-2 - 6 11 ∆t ∆v n+1 + 6 11 ∆t g(v n+1 ) = 0. Denote δu = u n+1 -v n+1 . We have δu - 6 11 ∆t ∆(δu) + 6 11 ∆t g(u n+1 ) -g(v n+1 ) = 0
On multiplying by δu in H, we get

|δu| 2 0 + 6 11 ∆t |δu| 2 1 + 6 11 ∆t g(u n+1 ) -g(v n+1 ), δu = 0
Using the one-sided Lipschitz condition (42), we obtain

g(u n+1 ) -g(v n+1 ), δu ≥ -c g |δu| 2 0 .
Hence,

|δu| 2 0 + 6 11 ∆t |δu| 2 1 ≤ 6 11 ∆tc g |δu| 2 0 ,
and by the Poincaré inequality (40), we have

λ 1 + 11 6∆t -c g |δu| 2 0 ≤ 0.
The proof is complete.

Remark 4.3. Notice that if c g ∆t ≤ 11/6, then uniqueness is ensured.

In order to achieve stability, the idea is to exhibit a Lyapunov functional for (45) in a way such that the set of equilibria set associated to the functional is related to the set of equilibria of the Allen-Cahn equation.

First recall the energy functional associated to the Allen-Cahn equation (cf. ( 2)),

E : V -→ R u -→ E(u) = 1 2 |u| 2 1 + (G(u), 1).
The growth assumption (37) on g and Sobolev imbeddings imply that E is well-defined and of class C 1 on V [START_REF] Kavian | Introduction à la théorie des points critiques et applications aux problèmes elliptiques[END_REF].

The set of critical points of E is defined as

S = {u ∈ V, -∆u + g(u) = 0 in V ′ }.
The Lyapunov functional E : V 3 -→ R associated to the BDF3 is defined by

E(u, v, w) = E(u) + 5 12∆t |v| 2 0 + 1 6∆t |w| 2 0 . (48)
Remark that the set of critical points of E is S × {0} × {0}, which, in a sense, is equivalent to S.

The following proposition shows the nonincreasing property. As a shortcut, we denote

E n = E(u n , ∂u n , ∂ 2 u n ) for all n ≥ 2.
Proposition 4.4 (Lyapunov stability). Assume that 5 6∆t ≥ c g and let (u n ) n≥2 be a sequence in V generated by (45). Then for all n ≥ 2,

E n+1 + 1 2 ∂u n+1 2 1 + ( 5 6∆t -c g ) ∂u n+1 2 0 + 1 4 ∆t ∂ 2 u n+1 2 0 + 1 6 ∆t ∂u n+1 -∂u n + ∂u n-1 2 0 ≤ E n . (49)
Proof. We write the BDF3 scheme (45) in the following form, valid for all n ≥ 2 (cf. (44)),

1 ∆t 3 j=1 1 j ∂ j u n+1 -∆u n+1 + g(u n+1 ) = 0. ( 50
)
On multiplying (50) by ∂u n+1 in H, we obtain

1 ∆t 3 j=1 1 j (∂ j u n+1 , ∂u n+1 ) +(∇u n+1 , ∇(∂u n+1 )) + (g(u n+1 ), ∂u n+1 ) = 0. Σ (51)
By the Taylor-Lagrange theorem and (38), we deduce that

∀a, b ∈ R G(b) -G(a) ≥ (b -a)g(a) - c g 2 (a -b) 2 ,
and so

1 ∆t Σ + (∇u n+1 , ∇(∂u n+1 )) ≤ (G(u n ), 1) -(G(u n+1 ), 1) + c g ∂u n+1 2 0 .
Using the well known identity

(A, A -B) = 1 2 (A, A) - 1 2 (B, B) + 1 2 (A -B, A -B), (52)
we obtain

1 ∆t Σ + 1 2 u n+1 2 1 -|u n | 2 1 + 1 2 ∂u n+1 2 1 ≤ (G(u n ), 1) -(G(u n+1 ), 1) + c g ∂u n+1 2 0 .
By mimicking identity [START_REF] Haraux | Systèmes dynamiques dissipatifs et applications[END_REF], we check that Σ can be written

Σ = 5 6 ∂u n+1 2 0 + 1 4 ∂ 2 u n+1 2 0 + + 1 6 ∂u n+1 -∂u n + ∂u n-1 2 0 + 5 12 ∂u n+1 2 0 -|∂u n | 2 0 + 1 6 ∂ 2 u n+1 2 0 -∂ 2 u n 2 0 . (53)
This is the key to the Lyapunov stability. Indeed, we find

E n+1 + 1 2 ∂u n+1 2 1 + 5 6∆t ∂u n+1 2 0 + 1 4 ∆t ∂ 2 u n+1 2 0 + 1 6 ∆t ∂u n+1 -∂u n + ∂u n-1 2 0 ≤ E n + c g ∂u n+1 2 0 ,
and the proposition is proved.

Convergence to equilibrium for the BDF3 scheme

Let (u n ) n be a sequence in V . We define its ω-limit set by

ω((u n ) n ) = {u ⋆ ∈ V, ∃n k -→ ∞ s.t. u nk -→ u ⋆ in V } .
Throughout this section, we assume that the stability condition holds, that is 5 6∆t ≥ c g . In particular, the uniqueness condition is also valid. Proposition 5.1. Let (u n ) n≥2 be a sequence in V generated by (45).

Then (u n ) n≥2 is bounded in V , ∂u n -→ 0 in V , ∂ j u n -→ 0 in H for j = 2, 3 and E is constant on ω((u n ) n ).
Proof. We deduce from (46) that

E(u, v, w) ≥ 1 2 (1 - k 1 λ 1 ) |u| 2 1 + 5 12∆t |v| 2 0 + 1 6∆t |w| 2 0 -k 2 |Ω|, (54)
for all (u, v, w) ∈ V × H × H. In particular, E n is bounded from below. Adding the fact that (E n ) n is nonincreasing by the stability inequality (49), we conclude that E n converges in R to some real number E ⋆ . From the boundedness of (E n ) and estimate (54), we deduce that (u n ) n≥2 is bounded in V . Using now the stability inequality (49), we have

1 2 ∂u n+1 2 1 + 1 4∆t ∂ 2 u n+1 2 0 + 1 6 ∆t ∂u n+1 -∂u n + ∂u n-1 2 0 ≤ E n -E n+1 , (55)
for all n ≥ 2. By summing on n, we find that

∞ n=2 ∂u n+1 2 1 + ∂ 2 u n+1 2 0 ≤ C ∞ n=2 E n -E n+1 < C E 2 -E ⋆ ,
where C = (min{1/2, 1/4∆t}) -1 . Since E 2 = E(u 2 , ∂u 2 , ∂ 2 u 2 ) < ∞, this implies that ∂u n -→ 0 in V and ∂ 2 u n -→ 0 in H. Thus, E(u n ) tends to the same limit E ⋆ as E n . This implies that E is constant and equal to E ⋆ on ω((u n ) n ).

In a similar way, we obtain from (55) that

∞ n=2 ∂u n+1 -∂u n + ∂u n-1 2 0 ≤ 6∆t(E 2 -E ⋆ ) < +∞.
Thus, ∂u n+1 -∂u n + ∂u n-1 → 0 in H and since

∂ 3 u n+1 = ∂u n+1 -2∂u n + ∂u n-1 = ∂u n+1 -∂u n + ∂u n-1 -∂u n ,
we see that ∂ 3 u n+1 → 0 in H. The proof is achieved.

Proposition 5.2. Let (u n ) n≥2 be a sequence in V generated by (45). Then the set ω((u n ) n ) is a nonempty compact and connected subset of V included in S.

Proof. We first claim that the set {u n : n ∈ N} is precompact in V . The idea is to use elliptic regularity in order to find a constant q > 1 such that (u n+1 ) is bounded in W 2,q (Ω) with W 2,q (Ω) compactly imbedded in V . Here, W 2,q (Ω) is the standard Sobolev space based on the L q (Ω) space. It is well-known [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF] that

W 2,q (Ω) is compactly imbedded in V if 2 - d q > 1 - d 2 ⇐⇒ q > 2d 2 + d . ( 56
)
From the BDF3 scheme (45), we have

∆u n+1 = l(u n+1 , u n , u n-1 , u n-2 ) + g(u n+1 ), (57)
where the functional l depends linearly on its variable. Thus, l(u n+1 , u n , u n-1 , u n-2 ) is bounded in L 2 (Ω) (and also in L q (Ω), for any 1 < q ≤ 2). By Proposition 5.1, (u n ) n is bounded in V. This implies that (g(u n+1 )) n is bounded in L q (Ω) for an appropriate choice of q. Now we discuss according to d.

First assume that d ≥ 3. For any q ∈ (1, 2], we have, using the growth condition (37),

||g(u n+1 )|| L q (Ω) ≤ C 1 + ||u n+1 || L (p+1)q (Ω) ,
for some constant C independent of n. By Sobolev imbeddings, (u n+1 ) is bounded in L 2 * (Ω) where 2 * = 2d d-2 . Since p(d -2) < 4, we have p + 1 < d+2 d-2 , and by choosing q = min{ 2 * p+1 , 2} ∈ (1, 2], relation (56) is satisfied and g(u n+1 ) is bounded in L q (Ω). By elliptic regularity [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF], from (57) we deduce that (u n+1 ) is bounded in W 2,q (Ω) and this concludes the case d ≥ 3.

If d = 1 or d = 2, we obtain from Sobolev embeddings that (g(u n+1 )) n is bounded in any L r (Ω), r < ∞, so we may choose q = 2 in (57) and we conclude as previously. This proves the claim concerning precompactness of the set ω((u n ) n ).

It is well known that ω((u n ) n ) is closed in V , so it is a nonempty compact subset of V . Using that ∂u n+1 1 = u n+1 -u n 1 converges to 0 (cf. Proposition 5.1), a standard argument shows that ω((u n ) n ) is also connected in V .

Let now u ⋆ ∈ ω((u n ) n ). There exists n k -→ ∞ such that u nk -→ u ⋆ strongly in V and (u nk ) complies with (50). By Proposition 5.1, ∂ j u nk converges to 0 for j = 1, 2 and 3. Hence, -∆u ⋆ + g(u ⋆ ) = 0, so u ⋆ belongs to S as claimed.

The results obtained so far are valid if g ∈ C 1 (R, R) only. The following lemma, which is the Lojasiewicz-Simon inequality, uses that g is real analytic. Lemma 5.3. Let u ⋆ ∈ S. Then there exist constants θ ∈ (0, 1 2 ) and δ > 0 depending on u * such that for any u ∈ V satisfying |u -u ⋆ | 1 < δ, we have

|E(u) -E(u ⋆ )| 1-θ ≤ -∆u + g(u) V ′ . (58)
Proof. The inequality has been proved in [START_REF] Haraux | The convergence problem for dissipative autonomous systems[END_REF]Proposition 11.4.1].

Theorem 5.4. Let (u n ) n≥2 be a sequence in V generated by (45). Then the whole sequence converges to a steady state u ∞ in V , with u ∞ ∈ S.

Proof. The union of balls {B δ (u ⋆ ), u ⋆ ∈ ω((u n ) n )} forms an open covering of ω((u n ) n ). Due to the compactness of ω((u n ) n ) in V , we can extract a finite subcovering {B δi (u ⋆ i )} i=1,...,m in a way that the constants θ i and δ i corresponding to u * i in Lemma 5.3 are indexed by i. From the definition of ω((u n ) n ), we can always find a sufficiently large n 0 such that u n ∈ we have for all n ≥ n 0 ,

|E(u) -E ⋆ | 1-θ ≤ -∆u + g(u) V ′ ,
where E ⋆ is the constant value of E on ω((u n ) n ). We will also choose n 0 large enough so that ∂ j u n+1 0 ≤ 1 for j = 1, 2 and for all n ≥ n 0 . Denote φ n = E n -E ⋆ so that (φ n ) n is nonincreasing and bounded from below by 0.

Let n ≥ n 0 . Using the well-known inequality (a + b) 1-θ ≤ a 1-θ + b 1-θ satisfied by all nonnegative numbers a and b, we find that

(φ n+1 ) 1-θ ≤ |E(u n+1 ) -E ∞ | 1-θ + 5 12∆t ∂u n+1 2 0 + 1 6∆t ∂ 2 u n+1 2 0 1-θ ≤ -∆u n+1 + g(u n+1 ) V ′ + 5 12∆t ∂u n+1 2 0 + 1 6∆t ∂ 2 u n+1 2 0 1-θ ≤ 1 ∆t 3 j=1 ∂ j u n+1 j V ′ + 5 12∆t ∂u n+1 2 0 + 1 6∆t ∂ 2 u n+1 2 0 1-θ .
The last inequality is obtained by the version (50) of the BDF3 scheme. Next, we note that

3 j=1 ∂ j u n+1 j = ∂u n+1 + 1 2 ∂u n+1 -∂u n + 1 3 ∂u n+1 -2∂u n + ∂u n-1 = 3 2 ∂u n+1 - 5 6 ∂u n + 1 3 ∂u n+1 -∂u n + ∂u n-1 = 2 3 ∂u n+1 + 5 6 ∂u n+1 -∂u n + 1 3 ∂u n+1 -∂u n + ∂u n-1 = 2 3 ∂u n+1 + 5 6 ∂ 2 u n+1 + 1 3 ∂u n+1 -∂u n + ∂u n-1 .
Since H is continuously imbedded into V ′ (cf. (41)), this yields

(φ n+1 ) 1-θ ≤ 2λ -1/2 1 3∆t ∂u n+1 0 + 5λ -1/2 1 6∆t ∂ 2 u n+1 0 + λ -1/2 1 3∆t ∂u n+1 -∂u n + ∂u n-1 0 +( 5 12∆t ) 1-θ ∂u n+1 2(1-θ) 0 + ( 1 6∆t ) 1-θ ∂ 2 u n+1 2(1-θ) 0 ≤ c ∂u n+1 2 1 + ∂ 2 u n+1 2 0 + ∂u n+1 -∂u n + ∂u n-1 2 0 1 2 , ( 59
)
where c is a constant independent of n. Assume first that φ n+1 > φ n /2. Then we have

(φ n ) θ -(φ n+1 ) θ = φ n φ n+1 θ x θ-1 dx ≥ θ(φ n ) θ-1 φ n -φ n+1 ≥ θ2 θ-1 (φ n+1 ) θ-1 φ n -φ n+1 .
From Proposition 4.4 and 59, we obtain

(φ n ) θ -(φ n+1 ) θ ≥ C ∂u n+1 2 1 + ∂ 2 u n+1 2 0 + ∂u n+1 -∂u n + ∂u n-1 2 0 1 2 . ( 60
)
If φ n+1 ≤ φ n /2, by using Proposition 4.4 again, we find that for all n ≥ 2,

(φ n ) 1 2 -(φ n+1 ) 1 2 ≥ √ 2 -1 √ 2 (φ n ) 1 2 ≥ √ 2 -1 √ 2 (φ n -φ n+1 ) 1 2 ≥ C ′ ∂u n+1 2 1 + ∂ 2 u n+1 2 0 + ∂u n+1 -∂u n + ∂u n-1 2 0 1 2 . (61)
In both cases, inequalities (60) and (61) imply that for all n ≥ n 0 ,

∂u n+1 1 ≤ 1 C (φ n ) θ -(φ n+1 ) θ + 1 C ′ (φ n ) 1 2 -(φ n+1 ) 1 2
.

By summing on n, we find that

∞ n=n0 u n+1 -u n 1 ≤ 1 C (φ n0 ) θ + 1 C ′ (φ n0 ) 1 2 < ∞.
This implies that (u n ) n converges to some u ∞ in V . By Proposition 5.2, u ∞ belongs to S.

Remark 5.5. By using the gradient stability proved in Section 3 and by arguing as for the BDF3 scheme, we could obtain a result similar to Theorem 5.4 for the BDF4 and BDF5 methods.

Numerical simulations

We present some numerical simulations to illustrate the theoretical results. We consider the 1d Allen-Cahn equation endowed with homogeneous Dirichlet boundary conditions,

     u t -∆u + g(u) = 0, (x, t) ∈ (0, L) × (0, T ], u(0, t) = u(L, t) = 0, t ∈ (0, T ], u(x, 0) = u 0 (x), x ∈ (0, L), (62)
where g(s) = s 3 -s. Note that c g = 1 (cf. ( 38)).

For the space discretization, we use P 1 continuous finite elements on a regular subdivision of [0, L] with mesh size h = 1/(M + 1). Due to the Dirichlet boundary conditions, the dimension of the finite element space V h is M . By using an appropriate basis (namely, an orthonormal basis of V h for the L 2 (0, L) scalar product), it is easily seen that this space discrete version of (62) is a gradient flow in R M of the form (3), for a function F which satisfies (4)- [START_REF] Attouch | On the convergence of the proximal algorithm for nonsmooth functions involving analytic features[END_REF] (see e.g. [START_REF] Alaa | Convergence to equilibrium for discretized gradient-like systems with analytic features[END_REF]). Moreover, F is a polynomial of M variables (and therefore real analytic), so that convergence to equilibrium holds for this gradient flow.

For the time discretization, we use the BDF3 scheme. By gradient stability [START_REF] Alaa | Convergence to equilibrium for discretized gradient-like systems with analytic features[END_REF][START_REF] Stuart | Dynamical systems and numerical analysis[END_REF], convergence to equilibrium also holds for this fully discrete version of (62).

Regarding the numerical resolution, we use the standard nodal basis of V h ("hat functions"), and in the following, we identify V h to R M through this basis. We obtain the scheme

B h 11 U n+3 -18 U n+2 + 9 U n+1 -2 U n 6∆t + A h U n+3 + G h (U n+3 ) = 0, (63)
where A h and B h are M ×M symmetric positive-definite matrix, G h : R M → R M corresponds to the nonlinearity g, and U n ∈ R M is the unknown vector (see e.g. [START_REF] Elliott | The global dynamics of discrete semilinear parabolic equations[END_REF]). The matrix A h is a discretization of -∆ and B h is the mass matrix. This scheme is fully implicit. When U n , U n+1 and U n+2 are known, U n+3 can be obtained as a solution of the nonlinear problem,

H n (V ) = 0, ( 64
)
where H n : R M -→ R M is the function,

H n (V ) = B h 11V -18U n+2 + 9U n+1 -2U n 6∆t + A h V + G h (V ).
For the numerical resolution of (64), one possibility is to use a Newton method. This reads: let V 0 = U n+2 and for k = 0, 1, 2, . . . , until convergence, let V k solve the linearized equation

H n (V k ) + dH n (V k )(V k+1 -V k ) = 0, ( 65
)
where dH n (V ) = B h /∆t + A h + dG h (V ) is a M × M symmetric matrix. Yet, the computation of dG(V k ) is costly, so we actually use an approximate Newton method, i.e. an approximation of dG(V k ).

B(U n+1 -U n )/∆t + A(U n+1 + U n ))/2 + [G(U n+1 ) + G(U n )]/2 = 0.
The choice of this second order scheme for U 1 and U 2 guarantees the third order accuracy of the BDF3 scheme [START_REF] Stuart | Dynamical systems and numerical analysis[END_REF].

We see that the discrete solution rapidly converges to a steady state, with values close to ±1, except for three transition layers (one in the middle, and two near the boundaries 0 and L. The equilibrium is reached graphically after only 14 iterations (i.e. with an accuracy of order 0.01 in the max norm of U n ).

Another simulation that shows the rapidity of convergence to equilibrium is given in Figure 2. It is a plot of the pseudo-energies E n k versus time iteration for the BDF methods of order k = 1, 2 and 3. The pseudo-energy here is the Lyapunov function associated to the scheme [START_REF] Stuart | Dynamical systems and numerical analysis[END_REF]. In particular, E n 3 is given by E(u n h ) where u n h ∈ V h corresponds to U n . Since ∆t is small enough, the pseudo-energies are decreasing, as expected. We note that they all converge to the same value, which is the energy level of the discrete steady state. In order to check numerically the third order accuracy of the BDF3 scheme, we have computed in Table 2 the error for the time discretization,

err ∆t = max 0≤n≤N u ∆t h (t n ) -u ∆tsol h (t n ) L 2 (0,L)
.

The maximum is evaluated on the coarse grid t n = n∆t, n = 0, 1, . . . , N , with ∆t = T /N . The notation u ∆t h (t n ) represents the function in V h corresponding to U n at time t n . Since the exact solution is not known, we used instead the discrete solution u ∆tsol h on a fine grid with total time solution T = 2 and time stepsize ∆t sol = 1/512 corresponding to 1024 iterations. The values of M = 200 and u 0 were chosen as above. We see that the ratio err ∆t /err ∆t 2 of consecutive errors is close to 8 = 2 3 , and this is consistent with the third order accuracy of the BDF3 scheme.

  m i=1 B δi (u ⋆ i ) for all n ≥ n 0 . Consider θ = min i=1,...,m {θ i }, then by Lemma 5.3 and Proposition 5.1,

Figure 1 .

 1 Figure 1. Evolution of Un to equilibrium from a regular initial condition.

Figure 1

 1 Figure 1 represents the evolution of the sequence (U n ) generated by (63) from a regular initial datum u 0 (x) = 1 5 sin( 2πx L ). The spatial discretization was chosen small enough to ensure good accuracy of the finite elements method by setting L = 20, ∆x = 0.0995 corresponding to M = 200 space nodes. The total time of the simulation is T = 10 with a time stepsize of ∆t = 0.5 ≤ 5 6 cg corresponding to n = 20 iterations. The first two iterations U 1 and U 2 are approximated using the Crank-Nicolson scheme, i.e for n = 0, 1

Figure 2 .

 2 Figure 2. Peudo-energy vs time iteration for BDF1, BDF2 and BDF3.

  An easy way to choose the coefficients of A r4 is to take the elements of a given diagonal all equal. This gives A r4 , and in turn, A q4 and A q4 . Using Sylvester's criterion[START_REF] Kavian | Introduction à la théorie des points critiques et applications aux problèmes elliptiques[END_REF], we easily check that A q4 and A r4 , given by[START_REF] Shen | Numerical approximations of Allen-Cahn and Cahn-Hilliard equations[END_REF], are positive definite. Indeed, the principal minors of A

	/12 -23/24 13/24 -1/8 -23/24 0 0 0 13/24 0 0 0	   
	-1/8	0	0	0

q4 are 225/144, 25 286/144 2 and 1 522 296/144 3 , whereas the principal minors of A r4 are 75/144, 3 509/144 2 , 155 448/144 3 and 6 460 912/144 4 . Proposition 2.8. The 5-step BDF method is quadratic-stable. Namely,

Table 1 .

 1 Uniqueness and gradient stability conditions for BDF methods

	k	1	2	3	4	5	6	≥ 7
	α k	1	3/2	11/6	25/12	137/60	49/20	k j=1 1/j
	β k	2 or ∞	2	5/3	41/72	1/10	?	0 (not g. s.)

Table 2 .

 2 Convergence error of the

	BDF3 scheme N ∆t	err ∆t	ratio
	16	1/8	1.8483e-03	-
	32	1/16	2.4165e-04	7.6487
	64	1/32	3.0741e-05	7.861
	128	1/64	3.8251e-06	8.0367
	256	1/256	4.2679e-07	8.9624

Acknowledgement

The authors are thankful to Frédéric Bosio for his help on the proof of Proposition 2.7.


Saman Noorzadeh

Pierre Maurel

Thomas Oberlin

Rémi Gribonval

Christian Barillot

Multi-modal EEG and fMRI Source Estimation Using Sparse Constraints

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The integration of different modalities has brought up interest to understand the complex procedure of the brain. Two of these modalities are functional Magnetic Resonance Imaging (fMRI) and electroencephalography (EEG). fMRI measures the oxygenation of the blood flow, which is closely correlated to the neuronal activity. The oxygenated blood flows to the active neurons, and this leads to contrast images called Blood Oxygenation Level Dependent (BOLD). Although fMRI has a high spatial resolution, the images taken from the volume of the brain have limited temporal precision. On the other hand, EEG measures the neuronal electrical activity of the brain in the scale of milliseconds. However, it does not directly provide the spatial location of active neurons, and the localization of the active electric sources from EEG is only possible by solving an inverse problem which is ill-posed and undetermined [START_REF] Vogel | Computational methods for inverse problems[END_REF]. A review of the techniques for solving the EEG inverse problem can be found in [START_REF] Becker | A performance study of various brain source imaging approaches[END_REF].

The use of both EEG and fMRI in one multi-modal framework can result in an estimation of neuronal activities with a high resolution in time and also in space. The complementary strengths that these two modalities can provide, has stimulated a variety of studies [START_REF] Rosa | EEG-fMRI integration: a critical review of biophysical modeling and data analysis approaches[END_REF][START_REF] Karahan | Tensor analysis and fusion of multimodal brain images[END_REF]. These methods could be fMRI-constrained, or EEG-informed where models driven on one modality act as the spatial (or temporal) prior on the solutions of the other modality [START_REF] Brookings | Using ICA and realistic bold models to obtain joint EEG/fMRI solutions to the problem of source localization[END_REF][START_REF] Bagshaw | Analysis of the EEG-fMRI response to prolonged bursts of interictal epileptiform activity[END_REF]. EEG and fMRI can also be used in symmetrical approaches, so that the information fusion is based on both data at the same time. For example, joint Independent Component Analysis (ICA) is used to analyze both data in a joint space as in [START_REF] Moosmann | Joint independent component analysis for simultaneous EEG-fMRI: principle and simulation[END_REF], or Bayesian approaches are used in this regard [START_REF] Deneux | EEG-fMRI fusion of paradigm-free activity using kalman filtering[END_REF]. According to the literature, the relation of neuronal activity, hemodynamics, and fMRI is still unclear [START_REF] Logothetis | Neurophysiological investigation of the basis of the fMRI signal[END_REF]; however, the existing methods model this relationship based on the balloon model or neural mass models [START_REF] Babajani | Integrated MEG/EEG and fMRI model based on neural masses[END_REF]. Although the integration of these two modalities has been studied a lot, it is mainly validated according to simulated data. A symmetrical approach is proposed in [START_REF] Oberlin | Symmetrical EEG-fMRI imaging by sparse regularization[END_REF] according to the balloon model, and the current paper adapts the proposed linear approaches to validate it on real data.

In the current study we propose to adapt a symmetrical approach based on EEG and fMRI for the estimation of the brain sources and validate it on real data. The joint algorithm is composed of two main parts: the simple linear model of both data, and the sparsity of the sources. The linear model is a first approximation and in this study we intend to show that despite the possible limits, the sparse and joint model can provide better spatio-temporal results in source reconstruction for real data, so that it can justify developing more complex models based on multi-modality.

In section 2 EEG and BOLD inverse problems are explained, based on which the joint model is introduced with the model parameters. The sources are then estimated with the proximal algorithm. Section 3 explains the experiment, the multi-modal data, and its analysis, and the results based on the real data, shown in 4, validates the efficiency of the joint approach, and justifies the calibration of the model parameters. Finally the conclusions are discussed in section 5.

Method

Model EEG Inverse Problem

The aim is to estimate the neural sources at time instants of EEG. According to the forward problem (see section 3), the J × L leadfield matrix, G, is calculated, where J and L are the number of electrodes and sources, respectively. This matrix contains coefficients that linearly link the source amplitudes to the electrode potentials. The estimation of the source matrix, S, from the following equation, is the solution of the EEG inverse problem:

E = GS + N e , (1)
where E is a J ×T e matrix corresponding to the signal obtained from J electrodes during T e instants. As noted, there is no unique solution of S because EEG has much fewer channels than the number of sources (L J). That is why the hypothesis of the sparsity is made about the sources.

BOLD Inverse Problem

In the framework of balloon model [START_REF] Friston | Nonlinear responses in fMRI: the balloon model, volterra kernels, and other hemodynamics[END_REF], the BOLD response is written as the linear convolution of the neural activity and the hemodynamic response, h(t). This is depicted in Fig. 1. Here, the electrical source locations are corresponding to voxels. So, considering B a L × T b matrix as the observed BOLD signal of T b instances for L voxels, the following equation holds:

B = SH + N b . (2)
where N b is the Gaussian noise, and H is a T e ×T b matrix whose (p, q) th element is h(t p -u q ), where t p and u q correspond to the time samples of EEG and BOLD, respectively. Note that T e T b makes the BOLD inverse problem undetermined in time, that is why the sparsity of sources is important in finding the solution.

Joint Inverse Problem Two minimization problems can be induced from the inverse problems of equations ( 1) and ( 2) to estimate sources S, and these two can be combined symmetrically in a single model. Assuming a Gaussian noise and sparse sources we get a least square problem with L 1 regularization. The sparsity is a logical hypothesis since a special mental task would activate only a certain number of sources in functional brain regions. Considering the penalty term φ(S) = ||S|| 1 = i k |S i,k |, the sources can be finally estimated as:

Ŝ = argmin S α||GS -E|| 2 2 + (1 -α)||SH -B|| 2 2 + λ j ||S|| 1 , (3)
Where ||.|| 2 denotes the L 2 norm, 0 < λ j , and 0 ≤ α ≤ 1. λ j adapts the solution as a good tradeoff between the joint data-fidelity term and the penalty term. α balances the relative importance of each modality in the joint approach.

Proximal Algorithm

Proximal algorithms are for obtaining optimization solutions for composite functions, which might not be differentiable. The joint problem has already been formulated as

Ŝ = argmin S f (S) + λ j φ(S) , (4)
where f is the data-fidelity term, a proper convex function whose gradient is Lipschitz continuous, while the defined φ is not differentiable. First consider the general definition of the proximity operator prox λ,φ (y) = argmin y { 1 2 ||x -y|| 2 2 + λφ(y)}. Now, an iterative method, called Iterative Soft Thresholding (IST) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], is used with a small gradient descent step µ to find the solution as follows:

S (k+1) = prox µ,φ S (k) -µ∇f (S (k) ) . (5)
It is proven that this algorithm converges to a minimum [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF].

Other norms (i.e. L 2 or mixed L 12 ) could be used for the φ regularization to avoid the estimation of scattered sources [START_REF] Oberlin | Symmetrical EEG-fMRI imaging by sparse regularization[END_REF]; but, comparing different penalty terms is beyond the scope of this paper. Interested readers can refer to [START_REF] Oberlin | Symmetrical EEG-fMRI imaging by sparse regularization[END_REF] for a theoretical investigation of the differences between those terms in that context. Having the L 1 norm penalty, the proximity operator on the elements of the source matrix, S ∈ R V ×T is given by the following equation:

[prox λ,||.||1 (S)] v,t = S v,t |S v,t | (|S v,t | -λ) (6)

Parameter Calibration

Consider again the joint model defined in equation [START_REF] Rosa | EEG-fMRI integration: a critical review of biophysical modeling and data analysis approaches[END_REF]. Indeed, α is hard to be tuned between the two sides of the data-fidelity term, since these two data-fidelity terms have no reason to have the same order of magnitude; so, a calibration of the respective roles of the two signals should be embedded. Therefore, λ e and λ b are separately defined for each data-fidelity term (corresponding to EEG and BOLD, respectively). The joint model can thus be re-written considering the new calibration parameters:

Ŝ = argmin S α λ e λ j ||GS -E|| 2 2 + (1 -α) λ b λ j ||SH -B|| 2 2 + φ(S) , α ∈ [0, 1]. (7)
Note that if α is 0 (or 1), the model will be uni-modal based on BOLD (or EEG). So, if α ∈ (0, 1) (ideally α = 0.5) gives better estimations than uni-modal approaches, this can validate the joint method.

Experimental Data and Implementations

The simultaneous EEG and fMRI are recorded from 8 healthy right-handed subjects while doing a motor task of clenching the right hand repeatedly. A study session consists of 10 blocks of 20s from which 5 are task each followed by a rest block. EEG is obtained via 64 electrodes based on the standard 10-20 system with f s = 200 Hz. The ECG (Electrocardiogram) and the gradient artifacts caused by MR system are separately eliminated. The sources are estimated for the first 15 seconds of a task block, since at the end of each block the subjects might not guarantee to do the task well because of fatigue. At the same time, the subject would enter a 3T MRI system, for the acquisition of the structural and functional MRI. Every 2 seconds a volume of the brain is being recorded. The fMRI data is then spatially corrected for the head motion. The head model, is made up of 3 nested concentric spheres with homogeneous and isotropic conductivity representing the brain, the skull, and the scalp [START_REF] Berg | A fast method for forward computation of multiple-shell spherical head models[END_REF]. The electrodes are placed on the scalp sphere surface defined in the MRI anatomical space (according to the 10-20 position). The sources are modeled by fixed dipoles located at the voxels' locations that are selected by brain extraction and segmentation (to consider the gray matters as the possible source locations). These sources are then down-sampled so that their number correspond to the physiological number of the current sources. A sphere, representing the brain, is then fitted to the sources. Finally, after defining another sphere for the skull, the subject-specific 3-layer spherical model is defined with different conductivities considered for each of the layers according to [START_REF] Berg | A fast method for forward computation of multiple-shell spherical head models[END_REF] and [START_REF] Albera | Brain source localization using a fourth-order deflation scheme[END_REF]. The forward problem is then solved to linearly correlate the sources and the observed EEG potentials in the lead-field matrix.

Results

The subjects were involved in the task of right hand movement, so we expect the activation on the motor cortex of the left hemisphere. The motor cortex is composed of the primary motor cortex (M1), premotor cortex (called PMA), the supplementary motor area (SMA), and the posterior parietal cortex (PPC). The brain space is aligned with the AAL (Automated Anatomical Labeling) atlas, based on which the brain region in which each source is placed is determined. When the sources S are estimated, the power p(v) of each source v is calculated as the L 1 norm of its activity during all task duration: p(v) = ||S(v, :)||. So, the power score of each region r is defined as

power(r) = v∈r p(v). (8)
Algorithm 1 Uni-modal parameter optimization function

1: function unimodal(α) 2: λj ← 1 3:
for λ ∈ [λ min , λ max ] do 4:

for each subject i do 5:

Compute S from (7) 6:

M (i) = Q(S) 7:
end for 8:

Q(λ) ← mean(M ) 9:

end for 10: for α ∈ (0, 1) do 5:

λ * ← argmax λ Q(λ
for each subject i do 6:

Compute S from (7) 7:

M (i) = Q(S) 8:
end for 9:

Q(α, λj) ← mean(M ) 10:

end for 11: end for 12: (α * , λ * j ) ← argmax

(α,λ j ) Q(α, λj)
Parameter Estimation The performance, noted Q, of any estimated source matrix is defined to be higher for an estimation with a power score higher in the left motor cortex and less in other regions:

Q(S) = r={M 1,P M A,P P C} lef t Power(r) - r ={M 1,P M A,P P C} lef t Power(r) (9)
So, the estimation and calibration of parameters is done based on the performance. λ b and λ e are first estimated in parallel by separate uni-modal source estimations when α is set to 0 and 1, respectively (with λ j = 1); afterwards, the joint source estimation tunes both α and λ j . This optimization algorithm is shown in algorithms 1 and 2 for better clarification. The performances of estimations for different parameters are then illustrated in Fig. 2. The detected values of λ * e and λ * b respectively equal to 3.4 × 10 -5 and 10 -3 , show the different order of magnitude of the two data-fidelity terms, and confirms our claim on the need of calibration parameters. Moreover, the value of α * = 0.6 maximizes the performance which shows that the estimation is almost equally balanced between the two EEG and BOLD data, and that the uni-modal parameters are well calibrated. Therefore, α can always be set to 0.5 and λ is the only parameter that remains to be estimated in such method. This calibration process can thus reduce the number of parameters and simplifies the procedure.

Discussions Here, some results of the joint method are depicted based on the estimated parameters. Fig. 3 shows the detected power score, defined in equation ( 8), of the estimation in the motor cortex. Since the task brings more attention toward the activation of PMA and M1, the activities of the left and right PMA and M1 for all subjects are compared. A higher difference between the power of the left and right regions is seen in the joint approach (p = 0.02 for a righttailed paired t-test), while this difference is not significant in the uni-modal ones (p = 0.1 for BOLD-only and p = 0.4 for EEG-only).

The spatial result of a subject is illustrated in Fig. 4. In fMRI-based method, increased activity is seen in the left motor cortex. The activities in unexpected regions can be either due to vascular activities or the errors of the BOLD linear model [START_REF] Monti | Statistical analysis of fMRI time-series: a critical review of the GLM approach[END_REF]. However, a precise localization of the active areas (SMA and M1) can be seen; while EEG detected activities in SMA but not the left motor cortex. The unexpected activities can be due the fact that every instantaneous neuronal activity is recorded in EEG because of its high temporal precision. The multi-modal approach is capable of avoiding the false positives of the uni-modal approaches. Besides, it has detected high activities in both SMA and the left motor cortex.

Finally, the temporal signal of a detected source is shown (Fig. 5). This result cannot be quantitatively validated. Analysis can be on spectral content of EEG. Here, the temporal resolution is only compared, and further analysis is left for future studies. Here, the detected source has low temporal resolution in fMRIonly approach comparing the EEG-only one. The parts detected as zero are caused by the sparsity constraint. The temporal resolution of the estimation detected by the joint approach is as high as the temporal resolution of EEG.

Conclusions and Perspectives

This paper showed that a symmetric multi-modal approach based on EEG and fMRI, can improve the spatio-temporal resolution of the estimation of neuronal sources. Linear models of EEG and fMRI are used, and the sparsity constraint is applied on the sources as a logical hypothesis which helps finding solutions to the undetermined EEG and BOLD inverse problems. This model is proved to be successful by the results on real data of 8 subjects doing a right hand motor task. Moreover, the detected values of the model parameters validated the efficiency of the joint estimation when both data are symmetrically and almost equally involved in the model. This work supports our assumption that a simple joint linear model can work on real data, so that in future studies more complex models can be taken into consideration.

Fig. 1 :

 1 Fig.1: On the left two example of neuronal stimuli. On the right, their corresponding HR. The BOLD response is the convolution of the sources and the HRF.

) 11 :

 11 return λ * 12: end function Algorithm 2 Multi-modal parameter optimization algorithm 1: λe ← unimodal(1) 2: λ b ← unimodal(0) 3: for λj ∈ [λ min j

Fig. 2 :

 2 Fig. 2: (a) selection of λ * e with EEG only performance; (b) selection of λ * b with fMRI only performance; (c) selection of λ * j and α * with EEG and fMRI jointly.

Fig. 3 :Fig. 4 :

 34 Fig. 3: Box plots over eight subject of the power score (8) for activations of the union of the PMA and M1 regions in both hemispheres, while doing a right-hand motor task.

Fig. 5 :

 5 Fig.5: Temporal activity of a source detected in SMA by uni-vs multi-modal approach.