Maurizio Grasselli

Morgan Pierre

ENERGY STABLE AND CONVERGENT FINITE ELEMENT SCHEMES FOR THE MODIFIED PHASE FIELD CRYSTAL EQUATION

Keywords: Finite elements, second-order schemes, gradient-like systems, Lojasiewicz inequality. MSC 2010: 65M60, 65P40, 74N05, 82C26

We propose a space semi-discrete and a fully discrete finite element scheme for the modified phase field crystal equation (MPFC). The space discretization is based on a splitting method and consists in a Galerkin approximation which allows low order (piecewise linear) finite elements. The time discretization is a second-order scheme which has been introduced by Gomez and Hughes for the Cahn-Hilliard equation. The fully discrete scheme is shown to be unconditionally energy stable and uniquely solvable for small time steps, with a smallness condition independent of the space step. Using energy estimates, we prove that in both cases, the discrete solution converges to the unique energy solution of the MPFC equation as the discretization parameters tend to 0. Using a Lojasiewicz inequality, we also establish that the discrete solution tends to a stationary solution as time goes to infinity. Numerical simulations illustrate the theoretical results.

Introduction

In this paper, we analyze finite element discretizations of the modified phase field crystal (MPFC) equation

βu tt + u t = ∆[∆ 2 u + 2∆u + f (u)],
in Ω × (0, +∞), (

with periodic boundary conditions on the d-parallelepiped Ω = Π d k=1 (0, L k ) (L k > 0 for k = 1, . . . , d) in space dimension d = 1, 2 or 3. Equation (1.1) is endowed with initial conditions u(0, x) = u 0 (x), u t (0, x) = v 0 (x), x ∈ Ω.

(1.

2)

The unknown function u is the phase function and β > 0 is a relaxation parameter. The nonlinearity f is the derivative of a polynomial potential F (see (2.1)-(2.2) for a precise definition). A relevant example in applications is given by

f (s) = s 3 + (1 -ε)s, (1.3)
where ε ∈ R is constant. When β = 0, equation (1.1) is known as the phase field crystal (PFC) equation: it has been employed to model and simulate the dynamics of crystalline materials, including crystal growth in a supercooled liquid, dendritic and eutectic solidification, epitaxial growth [START_REF] Elder | Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystal[END_REF][START_REF] Elder | Modeling elasticity in crystal growth[END_REF][START_REF] Emmeric | Selected issues of phase-field crystal simulations[END_REF][START_REF] Emmerich | Phase field-crystal models for condensed matter dynamics on atomic length and diffusive time scales: an overview[END_REF][START_REF] Ohnogi | Instability of spatially periodic patterns due to a zero mode in the phase-field crystal equation[END_REF][START_REF] Provatas | Using the phase-field crystal method in the multi-scale modeling of microstructure evolution[END_REF]. In the phase field approach, the number density of atoms is approximated by the phase function u. The parameter ε in (1.3) is proportional to the undercooling i.e. ε ∼ T e -T , T e being the equilibrium temperature at which the phase transition occurs. The PFC equation is a gradient flow for the free energy

E(u) = Ω 1 2
|∆u| 2 -|∇u| 2 + F (u) dx.

(1.4)

It preserves the total mass and can be viewed as an analog of the Swift-Hohenberg equation [START_REF] Swift | Hydrodynamic fluctuations at the convective instability[END_REF].

The MPFC equation (1.1) (with β > 0) was recently proposed in [START_REF] Stefanovic | Phase-field crystals with elastic interactions[END_REF] (cf. also [START_REF] Emmerich | Phase field-crystal models for condensed matter dynamics on atomic length and diffusive time scales: an overview[END_REF][START_REF] Galenko | Phase-field-crystal and Swift-Hohenberg equations with fast dynamics[END_REF][START_REF] Galenko | Marginal stability analysis of the phase field crystal model in one spatial dimension[END_REF][START_REF] Stefanovic | Phase-field crystal study of deformation and plasticity in nanocrystalline materials[END_REF]) in order to account for elastic interactions. Equations like (1.1) have also been derived in [START_REF] Galenko | Diffuse-interface model for rapid phase transformations in nonequilibrium systems[END_REF] to take large deviations from thermodynamic equilibrium into account.The MPFC equation is no longer a gradient flow for (1.4), but it is possible to associate to a solution (u, u t ) a "pseudo-energy", obtained by adding to E a kinetic energy term (see (2.7)). This leads in a natural way to the notion of energy solution introduced by Grasselli and Wu [START_REF] Grasselli | Well-posedness and longtime behavior for the modified phase-field crystal equation[END_REF] for the MPFC equation, following a terminology used for the modified Cahn-Hilliard equation [START_REF] Grasselli | On the 3D Cahn-Hilliard equation with inertial term[END_REF][START_REF] Grasselli | On the 2D Cahn-Hilliard equation with inertial term[END_REF][START_REF] Grasselli | Trajectory and smooth attractors for Cahn-Hilliard equations with inertial term[END_REF]. Existence of a unique energy solution was proved by Grasselli and Wu in [START_REF] Grasselli | Well-posedness and longtime behavior for the modified phase-field crystal equation[END_REF] as well as the convergence of single trajectories to single stationary states. The analysis of global dynamics (that is, existence of global and exponential attractors) for the MPFC equation was carried out in [START_REF] Grasselli | Well-posedness and longtime behavior for the modified phase-field crystal equation[END_REF][START_REF] Grasselli | Robust exponential attractors for the modified phase-field crystal equation[END_REF].

From the numerical analysis point of view, the MPFC equation has been studied in [START_REF] Baskaran | Energy stable and efficient finite-difference nonlinear multigrid schemes for the modified phase field crystal equation[END_REF][START_REF] Baskaran | Convergence analysis of a second order convex splitting scheme for the modified phase field crystal equation[END_REF][START_REF] Galenko | Unconditionally stable and numerical solution of the hyperbolic phase-field crystal equation[END_REF][START_REF] Wang | Global smooth solutions of the three-dimensional modified phase field crystal equation[END_REF][START_REF] Wang | An energy stable and convergent finite-difference scheme for the modified phase field crystal equation[END_REF], while the literature on the PFC equation is more abundant (see, e.g., [START_REF] Backofen | Nucleation and growth by a phase field crystal (PFC) model[END_REF][START_REF] Cheng | An efficient algorithm for solving the phase field crystal model[END_REF][START_REF] Elsey | A simple and efficient scheme for phase field crystal simulation[END_REF][START_REF] Gomez | An unconditionally energy-stable method for the phase field crystal equation[END_REF][START_REF] Hu | Stable and efficient finite-difference nonlinearmultigrid schemes for the phase-field crystal equation[END_REF][START_REF] Wise | An energy-stable and convergent finite-difference scheme for the -field crystal equation[END_REF]). In [START_REF] Baskaran | Energy stable and efficient finite-difference nonlinear multigrid schemes for the modified phase field crystal equation[END_REF][START_REF] Baskaran | Convergence analysis of a second order convex splitting scheme for the modified phase field crystal equation[END_REF][START_REF] Wang | An energy stable and convergent finite-difference scheme for the modified phase field crystal equation[END_REF], the authors proposed unconditionally energy stable and unconditionally uniquely solvable finite difference schemes. The time discretization was based on a convex splitting of the pseudo-energy and was either first order or second order. A priori error estimates were proved assuming enough regularity on the solution. A time semi-discrete scheme was used in [START_REF] Wang | Global smooth solutions of the three-dimensional modified phase field crystal equation[END_REF] to establish the existence of a weak solution and of a unique strong solution to the MPFC equation up to any positive final time T > 0. In [START_REF] Galenko | Unconditionally stable and numerical solution of the hyperbolic phase-field crystal equation[END_REF], an unconditionally energy stable finite element scheme was derived, but no proof of convergence was given.

The main purpose of this paper is to derive and analyze a second order (in time) fully discrete finite element scheme for the MPFC equation. For the space discretization we use a splitting approach which is well known in the context of phase field models (cf., e.g., [START_REF] Elliott | The Cahn-Hilliard model for the kinetics of phase separation[END_REF][START_REF] Grasselli | A splitting method for the Cahn-Hilliard equation with inertial term[END_REF]). This argument allows to consider low order (piecewise linear) finite elements, although the analysis is carried out in a more general setting, namely a Galerkin approximation. For the time discretization, we use a modified Crank-Nicolson scheme which was introduced by Gomez and Hughes for the Cahn-Hilliard equation [START_REF] Gomez | Provably unconditionally stable, second-order time-accurate, mixed variational methods for phase-field models[END_REF]; their approach represents an interesting alternative to secant schemes (cf. the discussion in [START_REF] Wu | Stabilized second-order convex splitting schemes for Cahn-Hilliard models with application to diffuse-interface tumor-growth models[END_REF]).

We prove that our scheme is unconditionally energy stable, solvable for any time step, and uniquely solvable for small enough time steps, with a smallness assumption independent of the space step. Using the energy estimate, and assuming only some natural conditions on the initial conditions, we show that the solution of the fully discrete scheme converges to the energy solution of problem (1.1)-(1.2) as the time step τ and the space step h tend to 0. Finally, we prove that the discrete solution tends to a stationary solution as time goes to infinity. This last issue is not trivial since the set of steady states can be very complicated (see [START_REF] Pierre | Stationary solutions to phase field crystal equations[END_REF] for an analysis of the one dimensional stationary problem, cf. also [START_REF] Ella | Steady state bifurcations for phase field crystal equations with underlying two dimensional kernel[END_REF]).

For equations involving a second order time derivative such as (1.1), second order time discretizations are very interesting because they do not regularize in finite time, unlike first order schemes: a fundamental property of the continuous model is therefore reproduced at the discrete level. In contrast with the second order two-step scheme in [START_REF] Baskaran | Energy stable and efficient finite-difference nonlinear multigrid schemes for the modified phase field crystal equation[END_REF][START_REF] Baskaran | Convergence analysis of a second order convex splitting scheme for the modified phase field crystal equation[END_REF], we loose the unconditional unique solvability, but one advantage is that our one-step scheme can be used with variable time steps. Moreover, we do not assume ε < 1 in (1.3), and we do not have any restriction on the initial value of u t . Since energy solutions have a weaker regularity than the weak solutions, our convergence result as (h, τ ) → (0, 0) holds with the minimal regularity on the solution.

Our proofs are crucially based on the energy estimates. In order to establish the convergence to equilibrium, we also use the gradient-like structure of the problem and a Lojasiewicz inequality [START_REF] Lojasiewicz | Ensembles semi-analytiques[END_REF], as in the continuous case [START_REF] Grasselli | Well-posedness and longtime behavior for the modified phase-field crystal equation[END_REF]. Related results have been obtained for first order schemes applied to second-order gradient-like equations in [START_REF] Alaa | Convergence to equilibrium for discretized gradient-like systems with analytic features[END_REF][START_REF] Grasselli | Convergence to equilibrium of solutions of the backward Euler scheme for asymptotically autonomous second-order gradient-like systems[END_REF], but the case of the second-order scheme here is more involved.

The paper is structured as follows. In Section 2 we introduce the functional setting and we recall useful results concerning the continuous problem. In Section 3 we consider the space semi-discrete problem. We show its well-posedness and establish energy estimates which allow us to prove the convergence of the semi-discrete solution to the energy solution of the MPFC problem (1.1)- (1.2). This gives a framework for the fully discrete problem which is treated in Section 4. Section 5 is concerned with the convergence to equilibrium for the fully discrete problem. In Section 6, numerical simulations in one and two space dimension illustrate the theoretical results.

The continuous problem

2.1. Notation and functional spaces. For a real Banach space V with dual V ⋆ , we indicate by •, • V ⋆ ,V the duality product between V and V ⋆ . We denote by H m per , m ∈ N, the space of functions that are in H m loc (R d ) and periodic with period Ω. For any m ∈ N, H m per is a Hilbert space for the scalar product

((u, v)) m = |κ|≤m Ω D κ u(x)D κ v(x)dx
(κ being a multi-index) and its associated norm u m = ((u, u)) We denote by Ḣm per = {u ∈ H m per : u = 0} (m ∈ Z) the Sobolev spaces for functions with zero mean. We will frequently use the fact that H m per is isomorphic to R × Ḣm per (m ∈ Z) through the decomposition u = u + u.

1/2 m . For m = 0, H 0 per = L 2 (Ω),
Using the dense and continuous inclusions

H 1 per ⊂ L 2 (Ω) ⊂ H -1 per , the semi-scalar product on H 1 per , (u, v) → (∇u, ∇v), defines a linear operator A = -∆ : D(A) → L 2 (Ω) with domain D(A) = H 2 per . We denote Ȧ = -∆ : D( Ȧ) → L 2 (Ω) the restriction of A to L2 (Ω), with domain D( Ȧ) = Ḣ2
per . We observe that Ȧ is a positive self-adjoint operator with compact resolvent so that its powers Ȧs (s ∈ R) are well defined and it is possible to prove that Ḣm per = D( Ȧm/2 ) (m ∈ Z). For m = -1, we introduce an equivalent and more convenient norm

| • | -1 on Ḣ-1 per associated with the inner product ( u, v) -1 = ( Ȧ-1/2 u, Ȧ-1/2 v), so that for any u ∈ Ḣ-1 per , we have | u| -1 = Ȧ-1/2 u = ∇ Ȧ-1 u .
Similarly, for m = 1, we will sometimes use the equivalent norm

|•| 1 in H 1 per associated with the inner product (u, v) 1 = u v + (∇u, ∇v).

Moreover, Ȧ defines a continuous bijection from Ḣm

per onto Ḣm-2 per . In particular, for s = -1,

( u, v) -1 = A -1 u, v H 1 per ,H -1 per = u, A -1 v H -1 per ,H 1 per = (A -1/2 u, A -1/2 v).
2.2. Energy solutions. The nonlinearity f is a polynomial of odd degree whose leading coefficient is positive and which vanishes at 0:

f (s) = 2p+1 i=1 a i s i ∀s ∈ R (a 2p+1 > 0), (2.1)
with p ∈ N ⋆ if d = 1 or d = 2 and with p ∈ {1, 2} if d = 3. We denote F the antiderivative of f which vanishes at 0, i.e.,

F (s) = 2p+2 i=2 a i-1 i s i ∀s ∈ R. (2.2)
We will make use of the Sobolev injection H 1 per ⊂ L 2p+2 (Ω). In particular, there is a constant

C S = C S (Ω) such that u L 2p+2 (Ω) ≤ C S |u| 1 , ∀u ∈ H 1 per , (2.3
)

and the map v → f (v) is Lipschitz continuous on bounded sets of H 1 per with values into L (2p+2)/(2p+1) (Ω) ⊂ H -1
per . We also have H 2 per ⊂ C 0 (Ω) with continuous injection. Finally, we note that there exist constants c 1 ≥ 0, c 2 > 0 and c 3 ≥ 0 such that

F (s) ≥ 2s 2 -c 1 ∀s ∈ R, (2.4) |f (s)| ≤ c 2 F (s) + c 3 ∀s ∈ R.
(2.5)

We point out that the expression (2.2) includes the standard quartic potential (obtained with p = 1)

F (s) = s 4 4 + (1 -ε) 2 s 2 . (2.6)
In contrast to some authors, we do not assume that ε < 1 (we simply have ε ∈ R).

In [START_REF] Grasselli | Well-posedness and longtime behavior for the modified phase-field crystal equation[END_REF], a notion of energy solution was introduced. This is based on the following pseudo-energy:

E(u, v) = E(u) + β 2 | v| 2 -1 , (2.7)
which is well defined for any (u,

v) ∈ H 2 per × H -1 per . Definition 2.1. A pair (u, u t ) is called an energy solution to problem (1.1)-(1.2) if (u, u t ) ∈ L ∞ (R + ; H 2 per × H -1 per ), u tt ∈ L ∞ (R + ; H -4 per ), (2.8)
E(u, u t ) ∈ L ∞ (R + ), (2.9)
and the following relations hold:

βu tt + u t = 0, (2.10) Ȧ-1 (βu tt + u t ) + A 2 u -2Au + f (u) -f (u) = 0, in Ḣ-2 per , a.e. in R + , (2.11)
u(0) = u 0 in H 2 per , u t (0) = v 0 in H -1 per . (2.12)
Equation (2.10) can be interpreted as a conservation law for the mass. With this definition, we have Theorem 2.2 ( [START_REF] Grasselli | Well-posedness and longtime behavior for the modified phase-field crystal equation[END_REF]). For any intial data

(u 0 , v 0 ) ∈ H 2 per × H -1 per , problem (1.1)-(1.
2) has a unique global energy solution (u, u t ). Moreover, any energy solution satifies the strong time continuity property

u ∈ C 2 (R + ; H -4 per ) ∩ C 1 (R + ; H -1 per ) ∩ C(R + ; H 2 per )
, as well as the following energy identity, for all s, t ∈ R + with s < t,

E(u(t), u t (t)) = E(u(s), u t (s)) - t s | ut (τ )| 2 -1 dτ + t s v 0 e -τ /β Ω f (u(τ ))dxdτ.
(2.13)

In particular, whenever v 0 = 0 then the pseudo-energy is nonincreasing.

3. The space semi-discrete problem 3.1. The space semi-discrete scheme. Our space discretization is based on two ideas: first, in view of the time discretization, we write the PDE (1.1) as a first order system; second, in order to use low order finite elements, we split the tri-Laplacian into three terms, in the spirit of a well-known splitting approach in Cahn-Hilliard type equations (see, e.g., [START_REF] Elliott | The Cahn-Hilliard model for the kinetics of phase separation[END_REF][START_REF] Grasselli | A splitting method for the Cahn-Hilliard equation with inertial term[END_REF][START_REF] Galenko | Unconditionally stable and numerical solution of the hyperbolic phase-field crystal equation[END_REF]). We obtain the following system, which is (formally) equivalent to (1.1):

         u t = v βv t = -v + ∆w z = -∆u w = -∆z + 2∆u + f (u).
Now, let V h denote a finite-dimensional subspace of H 1 per which contains the constants. In applications, V h will be a space of conforming finite elements (see Section 6). The space V h can also be obtained with a spectral basis.

The space semi-discrete scheme reads: find u h , v h , z h , w h : R

+ → V h such that          (∂ t u h , ϕ h ) = (v h , ϕ h ) β(∂ t v h , ψ h ) = -(v h , ψ h ) -(∇w h , ∇ψ h ) (z h , ζ h ) = (∇u h , ∇ζ h ) (w h , ξ) = (∇z h , ∇ξ h ) -2(∇u h , ∇ξ h ) + (f (u h ), ξ h ), (3.1)
for all ϕ h , ψ h , ζ h , ξ h in V h . This problem is completed with initial conditions

u h (0) = u 0 h , v h (0) = v 0 h , (3.2)
where u 0 h and v 0 h are given in V h . It will be convenient to work with an appropriate basis of V h . For this purpose, let (e i h ) 1≤i≤N h denote an orthonormal basis of V h for the L 2 (Ω) scalar product, such that e 1 h ≡ 1. The integer N h is the dimension of V h . To every function r h = N h i=1 r i e i h ∈ V h corresponds a unique (column) vector R = (r 1 , . . . , r N h ) t , represented by the corresponding capital letter. We seek

u h (t) = N h i=1 u i (t)e i h ≃ (u 1 (t), . . . , u N h (t)) t = U (t), v h ≃ V, z h ≃ Z, w h ≃ W. Define A = (A ij ) 1≤i,j≤N h ,
where

A ij = (∇e i h , ∇e j h ), 1 ≤ i, j ≤ N h , (3.3)
and let

F h (U ) = (F (u h ), 1), so that ∇F h (U ) = (f (u h ), e 1 h ), . . . , (f (u h ), e N h h ) t .
By choosing the test functions ϕ h , ψ h , ζ h , ξ h in (3.1) as the basis functions e i h , we obtain the following equivalent system:

         U t = V βV t = -V -AW Z = AU W = AZ -2AU + ∇F h (U ). (3.4)
Eliminating V , Z and W , we see that (3.4) is equivalent to

βU tt + U t = -A[A 2 U -2AU + ∇F h (U )], t ≥ 0. (3.5)
Since A is a discretization of -∆, this is natural space semi-discrete version of (1.1). Let U denote a solution of (3.5). We notice that the first line and the first column of A are filled with zeros (recall e 1 h ≡ 0, so that ∇e 1 h ≡ 0). Thus, the first component of U , u 1 (t) = (u h (t), 1), satisfies

β∂ tt u 1 + ∂ t u 1 = 0, t ≥ 0. (3.6) Solving (3.6) with initial conditions u 1 (0) = (u 0 h , 1) =: u 0 1 and ∂ t u 1 (0) = (v 0 h , 1) =: v 0 1 yields ∂ t u 1 (t) = v 0 1 e -t/β =: a(t) and u 1 (t) = M -βa(t), with M = βv 0 1 + u 0 1 . (3.7)
For every vector R = (r 1 , . . . , r

N h ) t ∈ R N h , we denote Ṙ = (r 2 , . . . , r N h ) t ∈ R N h -1 . Then U satisfies β Utt + Ut = -Ȧ[ Ȧ2 U -2 Ȧ U + ∇F h (U )], t ≥ 0, (3.8)
where Ȧ is the submatrix Ȧ = (A ij ) 2≤i,j≤N h , and

∇F h (U ) = (f (u h ), e 2 h ), . . . , (f (u h ), e N h h ) t .
We can also write ∇F h (U ) = Ṗ (∇F h (u 1 (t), U )), where Ṗ : R N h → R N h -1 is the projection on the components 2, . . . , N h . This shows that ∇F h (U ) is a "non autonomous" function of U (recall that u 1 (t) is determined by (3.7)). For later purpose, we note that by (3.3), Ȧ is symmetric positive definite: in particular, Ȧ is invertible. Conversely, any solution U of (3.7)-(3.8) satisfies (3.4), i.e. that the second equation of (3.4) is satisfied with V , Z and W given by the three other equations of the system (3.4).

3.2. Existence, uniqueness, and discrete energy estimate. The standard Euclidean norm in R N h or R N h -1 will be denoted |•|. We also use the following quadratic norm:

| Ṙ| -1 = Ṙt Ȧ-1 Ṙ 1/2 , (3.9)
defined for all Ṙ ∈ R N h -1 . Notice that |A s U | = | Ȧs U | (s > 0, U ∈ R N h ). We set E h (U ) = 1 2 |AU | 2 -|A 1/2 U | 2 + F h (U ), (3.10)
E h (U, V ) = E h (U ) + β 2 | V | 2 -1 . (3.11)
As a shortcut, for a solution (U, U t ) of (3.5), we will write

E h (t) = E h (U (t), U t (t)).
Notice that by the Cauchy-Schwarz inequality we have

|A 1/2 U | 2 = U t AU ≤ 1 4 |AU | 2 + |U | 2 .
(3.12)

Then, using also (2.4), we find that

E h (U ) ≥ 1 4 |AU | 2 + |U | 2 -c 1 . (3.13)
We first prove the following

Lemma 3.1. Any solution U ∈ C 2 ([0, T ); R N h ) of (3.5) satisfies the energy equality d dt E h (t) + | Ut | 2 -1 = v 0 1 e -t/β (f (u h ), 1), (3.14)
and the energy estimate

E h (t) + t 0 | Ut (s)| 2 -1 ds ≤ E h (0)e 2c 2 |v 0 1 |β + (c 1 c 2 + c 3 )|v 0 1 |βe 2c 2 |v 0 1 |β , (3.15)
for all t ∈ [0, T ), where c 1 , c 2 and c 3 depend only on f (see (2.4)-(2.5)), and where

v 0 1 = ∂ t u 1 (0). Proof. Recall that ∂ t u 1 (t) = v 0 1 e -t/
d dt [F h (U (t))] = (∇F h (U ), U t (t)) = N h i=1 ∂ t u i (t)(f (u h ), e i h ),
we find the energy equality (3.14). Estimate (2.5) yields

d dt E h (t) + | Ut | 2 -1 ≤ |v 0 1 |e -t/β (c 2 F h (U ) + c 3 )
. By the Cauchy-Schwarz inequality, we have

2|A 1/2 U | 2 = 2U t AU ≤ |AU | 2 + |U | 2 .

Thus we get

2E h (U ) = (|AU | 2 -2|A 1/2 U | 2 + |U | 2 ) + F h (U ) + (F h (U ) -|U | 2 ) ≥ F h (U ) -c 1 , so that F h (U ) ≤ 2E h (U, V ) + c 1 . (3.16) Therefore we obtain d dt E h (t) + | Ut | 2 -1 ≤ |v 0 1 |e -t/β (2c 2 E h (t) + c 1 c 2 + c 3 ) . Letting η(t) = t 0 2c 2 |v 0 1 |e -t/β and c ′ 3 = c 1 c 2 + c 3 , Gronwall's lemma yields E h (t) + t 0 | Ut (s)| 2 -1 e η(t)-η(s) ds ≤ E h (0)e η(t) + t 0 c ′ 3 |v 0 1 |e -s/β e η(t)-η(s) ds.
Since η(t) = 2c 2 |v 0 1 |β(1-e -t/β ) ≤ 2c 2 |v 0 1 |β, we deduce the energy estimate (3.15). Theorem 3.2. For every U 0 , V 0 in R N h , there exists a unique solution U ∈ C 2 (R + , R N h ) of (3.5) such that U (0) = U 0 and U t (0) = V 0 . Proof. By the Cauchy-Lipschitz theorem, there exists a unique maximal solution U ∈ C 2 ([0, T + ); R N h ) of (3.5) satisfying the given initial conditions. The energy estimate (3.15) shows that E h is uniformly bounded for t ≥ 0. By (3.13), |U | and | Ut | -1 are uniformly bounded for t ≥ 0. This, together with the estimate (3.7) on the mass, implies that T + = +∞.

3.3. Some additional notation. We assume now that (V h ) h>0 is a family of subspaces of H 1 per such that: (H1) for all h > 0, V h has finite dimension and contains all the constants; (H2) for any ϕ ∈ H 1 per , there exists ϕ h ∈ V h such that ϕ h → ϕ (strongly) in H 1 per , as h tends to 0. For the convergence result as h → 0, it will be useful to have h-dependent operators and norms. We denote A h : V h → V h the linear operator such that for any

q h ∈ V h , A h q h solves (A h q h , ζ h ) = (∇q h , ∇ζ h ), ∀ζ h ∈ V h . (3.17) The operator A h is a discrete Laplacian, A h ≃ -∆ h . Notice that if q h is constant, then A h q h =
0 so that A h is not invertible. In order to define a discrete version of Ȧ-1 , we introduce the subspace

Vh = {ϕ h ∈ V h : ϕ h = 0} .
The bilinear form (∇•, ∇•) is symmetric positive definite on Vh . We can define the operator Ṡh : Vh → Vh such that for any ṙh ∈ Vh , Ṡh ṙh is the unique solution of

(∇ Ṡh ṙh , ∇ φh ) = ( ṙh , φh ), ∀ φh ∈ Vh . (3.18)
By choosing ζ h ≡ 1 in (3.17), we see that A h (V h ) ⊂ Vh , so that the restriction Ȧh : Vh → Vh of A h is well defined. Using (3.17)-(3.18), it is easily seen that Ṡh = Ȧ-1 h . We also define the L 2 -orthogonal projector

P h : L 2 (Ω) → V h , i.e., (P h q, ϕ h ) = (q, ϕ h ), ∀q ∈ L 2 (Ω), ∀ϕ h ∈ V h .
By Pythagoras' theorem and assumption (H2), for any q ∈ L 2 (Ω), we have

q -P h q = inf r h ∈V h q -r h → 0, as h → 0. (3.19)
Since V h ⊂ H 1 per , the operator P h has a natural extension to H -1 per (also denoted P h ), by setting

P h q ∈ V h , (P h q, ϕ h ) = q, ϕ h H -1 per ,H 1 per , ∀q ∈ H -1 per , ∀ϕ h ∈ H 1 per .
The H 1 -orthogonal projector Π h : H 1 per → V h is defined as follows: for any q ∈ H 1 per , Π h q ∈ V h is uniquely defined by

Π h q = q and (∇Π h q, ∇ϕ h ) = (∇q, ∇ϕ h ), ∀q ∈ H 1 per , ∀ϕ h ∈ V h . (3.

20) By Pythagora's theorem and assumption (H2), for any q

∈ H 1 per , we get |q -Π h q| 1 = inf r h ∈V h |q -r h | 1 → 0, as h → 0. ( 3.21)
We point out that the space Vh is invariant by Π h and by P h . By using a L 2 -orthonormal basis (e i h ) 1≤i≤N h of V h as in Section 3.1, with e 1 h ≡ 1, we see that the matrix of A h is A, so that the energy E h from (3.10) can be rewritten

E h (u h ) = 1 2 A h u h 2 -| uh | 2 1 + (F (u h ), 1); the norm | • | -1 from (3.9) becomes for any element ṙh ∈ Vh | ṙh | -1,h := ( ṙh , Ȧ-1 h ṙh ) = | Ȧ-1 h ṙh | 1 .

It is an Euclidean norm for the following scalar product on

Vh ( ṙh , qh ) -1,h := ( ṙh , Ȧ-1 h qh ) = ( Ȧ-1 h ṙr , Ȧ-1 h qh ) 1 . Thus, the discrete pseudo-energy (3.11) reads E h (u h , v h ) := E h (u h ) + β 2 | vh | 2 -1,h , ∀u h , v h ∈ V h . (3.22)
The Cauchy-Schwarz (3.12) inequality gives

| u| 2 1 ≤ 1 4 Ȧh u h 2 + u h 2 , ∀u h ∈ V h , (3.23)
and estimate (3.13) becomes

E h (u h ) ≥ 1 4 A h u h 2 + u h 2 -c 1 , ∀u h ∈ V h . (3.24) 3.4. Convergence as h → 0. Theorem 3.3. Let (u 0 , v 0 ) ∈ H 2 per × H -1 per . Assume that (u 0 h , v 0 h ) h>0 is a family of functions in V h × V h such that u 0 h → u 0 in H 1 per , A h u 0 h → Au 0 in L 2 (Ω), (3.25)
v 0 h → v 0 in R, Ȧ-1 h v0 h → Ȧ-1 v0 in Ḣ1 per , (3.26)
as h → 0. Then the solution (u h , ∂ t u h ) of the space semi-discrete scheme (3.1)-(3.2) tends to the energy solution (u, u t ) of (1.1)-(1.2) in the following sense:

u h → u weakly ⋆ in L ∞ (R + ; H 1 per ), u h → u strongly in C([0, T ]; L 2 (Ω)), for all T > 0, A h u h → Au weakly ⋆ in L ∞ (R + ; L 2 (Ω)), Ȧ-1 h ∂ t uh → Ȧ-1 ∂ t u weakly ⋆ in L ∞ (R + ; Ḣ1
per ) and weakly in L 2 (R + ; Ḣ1 per ). Proof. The idea is to use a priori estimates on the mass and on the discrete energy, and to pass to the limit in the equation by a compactness argument. We first consider the conservation law for the mass. By (3.7), we get

a h (t) := ∂ t u h (t) = v 0 h e -t/β , t ≥ 0, and u h (t) = β v 0 h + u 0 h -βa h (t), t ≥ 0. (3.27) By assumption, u 0 h → u 0 and v 0 h → v 0 in R,
so that a h converges uniformly on R + to the function a(t) := v 0 e -t/β , and u h converges uniformly on R + to the function β v 0 + u 0 -βa(t). The estimates below show that (u h ) h>0 is bounded in L ∞ (R + ; L 2 (Ω)), so that, up to a subsequence, u h converges weakly ⋆ in L ∞ (R + ; L 2 (Ω)) to some u and so u h → u weakly ⋆ in L ∞ (R + ). By uniqueness of the limit, we find u = β v 0 + u 0 -βa(t). By differentiating, we recover the conservation law for the mass:

β∂ tt u + ∂ t u = 0, t ≥ 0.
We now turn to the energy estimate. As pointed out in Section 3.1, the (unique)

solution (u h , ∂ t u h ) of (3.1)-(3.2) is in fact a solution (u h , v h , z h , w h ) of (3.1)-(3.2). In particular, v h = ∂ t u h and z h = A h u h . We have (recall (3.22)): E(u 0 h , v 0 h ) = 1 2 A h u 0 h 2 -| u0 h | 2 1 + (F (u 0 h ), 1) + β 2 | Ȧ-1 h v0 h | 2 1 .

By using assumptions (3.25)-(3.26) and the Sobolev injection

H 1 per ֒→ L 2p+2 (Ω), we see that E(u 0 h , v 0 h
) is uniformly bounded as h tends to 0. The energy estimate (3.15) shows that there exists a constant C independent of h such that

E h (u h (t), ∂ t u h (t)) + t 0 |∂ t uh (s)| 2 -1,h ds ≤ C,
for all t ≥ 0. By (3.24) and (3.23), we obtain that z h = A h u h and u h are uniformly bounded in L 2 (Ω), that uh and ṙh := Ȧ-1 h ∂ t uh are uniformly bounded in Ḣ1 per , and that

∞ 0 | Ȧ-1 h ∂ t uh (t)| 2 1 dt ≤ C.
This implies that (u h ) h>0 is precompact in the space C([0, T ]; L 2 (Ω)), for all T > 0, as a consequence of the Ascoli-Arzelà Theorem. Indeed, let T > 0. The family (u h ) h>0 is uniformly bounded from [0, T ] with values in H 1 per , and H 1 per is compactly embedded into L 2 (Ω) by Rellich's Theorem. Moreover, for all 0 ≤ s ≤ t ≤ T , we have

uh (t) -uh (s) 2 = 2 t s (∂ t uh (σ), uh (σ) -uh (s))dσ = 2 t s (∇ ṙh (σ), ∇[ uh (σ) -uh (s)])dσ ≤ 4 ṙh L ∞ (R + ;H 1 per ) uh L ∞ (R + ;H 1 per ) |t -s|. (3.28)
Moreover, by (3.27) and the mean value theorem, we find

| u h (t) -u h (s) | ≤ | v 0 h ||t -s|.
Thus, (u h ) h>0 is uniformly equicontinuous from [0, T ] with values into L 2 (Ω), and therefore pre-compact in C([0, T ]; L 2 (Ω)), as claimed. Up to a subsequence, we have the following convergence results

u h → u weakly ⋆ in L ∞ (R + ; H 1 per ), u h → u strongly in C([0, T ]; L 2 (Ω)), for all T > 0, u h → u a.e. in Ω × R + , f (u h ) → f (u) weakly in L q (0, T ; L q (Ω)), for all T > 0, z h → z weakly ⋆ in L ∞ (R + ; L 2 (Ω)), ṙh → ṙ weakly ⋆ in L ∞ (R + ; Ḣ1 per ), ṙh → ṙ weakly in L 2 (R + ; Ḣ1 per ),
where q = (2p + 2)/(2p + 1) > 1. Let now ψ ∈ Ḣ1 per and let ψh = Π h ( ψ) so that ψh → ψ strongly in Ḣ1

per . We have

(∂ t uh , ψh ) = ( Ȧh ṙh , ψh ) = (∇ ṙh , ∇ ψh ) → ( ṙ, ψ) weakly ⋆ in L ∞ (R + ). On the other hand, ∂ t ( uh , ψh ) → ∂ t ( u, ψ) in D ′ (0, ∞) (i.e. in the sense of distributions), since ( uh , ψh ) → ( u, ψ) in L ∞ (R + ) weakly ⋆. Thus, ∂ t ( u, ψ) = (∇ ṙ, ∇ ψ) = v, ψ H -1 per ,H 1 per , (3.29
)

with v = Ȧ ṙ ∈ L ∞ (R + ; Ḣ-1 per ). This shows that ∂ t u = v ∈ L ∞ (R + ; Ḣ-1 per ). Next, we set ϕ ∈ H 2 per and we let ϕ h = Π h (ϕ), so that ϕ h → ϕ strongly in H 1 per . Let (e i h ) 1≤i≤N h be an orthonormal basis of V h with e 1 h ≡ 1. We let ϕ h = N h
i=1 ϕ i e i h and Φ = (ϕ 1 , . . . , ϕ N h ) t be the vector associated to ϕ h , as in Section 3.1. On multiplying (3.8) by Φt Ȧ-1 and using z h = A h u h , we find

β(∂ tt uh , φh ) -1,h + (∂ t uh , φh ) -1,h + (∇z h , ∇ϕ h ) -2(∇u h , ∇ϕ h ) + (f (u h ), φh ) = 0,
(3.30) for all t ≥ 0. We have that

(∂ t uh , φh ) -1,h = ( ṙh , φh ) → ( ṙ, φ) = (∂ t u, φ) -1 , by (3.29). The convergence above holds in L ∞ (R + ) weak ⋆, so that (∂ tt uh , φh ) -1,h = ∂ t (∂ t uh , φh ) -1,h → ∂ t (∂ t u, φ) -1 in D ′ (0, ∞). Moreover, (∇u h , ∇ϕ h ) → (∇u, ∇ϕ) in L ∞ (R + ) weak ⋆ . Since ϕ h = Π h ϕ, we have (∇z h , ∇ϕ h ) = (∇z h , ∇ϕ) = (z h , Aϕ) → (z, Aϕ), in L ∞ (R + ) weak ⋆. Concerning the last term in (3.30), we have (f (u h ), φh ) → (f (u), φ) weakly in L q (0, T ), ∀T > 0.
Summing up, we have proved that

β∂ t (∂ t u, φ) -1 + (∂ t u, φ) -1 + (z, Aϕ) -2(∇u, ∇ϕ) + (f (u), ϕ) = f (u) ϕ . (3.31)
The equality holds in

D ′ (0, ∞), for all ϕ ∈ H 2 per . Moreover, ∂ t u ∈ L ∞ (R + ; H -1 per ), z ∈ L ∞ (R + ; L 2 (Ω)), u ∈ L ∞ (R + ; H 1 per ) and f (u) ∈ L ∞ (R + ; H -1 per ) so that ∂ tt u belongs to L ∞ (R + ; H -4 per ) and Ȧ-1 (β∂ tt u + ∂ t u) + Az -2Au + f (u) -f (u) = 0 in Ḣ-2 per , a.e. in R + . Now, recall that z h = A h u h . Let ζ ∈ H 1 per and ζ h = Π h ζ, so that ζ h → ζ strongly in H 1 per . We have (z h , ζ h ) = (A h u h , ζ h ) = (∇u h , ∇ζ h ) → (∇u, ∇ζ),
on one hand, and (z h , ζ h ) → (z, ζ), on the other hand. Thus, we deduce (z, ζ) = (∇u, ∇ζ), in L ∞ (R + ). The equality holds for every ζ ∈ H 1 per , so z = Au, u ∈ L ∞ (R + ; H 2 per ) and (u, u t ) is an energy solution of (1.1)-(1.2). By uniqueness of the limit (u, u t ), the whole family (u h , ∂ t u h ) converges to (u, u t ).

Remark 3.4. Let (u 0 , v 0 ) ∈ H 2 per × H -1 per . If u 0 h = Π h (u 0 ) and v 0 h = P h (v 0 ), then assumptions (3.25)-(3.26) are satisfied. Indeed, for all ϕ h ∈ V h , (Au 0 , ϕ h ) = (∇u 0 , ∇ϕ h ) = (∇Π h u 0 , ∇ϕ h ) = (A h (Π h u 0 ), ϕ h ).
Then A h (Π h u 0 ) = P h (Au 0 ) and by (3.19) we obtain

A h (Π h u 0 ) → Au 0 in L 2 (Ω), as h → 0.

By definition, observe that

P h (v 0 ) = (P h (v 0 ), 1) = v 0 , 1 H -1 per ,H 1 per = v 0 .
Moreover, for all ϕ h ∈ V h , we have

(P h ( v0 ), φh ) = v0 , φh H -1 per ,H 1 per = (∇ Ȧ-1 v0 , ∇ φh ) = (∇Π h ( Ȧ-1 v0 ), ∇ φh ) = ( Ȧh (Π h ( Ȧ-1 v0 )), φh ), so that P h ( v0 ) = Ȧh (Π h ( Ȧ-1 v0 )). Thus, thanks to (3.21), we deduce Ȧ-1 h (P h ( v0 )) = Π h ( Ȧ-1 v0 ) → Ȧ-1 v0 in H 1 per .

The fully discrete problem

4.1. The fully discrete scheme. For the time discretization, we use the following decomposition:

(H3) F = F + + F -, where F + and F -are polynomials such that As a consequence, deg(F + ) = deg(F ) and F + , F have the same leading coefficient. We denote f = F ′ = f + +f -, where f + = F ′ + and f -= F ′ -. For the energy estimate, we will use the fact that there exist two constants c 5 > 0 and c 6 ≥ 0 which depend only on f and on the decomposition

F (iv) + ≥ 0, F (iv)
f = f + + f -such that 1 2 (|f (r)| + |f (s)|) + 1 12 (s -r) 2 (|f ′′ + (r)| + |f ′′ -(s)|) ≤ c 5 (F (r) + F (s)) + c 6 , (4.1)
for all r, s ∈ R.

Remark 4.1. A decomposition (H3) is always possible for a polynomial potential such as (2.2). Indeed, for a quartic polynomial (for instance (2.6)), we can always choose F + = F and F -= 0. For a polynomial with higher degree, we notice that F (iv) , being a polynomial of even degree with strictly positive leading coefficient, is bounded from below, i.e.

F (iv) (s) ≥ -c 4 ∀s ∈ R,
for some constant c 4 ≥ 0. A possible (but not unique !) choice is then

F (iv) + = F (iv) + c 4 and F (iv) - = -c 4 , i.e. F + (s) = F (s) + c 4 s 4 /24 and F -(s) = -c 4 s 4 /24.
We use the same notation as in Section 3. In particular, V h is a family of finitedimensional subspaces of H 1 per which satisfies assumptions (H1)-(H2). Let τ > 0 denote the time step, and (u 0 h , v 0 h ) in V h × V h the initial datum. The fully discrete scheme reads: for n ≥ 0, find (u n+1 , v n+1 , z n+1 , w n+1 )

∈ (V h ) 4 such that                    ((u n+1 h -u n h )/τ, ϕ h ) = (v n+1/2 h , ϕ h ) β((v n+1 h -v n h )/τ, ψ h ) = -(v n+1/2 h , ψ h ) -(∇w n+1 h , ∇ψ h ) (z n+1 h , ζ h ) = (∇u n+1/2 h , ∇ζ h ) (w n+1 h , ξ h ) = (∇z n+1 h , ∇ξ h ) -2(∇u n+1/2 h , ∇ξ h ) + ((f (u n h ) + f (u n+1 h ))/2, ξ h ) - 1 12 (u n+1 h -u n h ) 2 f ′′ + (u n h ) + f ′′ -(u n+1 h ) , ξ h , (4.2) for all ϕ h , ψ h , ζ h , ξ h in V h .
Here, we have denoted

u n+1/2 h = (u n+1 h + u n h )/2 and v n+1/2 h = (v n+1 h + v n h )/2.
Notice that z 0 h and w 0 h do not need to be defined. In fact, z n+1 h (resp. (w n+1 h )) is a second-order (in time) approximation of z h (t n+1/2 ) (resp. w h (t n+1/2 )), where t n+1/2 = (n + 1/2)τ .

Let (e i h ) 1≤i≤N h be an L 2 -orthonormal basis of V h , with e 1 h ≡ 1, so that we have the identification

V h ∋ u h ≃ U ∈ R N h . In R N h , the scheme reads: let U 0 , V 0 in R N h and for n ≥ 0 find (U n+1 , V n+1 , Z n+1 , W n+1 ) ∈ (R N h ) 4 which solves                (U n+1 -U n )/τ = V n+1/2 β(V n+1 -V n )/τ = -V n+1/2 -AW n+1 Z n+1 = AU n+1/2 W n+1 = AZ n+1 -2AU n+1/2 + ∇F h (U n ) + ∇F h (U n+1 ) /2 -G(U n , U n+1 ), (4.3)
where

G(U n , U n+1 ) = 1 12 (u n+1 h -u n h ) 2 f ′′ + (u n h ) + f ′′ -(u n+1 h ) , e i h 1≤i≤N h . (4.4)
On eliminating Z n+1 and W n+1 , the scheme becomes

     (U n+1 -U n )/τ -V n+1/2 = 0 β(V n+1 -V n )/τ + V n+1/2 + A A 2 U n+1/2 -2AU n+1/2 +A ∇F h (U n ) + ∇F h (U n+1 ) /2 -G(U n , U n+1 ) = 0. (4.5)
In Section 6.1, a numerical example indicates that our fully discrete scheme (4.2) has a second order convergence error in time (and also in space if V h is the space of P 1 finite elements). By arguing as in Gomez and Hughes [START_REF] Gomez | Provably unconditionally stable, second-order time-accurate, mixed variational methods for phase-field models[END_REF], we check here that:

Proposition 4.2. The scheme has a second order consistency error in time, i.e. that any solution of the space semi-discrete problem (3.4) satisfies the fully discrete scheme (4.5) with order O(τ 2 ).

Proof. Let (U, V ) be a solution of (3.4) on a finite time interval [0, T ]. Since f is a polynomial, by a bootstrap argument, we know that (U, V ) ∈ C ∞ ([0, T ]; V h × V h ). In the time discrete scheme (4.5), we replace

U n , U n+1 , V n , V n+1 by U (t n ), U (t n+1 ), V (t n ), V (t n+1
) respectively. The purpose of this replacement is to find a local truncation error (or consistency error) in O(τ 2 ) in the right-hand side of (4.5), instead of (0, 0) t . Now consider the midpoint scheme, which is the same as (4.5) without the term G(U n , U n+1 ). Standard results show that the midpoint scheme has a second order consistency error, so that it is sufficient to show that G(U

(t n ), U (t n+1 )) = O(τ 2 )
. This is obvious by definition (4.4). Indeed, using the assumption on p and appropriate Hölder inequalities and Sobolev injections, we have

|G(U (t n ), U (t n+1 ))| ≤ C u(t n+1 ) -u(t n ) 2 1 . Moreover, by Taylor expansion, (u h (t n+1 )-u h (t n )) 2 1 = O(τ 2
), so the local truncation error is O(τ 2 ) as well. Notice that the constant in the consistency error depends on h, on T , and on maximum norms of derivatives of U , V up to order 3. 4 which complies with (4.2). Proof. We work with the R N h version (4.5). We will show that this problem is variational, and that we can find U n+1 by a minimization procedure. Let (U n , V n ) be fixed in R N h . Consider the polynomial of two variables

Existence, discrete energy estimate and uniqueness. Let us prove the following

(u 0 h , v 0 h ) ∈ V h × V h , there exists at least one sequence (u n h , v n h , z n h , w n h ) n≥1 in (V h )
g(r, s) = 1 12 (s -r) 2 (f ′′ + (r) + f ′′ -(s)) (r, s ∈ R).
By assumption (H3), we have deg(f -) < deg(f ) and deg(f

+ ) = deg(f ) so deg(f ′′ -) < deg(f ) -2, and deg(f ′′ + ) = deg(f ) -2.
Thus, g is a polynomial of total degree less than or equal to 2p + 1, and its partial degree with respect to the variable s is strictly less than 2p + 1. We can write

g(r, s) = k,l b k,l r k s l , (4.6)
for coefficients b k,l ∈ R, where 0 ≤ k ≤ 2p + 1 and 0 ≤ l < 2p + 1. Let us set now

h(r, s) = k,l b k,l r k s l+1 l + 1 , (4.7)
so that ∂ s h(r, s) = g(r, s). We define

H n h (U ) = (h(u n h , u h ), 1) with u h ≃ U , so that ∇H n h (U ) = (g(u n h , u h ), e i ) 1≤i≤N h = G(U n , U )
. By (4.7) and Hölder's inequality, we get

|H n h (U )| ≤ C n u h 2p+1 L 2p+2 (Ω) + 1 ∀U ∈ R N h , (4.8)
where the constant C n depends on u n h L 2p+2 (Ω) . Now, by eliminating V n+1 , we find that (4.5) is equivalent to

β τ 2(U n+1 -U n ) τ -V n + U n+1 -U n τ + A A 2 (U n+1 + U n ) 2 -2A (U n+1 + U n ) 2 + ∇F h (U n ) + ∇F h (U n+1 ) 2 -∇H n h (U n+1 ) = 0.
Writing U = (u 1 , U ), we see that this is equivalent to

β τ 2(u n+1 1 -u n 1 ) τ -v n 1 + u n+1 1 -u n 1 τ = 0, (4.9) β Ȧ-1 τ 2( U n+1 -U n ) τ -V n + Ȧ-1 ( U n+1 -U n ) τ + Ȧ2 ( U n+1 + U n ) 2 -2 Ȧ ( U n+1 + U n ) 2 + ∇F h (U n ) + ∇F h (U n+1 ) 2 -∇H n h (U n+1 ) = 0. (4.10)
The first equation determines u n+1 1 uniquely. The second equation can be solved by letting U n+1 be a minimizer on R N h -1 of the function:

G : U → β τ 2 | U -U n | 2 -1 - β τ ( V n ) t Ȧ-1 U + 1 2τ | U -U n | 2 -1 + 1 4 | Ȧ( U + U n )| 2 - 1 2 | Ȧ1/2 ( U + U n )| 2 + ( ∇F h (U n )) t 2 U + F n h ( U ) 2 -Hn h ( U ), where F n h ( U ) = F h (u n+1 1 , U ), Hn h ( U ) = H n h (u n+1 1 , U ).
By (2.2), we deduce

F n h ( U ) ≥ a 2p+1 2p + 2 u h 2p+2 L 2p+2 (Ω) -C ′ p ( u h 2p+1 L 2p+2 (Ω) + 1), ∀ U ∈ R N h -1 ,
where the constant C ′ p depends only on the coefficients of F . Thus, by (4.8), we find

G( U ) ≥ a 2p+1 2(2p + 2) u h 2p+2 L 2p+2 (Ω) -(C n + C ′ p 2 ) u h 2p+1 L 2p+2 (Ω) + c uh 2 -c ′ uh -c ′′ ,
where c > 0 and c ′ , c ′′ ≥ 0 depend on h, F and u n h . For the quadratic term, we used that all norms are equivalent in R N h -1 . As a consequence, G( U ) → +∞ as | U | → +∞: the continuous function G has a minimizer in R N h -1 , and the proof is complete.

The behavior of (u n 1 , v n 1 ) is straightforward, thanks to a discrete conservation law for the mass. Indeed, choosing ψ = 1 in the second equation of (4.2), we find

v n+1 1 = qv n 1 , with q = q(β, τ ) = 2β -τ 2β + τ . (4.11)
Thus, we obtain v n 1 = q n v 0 1 . (4.12) We also have v n+1/2 = qv n-1/2 , so that v n+1/2 = q n v 1/2 . Notice that |q| < 1, since β > 0 and τ > 0, so that v n 1 → 0. If τ > 2β, then q < 0. On choosing ϕ = 1 in the first equation of (4.2), we find

u n+1 1 = u n 1 + τ v n+1/2 1
.

By induction, we deduce

u n 1 = u 0 1 + τ n-1 k=0 q k v 1/2 1 = u 0 1 + τ 1 -q n 1 -q v 1/2 1 . (4.13)
For the energy estimate, we will need a technical lemma, adapted from [START_REF] Gomez | Provably unconditionally stable, second-order time-accurate, mixed variational methods for phase-field models[END_REF]: k + 2 (σ)g (3) (σ)dσ, (4.14)

1 0 g(s)ds = 1 2 (g(0) + g(1)) - 1 12 g ′′ (1) + 1 2 1 0 k - 2 (σ)g (3) (σ)dσ, (4.15)
where k + 2 (σ) = (1-σ) 2 (2σ+1)/6 and k - 2 (σ) = σ 2 (3-2σ)/6. In particular, k + 2 (σ) ≥ 0 and k - 2 (σ) ≥ 0 for all σ ∈ [0, 1]. Proof. We prove (4.15) (the proof of (4.14) is similar). For a function ϕ ∈ C 2 ([0, 1]), let

Err(ϕ) = 1 0 ϕ(s)ds - 1 2 (ϕ(0) + ϕ(1)) - 1 12 ϕ ′′ (1)
denote the error of the quadrature formula. If ϕ is a polynomial of degree ≤ 2, a direct computation shows that Err(ϕ) = 0. Now, let g ∈ C 3 ([0, 1]). The Taylor formula of order 2 at s = 0 reads

g(s) = p 2 (s) + 1 2 1 0 (s -σ) 2 + g (3) (σ)dσ, with p 2 (s) = g(0) + sg ′ (0) + s 2 g ′′ (0)/2 and (s -σ) + = s -σ if s ≥ σ 0 if s ≤ σ.
In particular, Err(p 2 ) = 0. By linearity of Err,

Err(g) = Err(p 2 ) + Err s → 1 2 1 0 (s -σ) 2 + g (3) (σ)dσ .
By inverting the integration signs, we obtain

Err(g) = 1 2 1 0 k - 2 (σ)g (3) (σ)dσ,
where k - 2 (s) = Err(s → (s -σ) 2 + ). Using the definition of Err, we find that for σ ∈ [0, 1),

k - 2 (σ) = 1 0 (s -σ) 2 + ds - 1 2 [(0 -σ) 2 + + (1 -σ) 2 + ] - 2 12 , = σ 2 /2 -σ 3 /3.
The claim is proved.

We have (compare with Lemma 3.1):

Lemma 4.5 (Energy estimate for any τ ). 4 which complies with (4.3), then for all n ≥ 0,

If (U n , V n , Z n , W n ) n≥1 is a sequence in (R N h )
E h (U n+1 , V n+1 ) -E h (U n , V n ) + τ | V n+1/2 | 2 -1 ≤ τ v n+1/2 1 w n+1 1 . (4.16)
As a consequence, for all k ≥ 0, we have

E h (U N 0 +k , V N 0 +k ) + k-1 j=0 τ | V N 0 +j+1/2 | 2 -1 ≤ exp 16c 5 τ |q| N 0 1 -|q| |v 1/2 1 | E h (U N 0 , V N 0 ) + c 7 τ |q| N 0 1 -|q| |v 1/2 1 | , (4.17)
where

N 0 = N 0 (β, c 5 , τ, |v 0 1 |) ∈ N is such that 2c 5 τ |q| N 0 |v 1/2 1 | ≤ 1/2, (4.18)
c 7 = 2c 1 c 5 +c 6 depends only on f , f + , f -(see (2.4), (4.1)), and q is defined by (4.11).

Proof. Let δu

n h = u n+1 h -u n h . Since f + = F ′ + and f -= F ′ -, we have F + (u n+1 h ) -F + (u n h ) = δu n h 1 0 f + (u n h + sδu n h )ds, (4.19)
F -(u n+1 h ) -F -(u n h ) = δu n h 1 0 f -(u n h + sδu n h )ds. (4.20) Choosing g(s) = f + (u n h + sδu n h ) in (4.14), we find 1 0 f + (u n h + sδu n h )ds = 1 2 (f + (u n h ) + f + (u n+1 h )) - (δu n h ) 2 12 f ′′ + (u n h ) - (δu n h ) 3 2 1 0 k + 2 (σ)f ′′′ + (u n h + σδu n h )dσ. (4.21) Setting g(s) = f -(u n h + sδu n h ) in (4.15), we find 1 0 f -(u n h + sδu n h )ds = 1 2 (f -(u n h ) + f -(u n+1 h )) - (δu n h ) 2 12 f ′′ -(u n+1 h ) + (δu n h ) 3 2 1 0 k - 2 (σ)f ′′′ -(u n h + σδu n h )dσ. ( 4
F (u n+1 h ) -F (u n h ) = δu n h 1 2 (f (u n h ) + f (u n+1 h )) - (δu n h ) 2 12 (f ′′ + (u n h ) + f ′′ -(u n+1 h )) -α n ,
where

α n = (δu n h ) 4 2 1 0 k + 2 (σ)f ′′′ + (u n h + σδu n h )dσ - 1 0 k - 2 (σ)f ′′′ -(u n h + σδu n h )dσ ≥ 0.
by assumption (H3) on the decomposition. Next, we choose ξ h = δu n h in the last equation of (4.2). This gives

n+1 h , δu n h ) -(∇z n+1 h , ∇δu n h ) + 2(∇u n+1/2 h , ∇δu n h ) = F (u n+1 h ) -F (u n h ) + (α n , 1).
Using the vector form with δU n = U n+1 -U n , and eliminating z n+1 h , we obtain

F h (U n+1 ) -F h (U n ) + (α n , 1) = (W n+1 ) t δU n - 1 2 (|AU n+1 | 2 -|AU n | 2 ) +|A 1/2 U n+1 | 2 -|A 1/2 U n | 2 . (4.23)
The second equation in (4.2) implies

-Ẇ = Ȧ-1 β ( V n+1 -V n ) τ + V n+1/2 .
Plugging this in (4.23), together with δU n = τ V n+1/2 , we get

E h (U n+1 ) -E h (U n ) + (α n , 1) + β 2 | V n+1 | 2 -1 -| V n | 2 -1 + τ | V n+1/2 | 2 -1 = τ v n+1/2 1 w n+1 1 .
This yields the energy estimate (4.16).

Choosing r = 1 in the last equation of (4.2), we find

w n+1 1 = 1 2 (f (u n h ) + f (u n+1 h ), 1) - 1 12 (u n+1 h -u n h ) 2 f ′′ + (u n h ) + f ′′ -(u n+1 h ) , 1 .
Thus, by (4.1), we deduce

|w n+1 1 | ≤ 1 2 (|f (u n h )| + |f (u n+1 h )|) + 1 12 ((u n+1 h -u n h ) 2 (|f ′′ + (u n h )| + |f ′′ -(u n+1 h )|), 1) ≤ c 5 (F (u n h ) + F (u n+1 h )) + c 6 , 1 ,
As a consequence, by (3.16), we get

E h (U n+1 , V n+1 ) -E h (U n , V n ) + τ | V n+1/2 | 2 -1 ≤ τ |v n+1/2 1 | c 5 (F (u n h ), 1) + c 5 (F (u n+1 h ), 1) + c 6 ≤ τ |v n+1/2 1 | 2c 5 E h (U n+1 ) + 2c 5 E h (U n ) + 2c 5 c 1 + c 6

Let us set

E n h = E h (U n ) + β 2 | V n | 2 -1 .
So far, we have proved that

E n+1 h + τ | V n+1/2 | 2 -1 ≤ E n h + τ |q| n |v 1/2 1 | 2c 5 E n h + 2c 5 E n+1 h + c 7 ,
where

c 7 = 2c 5 c 1 + c 6 . Let N 0 = N 0 (β, c 5 , τ, |v 0 1 |) ∈ N satisfy (4.18). Then for n ≥ N 0 , we have (1 -2c 5 τ |q| n |v 1/2 1 |)E n+1 h + τ | V n+1/2 | 2 -1 ≤ (1 + 2c 5 τ |q| n |v 1/2 1 |)E n h + c 7 τ |q| n |v 1/2 1 |.
We divide by this inequality (1

-2c 5 τ |q| n |v 1/2
1 |) and we use that (by the mean value inequality) for all x ∈ [0, 1/2],

1 ≤ 1 1 -x and 1 + x 1 -x ≤ 1 + 8x ≤ exp(8x).
We obtain

E n+1 h + τ | V n+1/2 | 2 -1 ≤ exp(16c 5 τ |q| n |v 1/2 1 |) E n h + c 7 τ |q| n |v 1/2 1 | ,
for all n ≥ N 0 . By induction, for all k ∈ N, we deduce

E N 0 +k h + k-1 j=0 τ | V N 0 +j+1/2 | 2 -1 ≤ exp   16c 5 τ |v 1/2 1 | k-1 j=0 |q| N 0 +j   E N 0 h + k-1 j=0 exp 16c 5 τ |v 1/2 1 |(|q| N 0 + • • • + |q| N 0 +k-1-j ) c 7 τ |q| N 0 +j |v 1/2 1 |.
Estimate (4.17) follows by using the inequality k-1 j=0 |q| N 0 +j ≤ |q| N 0 /(1 -|q|).

Theorem 4.6 (Uniqueness for small τ ). For any

(u 0 h , v 0 h ) ∈ V h ×V h , there exists τ ⋆ = τ ⋆ (h) > 0 such that for any τ ∈ (0, τ ⋆ ), there is a unique sequence (u n h , v n h , z n h , w n h ) n≥1 which complies with (4.2). Moreover, τ ⋆ can be made independent of h if (u 0 h , v 0 h ) h>0 is a family such that | v 0 h | + E h (u 0 h , v 0 h ) ≤ C 1 , (4.24)
for some constant C 1 independent of h.

Proof. Assume that (u n h , v n h ) is uniquely determined for some n ≥ 0. We have seen that u n+1 [START_REF] Grasselli | Well-posedness and longtime behavior for the modified phase-field crystal equation[END_REF]) By (2.1), f ′ is a polynomial of even degree with strictly positive leading coefficient, so that f ′ is bounded from below. There exists an (optimal

U = U n+1 - U n+1 ≃ δ uh = δu h . We obtain 2β τ 2 |δ uh | 2 -1,h + 1 τ |δ uh | 2 -1,h + 1 2 A h δ uh 2 -|δ uh | 2 1 + 1 2 (f (u n+1 h ) -f (u n+1 h ), δu h ) -(g(u n h , u n+1 h ) -g(u n h , u n+1 h ), δu h ) = 0.(4.
) constant c f ≥ 0 such that f ′ (s) ≥ -c f ∀ s ∈ R. (4.26) By the mean value theorem, (f (u n+1 h ) -f (u n+1 h ), δu h ) ≥ -c f δ uh 2 .
On the other hand, by (4.6), we have

g(u n h , u n+1 h ) -g(u n h , u n+1 h ) = 0≤k+l≤2p+1 b k,l (u n h ) k [(u n+1 h ) l -(u n+1 h ) l ],
so that by Hölder's inequality,

(g(u n h , u n+1 h ) -g(u n h , u n+1 h ), δu h ) ≤ C ′ n δ uh 2 L 2p+2 (Ω) , where C ′ n = C ′ u n h L 2p+2 (Ω) , u n+1 h L 2p+2 (Ω) , u n+1 h L 2p+2 (Ω)
, and C ′ is a nondecreasing function of its arguments. Thus, by (2.3), equation ( 4

.25) implies 2β τ 2 |δ uh | 2 -1,h + 1 τ |δ uh | 2 -1,h + 1 2 A h δ uh 2 ≤ c f 2 δ uh 2 + (C ′ n C S + 1)|δ uh | 2 1 . (4.27)
Let (u 0 h , v 0 h ) be a given initial data and let τ = min{2β, (4c 5 |v 0 1 |) -1 }. Then for τ ∈ (0, τ ), q = (2β -τ )/(2β + τ ) ∈ (0, 1) and (4.18) is satisfied for N 0 = 0 since |v

1/2 1 | = |(1 + q)v 0 1 /2| ≤ |v 0 1 |. Moreover, τ 1 -q = β + τ 2 ≤ 2β.
By the energy estimate (4.17), C ′ n is bounded by a constant independent of n and τ . Since all norms are equivalent in V h , estimate (4.27) implies that for τ > 0 small enough (but dependent on h !), δu h = 0. Now, assume that the bound (4.24) is satisfied, and replace τ by

τ = min{2β, (4c 5 C 1 ) -1 }.
By the energy estimate (4.17), E h (U n , V n ) is bounded by a constant independent of h, n and τ . Thus, C ′ n is bounded by a constant C ′ independent of h and n. We apply Lemma 4.7 below with ε 1 = 1/(4(C ′ C S + 1)) and ε 2 = 1/(2c f ), and we obtain

2β τ 2 |δ uh | 2 -1,h + 1 τ |δ uh | 2 -1,h ≤ c f 2 1 4ε 2 2 + 1 4 + C ′ C S + 1 4ε 2 1 |δ uh | 2 -1,h .
We see that for τ > 0 small enough (independent of h now), δ uh = 0 and the proof is complete.

Lemma 4.7. Let ε 1 , ε 2 > 0. Then, for all uh ∈ Vh , there hold

| uh | 2 1 ≤ ε 1 Ȧh uh 2 + 1 4ε 2 1 | uh | 2 -1,h , (4.28) uh 2 ≤ ε 2 Ȧh uh 2 + 1 4ε 2 2 + 1 4 | uh | 2 -1,h . (4.29)
Proof. By arguing as in (3.12), we see that

| Ȧ1/2 U | = ( Ȧ U ) t U ≤ ε 1 2 | Ȧ U | 2 + 1 2ε 1 | U | 2 .
Let ε > 0. Similarly, we have

| U | 2 = |( Ȧ1/2 U ) t Ȧ-1/2 U | ≤ ε| Ȧ1/2 U | 2 + 1 4ε | Ȧ-1/2 U | 2 . (4.30)
Thus, we get

| Ȧ1/2 U | ≤ ε 1 2 | Ȧ U | 2 + 1 2ε 1 ε| Ȧ1/2 U | 2 + 1 4ε | Ȧ-1/2 U | 2 .
By choosing ε = ε 1 , we obtain (4.28). Next, we plug (4.28) into (4.30), with ε = 1 and ε 1 = ε 2 , and we deduce (4.29).

4.3.

Convergence as (h, τ ) → (0, 0). For a time step τ > 0, let (u n h , v n h , z n h , w n h ) n≥1 be a solution of the fully discrete scheme (4.2). We define the following functions from R + into V h :

u τ h (t) = ((n + 1) -t/τ )u n h + (t/τ -n)u n+1 h , t ∈ [nτ, (n + 1)τ ) (n ∈ N), u τ h (t) = u n+1 h , t ∈ [nτ, (n + 1)τ ) (n ∈ N), u τ h (t) = u n h , t ∈ [nτ, (n + 1)τ ) (n ∈ N), ûτ h (t) = (u n h + u n+1 h )/2, t ∈ [nτ, (n + 1)τ ) (n ∈ N).
We define similarly the functions v τ h , v τ h , v τ h , vτ h associated to the sequence (v n h ) n≥0 and the functions z τ h , w τ h . Notice that ûτ h = (u τ h + u τ h )/2 for all t ∈ R + and that

∂ t u τ h = (u n+1 h -u n h )/τ in D ′ (0, ∞); V h . The convergence results is as follows: Theorem 4.8. Let (u 0 , v 0 ) ∈ H 2 per × H -1 per . Assume that (u 0 h , v 0 h ) h>0 is a family of functions in V h × V h which satisfies assumptions (3.25)-(3.26) as h → 0. Then the solution (u τ h , v τ h
) associated to the fully discrete scheme (4.2) tends to the energy solution of problem (1.1)-(1.2) in the following sense, as (h, τ ) → (0, 0):

u τ h → u weakly ⋆ in L ∞ (R + ; H 1 per ), u τ h → u strongly in C([0, T ], L 2 (Ω)), for all T > 0, A h u τ h → Au weakly ⋆ in L ∞ (R + ; L 2 (Ω)), Ȧ-1 h ∂ t uτ h → Ȧ-1 ∂ t weakly ⋆ in L ∞ (R + ; Ḣ1 per )
and weakly in L 2 (R + ; Ḣ1 per ). Proof. We proceed as in the proof of Theorem 3.3. We first consider the conservation law for the mass. By choosing ϕ h = 1 and ψ h = 1 in (4.2), we find

(∂ t u τ h , 1) = vτ h β∂ t v τ h + vτ h = 0, (4.31)
in D ′ ((0, ∞)). The estimates below show that (u τ h ) h>0,τ >0 is bounded in L ∞ (R + ; H 1 per ), so that, up to a subsequence, u τ h → u in L ∞ (R + ; H 1 per ) weakly ⋆, and so (∂ t u τ h , 1) → (∂ t u, 1) in D ′ (0, ∞) , as (h, τ ) → (0, 0). By (4.12), |v n 1 | ≤ |v 0 1 | for all n. Thus, vτ h is bounded in L ∞ (R + ), and so, up to a subsequence, vτ h converges weakly ⋆ in L ∞ (R + ) to some function a ∈ L ∞ (R + ). Moreover, by (4.12), we have

v n+1 1 -v n 1 = |q| n |1 -q||v 0 1 | = |q| n 2τ 2β + τ |v 0 1 | ≤ τ β |v 0 1 |, Observe now that v τ h -vτ h = (t/τ -(n + 1/2))(v n+1 h -v n h ), t ∈ [nτ, (n + 1)τ ), (n ∈ N). (4.32) Therefore we get | v τ h -vτ h | ≤ τ |v 0 1 |/(2β) and so | v τ h -vτ h | converges uniformly to 0 in R + , as (h, τ ) → (0, 0). Hence v τ h converges to a weakly ⋆ in L ∞ (R + )
. We can pass to the limit in (4.31) in the sense of distributions on (0, ∞) and we find (∂ t u, 1) = a(t) β∂ t a(t) + a(t) = 0, which is the conservation law for the mass.

We now turn to the energy estimate. Let

τ ⋆ = min{2β, (4c 5 sup h>0 v 0 h ) -1 }.
If τ ≤ τ ⋆ , then (4.18) is satisfied for N 0 = 0 (and for all h > 0). Since τ 1 -q = β + τ /2 ≤ 2β, the energy estimate (4.17) implies

E h (u n h , v n h ) + n-1 k=0 τ |v k+1/2 h | 2 -1,h ≤ exp(32βc 5 |v 0 1 |) E h (u 0 h , v 0 h ) + 2βc 7 |v 0 1 | , (4.33)
for all n ≥ 0. Assumptions (3.25)- (3.26) imply that E h (u 0 h , v 0 h ) and |v 0 1 | are bounded by a constant independent of h. The right-hand side of (4.33) is bounded by a constant independent of h and τ . Thus,

u τ h is uniformly bounded in H 1 per , z τ h = A h ûτ h is uniformly bounded in L 2 (Ω), and
ṙτ h := Ȧ-1 h ∂ t uτ h = Ȧ-1 h vτ h is uniformly bounded in H 1 per .
By arguing as in (3.28), we see that for all 0 ≤ s ≤ t,

uτ h (t) -uτ h (s) 2 ≤ 4 ṙτ h L ∞ (R + ;H 1 per ) uτ h L ∞ (R + ;H 1 per ) |t -s|. Moreover, for all 0 ≤ s ≤ t, observe that | u τ h (t) -u τ h (s) | = t s vτ h (σ) dσ ≤ | v 0 h ||t -s|.
Thus, for all T > 0, there is a constant C T independent of (h, τ ) such that

u τ h (t) -u τ h (s) ≤ C T |t -s| 1/2 , (4.34)
for all 0 ≤ s ≤ t ≤ T . By the Ascoli-Arzelà Theorem, (u τ h ) is precompact in the space C([0, T ]; L 2 (Ω)), for all T > 0. Applying (4.34) with s = nτ and t = (n + 1)τ yields

u n+1 h -u n h ≤ C T τ 1/2 , so that u τ h -u τ h L ∞ (0,T ;L 2 (Ω)) → 0 and u τ h -u τ h L ∞ (0,T ;L 2 (Ω))
→ 0, as τ → 0. Up to a subsequence, we have

u τ h , ûτ h → u weakly ⋆ in L ∞ (R + ; H 1 per ), u τ h , u τ h → u strongly in C([0, T ]; L 2 (Ω)), for all T > 0, u τ h , u τ h → u a.e. in Ω × R + , z τ h → weakly ⋆ in L ∞ (R + ; L 2 (Ω)), ṙτ h → ṙ weakly ⋆ in L ∞ (R + ; Ḣ1 per ) and weakly in L 2 (R + ; Ḣ1 per ), as (h, τ ) → (0, 0). Let ϕ ∈ H 1 per and set ϕ h = Π h (ϕ), so that ϕ h → ϕ strongly in H 1 per . The first equation in (4.2) reads (∂ t u τ h , ϕ h ) = (v τ h , ϕ h )
. By arguing as in (3.29) and letting (h, τ ) → (0, 0), we obtain that

∂ t ( u, φ) = (∇ ṙ, ∇ φ) = v, φ H -1 per ,H 1 per in D ′ ((0, ∞)), with v = Ȧ ṙ. This shows that ∂ t u = v ∈ L ∞ (R + ; H 1 per ). Next, we set ψ ∈ H 2
per and we let ψ h = Π h (ψ) so that ψ h → ψ strongly in H 1 per . We have ψ h = N h i=1 Ψ i e i h and Ψ = (ψ 1 , . . . , ψ N h ) t is the vector associated to ψ h . On multiplying (4.5) by Ψt Ȧ-1 , we find

β(∂ t vτ h , ψh ) -1,h + ( vτ h , ψh ) -1,h + (∇z τ h , ∇ψ h ) -2(∇û τ h , ∇ψ h ) + 1 2 (f (u τ h ) + f (u τ h ), ψh ) - 1 12 ((u τ h -u τ h ) 2 (f ′′ + (u τ h ) + f ′′ -(u τ h )), ψh ) = 0. (4.35)
By arguing as in the proof of Theorem 3.3, we get

( vτ h , ψh ) -1,h → (v, ψ) -1 , (∇z τ h , ∇ψ h ) → (z, Aψ), (∇û τ h , ∇ψ h ) → (∇u, ∇ψ) in D ′ ((0, ∞)). Thanks to the Sobolev injection H 1 per ֒→ L 2p+2 (Ω)
, for all T > 0, the terms f (u τ h ) L q (0,T ;L q (Ω)) , f (u τ h ) L q (0,T ;L q (Ω)) , and

(u τ h -u τ h ) 2 (f ′′ + (u τ h ) + f ′′ -(u τ h
)) L q (0,T ;L q (Ω)) are bounded by a constant independent of h and τ , for q = (2p + 2)/(2p + 1) ∈ (1, 2). We can therefore pass to the limit in the nonlinear terms, and we find that

1 2 (f (u τ h ) + f (u τ h ), ψh ) → (f (u), ψ) and 1 12 ((u τ h -u τ h ) 2 (f ′′ + (u τ h ) + f ′′ -(u τ h )), ψh ) → 0 in D ′ ((0, ∞))
. Thus, the first term in equation (4.35) has a limit in D ′ ((0, ∞)),

β(∂ t vτ h , ψh ) -1,h → η ψ .
As a consequence,

( vτ h -vτ h , ψh ) -1,h = τ (∂ t vτ h , ψh ) -1,h → 0 in D ′ ((0, ∞))
. Thus, by (4.32), we have

( vτ h , ψh ) -1,h → ( v, ψ) -1 .
Summing up, we have proved that

β∂ t ( v, ψ) -1 + ( v, ψ) -1 + (z, Aψ) -2(∇u, ∇ψ) + (f (u), ψ) = 0,
with v = ∂ t u and z = Au. We conclude as in Theorem 3.3 that (u, u t ) is an energy solution of (1.1)-(1.2). Note that the whole family converges to (u, u t ) due to the uniqueness of the limit.

Convergence to equilibrium

In this section, we prove that any solution of the fully discrete scheme converges to a single equilibrium, for any time step τ > 0. The parameter h is fixed (so that assumption (H2) is not relevant). We adapt the proof from [START_REF] Grasselli | Well-posedness and longtime behavior for the modified phase-field crystal equation[END_REF] in a discrete setting. The main idea is to use the gradient-like flow structure of the problem and a suitable Lojasiewicz inequality. In three space dimensions, in addition to (H1) and (H3), we need the following assumption:

(H4) If d = 3, then either V h ⊂ L ∞ (Ω) or p = 1.

Theorem 5.1. τ > 0 denote the time step and let

(U n , V n ) n≥0 denote any sequence in R N h × R N h which complies with (4.5). Then (U n , V n ) converges to (U ∞ , 0), where U ∞ = (u ∞ 1 , U ∞ ) is a stationary solution with average constraint, i.e., u ∞ 1 = M = u 0 1 + βv 0 1 , Ȧ2 U ∞ -2 Ȧ U ∞ + ∇F h (U ∞ ) = 0.
(5.1)

We first prove the following Lemma 5.2. Let the assumptions of Theorem 5.1 hold. Then V n+1/2 → 0.

Proof. Since |q| < 1, estimate (4.18) is satisfied for N 0 large enough. By the energy estimate (4.17),

∞ n=N 0 | V n+1/2 | 2 -1 < ∞. In particular, V n+1/2 → 0 in R N h -1 . Moreover, v n+1/2 1 = q n v 1/2 1
by (4.12), so V n+1/2 → 0, as claimed.

For any M ∈ R, we introduce the auxiliary function F M (y) = F (M + y) and the following functionals

F M,h ( U ) = (F M ( uh ), 1), (5.2)
E M,h ( U ) = 1 2 | Ȧ U | 2 -| Ȧ1/2 U | 2 + F M,h ( U ), (5.3)
E M,h ( U , V ) = E M,h ( U ) + β 2 | V | 2 -1 , (5.4)
defined for every U ≃ uh and every V in R N h -1 .

For any M ∈ R, we also consider

S M = U ∈ R N h : U satisfies (5.1) . For any sequence (U n , V n ) n≥0 in R N h × R N h , we define its ω-limit set in R N h × R N h : ω ((U n , V n ) n≥0 ) = {(U ⋆ , V ⋆ ) : ∃n j ր ∞, (U n j , V n j ) → (U ⋆ , V ⋆ )|} .
Similarly, we set

ω (U n ) n≥0 = {U ⋆ : ∃n j ր ∞, U n j → U ⋆ } .
We have:

Proposition 5.3. Let the assumptions of Theorem 5.1 hold. Then ω ((U n , V n ) n≥0 ) is a nonempty compact and connected set such that

ω ((U n , V n ) n≥0 ) = ω (U n ) n≥0 × {0} ⊂ {(U ⋆ , 0) : U ⋆ ∈ S M } , with M = u 0 1 + βv 0 1 . Moreover, E M,h is constant on ω (U n ) n .
This result implies in particular that V n → 0, as proved below.

Proof. Since q = (2β -τ )/(2β + τ ), we can rewrite (4.13) as

u n 1 = u 0 1 + (1 -q n )βv 0 1 . (5.5)
Let M = u 0 1 + βv 0 1 . We introduce the auxiliary functions f M (y) = f (M + y) and fM (r, s) = f (M + r, M + s),

where f (r, s) = 1 2 (f (r) + f (s)) - 1 12 (s -r) 2 (f ′′ + (r) + f ′′ -(s)
). We also set Then we rewrite the second equation in (4.5) in the following form:

β( V n+1 -V n )/τ + V n+1/2 + Ȧ Ȧ2 U n+1/2 -2 Ȧ U n+1/2 + JM ( U n , U n+1 ) = Ȧ JM ( U n , U n+1 ) -J(U n , U n+1 ) , (5.6)
where

JM ( U n , U n+1 ) = ( fM ( un h , un+1 h ), e i ) 2≤i≤N h , J(U n , U n+1 ) = ( f (u n h , u n+1 h ), e i ) 2≤i≤N h . Multiplying (5.6) by ( U n+1 -U n ) t Ȧ-1 , using that (4.5) implies U n+1 -U n = τ V n+1/2 , (5.7)
we find

β 2 (| V n+1 | 2 -1 -| V n | 2 -1 ) + τ | V n+1/2 | 2 -1 + 1 2 (| Ȧ U n+1 | 2 -| Ȧ U n | 2 ) -(| Ȧ1/2 U n+1 | 2 -| Ȧ1/2 U n | 2 ) + JM ( U n , U n+1 ), U n+1 -U n = ( fM ( un h , un+1 h ) -f (u n h , u n+1 h ), un+1 h -un h ).
(5.8)

Using now Lemma 4.4 and arguing as in the proof of Lemma 4.5, we obtain

(F M ( un+1 h ), 1) -(F M ( un h ), 1) ≤ JM ( U n , U n+1 ), U n+1 -U n .
(5.9) By (5.5), for any solution u n h ≃ U n of (4.5), we have

f (u n h , u n+1 h ) = fM ( un h -βq n v 0 1 , un+1 h -βq n+1 v 0 1 ). Thus we get ( fM ( un h , un+1 h ) -f (u n h , u n+1 h ), un+1 h -un h ) = ( fM ( un h , un+1 h ) -fM ( un h -βq n v 0 1 , un+1 h -βq n+1 v 0 1 ), un+1 h -un h ) = -βq n v 0 1 1 0 ∂ r fM (ũ n h (s),
) -f (u n h , u n+1 h ), un+1 h -un h ) ≤ β|q| n |v 0 1 |C( u n h 1 , u n+1 h 1 ) un+1 h -un h 1
(5.10)

≤ 1 4τ | U n+1 -U n | 2 -1 + C 0 |q| 2n . (5.11)
Here and in the following, C k (k = 0, 1, . . . ) denotes a constant independent of n (but which may depend on τ , h and other parameters of the problem). In the last inequality we have used that all norms are equivalent in Vh . Adding up (5.8), (5.9) and (5.11), we find

E M,h ( U n+1 , V n+1 ) -E M,h ( U n , V n ) + 3τ 4 | V n+1/2 | 2 -1 ≤ C 0 |q| 2n , (5.12)
for all n ≥ 0. Set now

G n = Ȧ-1 V n , Ȧ-1 ( Ȧ2 U n -2 Ȧ U n + ∇F M,h ( U n ) -1 , where U , V -1 = U t Ȧ-1 V , for all U , V ∈ R N h -1 . We have G n+1 -G n = Ȧ-1 ( V n+1 -V n ), Ȧ-1 ( Ȧ2 U n -2 Ȧ U n + ∇F M,h ( U n ) -1 + Ȧ-1 V n+1 , Ȧ-1 (( Ȧ2 -2 Ȧ)( U n+1 -U n ) + ∇F M,h ( U n+1 ) -∇F M,h ( U n ) -1 .
Let δG n 1 denote the first term on the right-hand side of this equality and denote the second by δG n 2 . Using (5.6),

U n = U n+1/2 -( U n+1 -U n )/2 and JM ( U n , U n ) = ∇F M,h ( U n ), we obtain δG n 1 = - τ β Ȧ-1 V n+1/2 + Ṡn -JM ( U n , U n+1 ) -J(U n , U n+1 ) , Ȧ-1 Ṡn -1 τ 2β Ȧ-1 V n+1/2 + Ṡn -JM ( U n , U n+1 ) -J(U n , U n+1 ) , Ȧ-1 Ṫ n 1 -1 , where Ṡn = Ȧ2 U n+1/2 -2 Ȧ U n+1/2 + JM ( U n , U n+1 ), Ṫ n 1 = ( Ȧ2 -2 Ȧ)( U n+1 -U n ) + 2 JM ( U n , U n+1 ) -2 JM ( U n , U n ).
Using that all norms are equivalent in Vh , the Cauchy-Schwarz inequality, Young's inequality and (5.7), we deduce

δG n 1 + 3τ 4β | Ȧ-1 Ṡn | 2 ≤ C 1 τ | V n+1/2 | 2 -1 + C 2 JM ( U n , U n+1 ) -J(U n , U n+1 ) 2 +C 3 JM ( U n , U n+1 ) -JM ( U n , U n ) 2 .
By Bessel's inequality, we get

JM ( U n , U n+1 ) -J(U n , U n+1 ) 2 ≤ fM (u n h , u n+1 h ) -f (u n h , u n+1 h ) 2 .
Arguing as in (5.10), and using assumption (H4), we find that

JM ( U n , U n+1 ) -J(U n , U n+1 ) 2 ≤ C 4 |q| 2n .
(5.13)

Similarly, we have

JM ( U n , U n+1 ) -JM ( U n , U n ) 2 ≤ C 5 τ V n+1/2 2 -1
.

(5.14)

Summing up, we have proved

δG n 1 + 3τ 4β | Ȧ-1 Ṡn | 2 ≤ (C 1 + C 3 C 5 )τ | V n+1/2 | 2 -1 + C 2 C 4 |q| 2n ,
for all n ≥ 0. We now consider the term δG n 2 . Using V n+1 = V n+1/2 +(V n+1 -V n )/2, equation (5.6), and arguing as for δG n 1 , we obtain

δG n 2 ≤ τ 4β | Ȧ-1 Ṡn | 2 + C 6 τ | V n+1/2 | 2 -1 + C 7 |q| 2n .
Thus, we get

G n+1 -G n + τ 2β | Ȧ-1 Ṡn | 2 ≤ C 8 τ | V n+1/2 | 2 -1 + C 9 |q| 2n , (5.15)
for all n ≥ 0, with

C 8 = C 1 + C 3 C 5 + C 6 and C 9 = C 2 C 4 + C 7 .
Let us introduce the sequence

W n = 2E M,h ( U n , V n ) + νG n ,
where ν > 0 is sufficiently small so that νC 8 ≤ 1/2. From estimates (5.12) and (5.15), it follows that

W n+1 -W n + τ | V n+1/2 | 2 -1 + ντ 2β | Ȧ-1 Ṡn | 2 ≤ C 10 |q| 2n , (5.16)
with C 10 = 2C 0 + νC 9 . By the energy estimate (4.17), the sequence (U n , V n ) is bounded, so (W n ) n≥0 is bounded. This implies that W n converges to some real number W ∞ as n tends to ∞. Indeed, let

W n = W n + C 10 1 -|q| 2 |q| 2n .
Using (5.16), we see that W n+1 -W n ≤ 0, i.e. W n is nonincreasing. and since W is bounded, W n has a limit W ∞ = W ∞ .

Adding up estimate (5.16), we obtain that ∞ n=0 | Ȧ-1 Ṡn | 2 < ∞. In particular, Ṡn → 0. Moreover, by (5.13), we have

JM ( U n , U n+1 ) -J(U n , U n+1 ) → 0 From (5.16), it follows that ∞ k=n a 2 k ≤ W n -W ∞ + C 11 |q| 2n .
On the other hand, using the Lojasiewicz inequality (5.18) and the fact 1/(1-θ) < 2, we deduce that, for all n ≥ n 0 (changing n 0 into a larger integer if necessary),

|W n -W ∞ | ≤ 2|E M,h ( U n ) -E ∞ | + β| V n | 2 -1 + ν|G n | ≤ 2| Ȧ2 U n -2 Ȧ U n + ∇F M,h ( U n )| 1/(1-θ) + β| V n | 2 -1 +C 12 | V n | -1 | Ȧ2 U n -2 Ȧ U n + ∇F M,h ( U n )| ≤ C 13 | Ȧ2 U n -2 Ȧ U n + ∇F M,h ( U n )| 1/(1-θ) + | V n | 1/(1-θ) -1 . (5.20) Using V n = V n+1/2 -( V n+1 -V n )/2,
we deduce from (5.6) and (5.13) that

| V n | -1 ≤ C 13 | V n+1/2 | -1 + | Ȧ-1 Ṡn | + |q| n .
(5.21)

Similarly, from Ȧ2 U n -2 Ȧ U n + ∇F M,h ( U n ) = Ṡn + ( Ȧ2 -2 Ȧ)( U n -U n+1/2 ) + JM ( U n , U n ) -JM ( U n , U n+1 ), U n+1 -U n = τ V n+1/2 and (5.14), we infer that | Ȧ2 U n -2 Ȧ U n + ∇F M,h ( U n )| ≤ C 14 | Ȧ-1 Ṡn | + | V n+1/2 | -1 . (5.22)
Collecting (5.20), (5.21) and (5.22), we obtain

|W n -W ∞ | ≤ C 15 | V n+1/2 | 1/(1-θ) -1 + | Ȧ-1 Ṡn | 1/(1-θ) + (|q| n ) 1/(1-θ) ,
for all n ≥ n 0 . This gives

∞ k=n a 2 k ≤ C 16 a 1/(1-θ) n , ∀n ≥ n 0 .
From Lemma 5.5 below, we conclude that ∞ n=0 a n < ∞. In particular, we have

τ ∞ n=0 | V n+1/2 | -1 = ∞ n=0 | U n+1 -U n | < ∞.
This shows that the whole sequence ( U n ) n has a limit U ∞ as n → ∞. From (5.5), we know that u n 1 → M . Thus, (U n ) n tends to some U ∞ in R N h , and the proof is complete.

For the following lemma and its proof, we adapt Lemma 4.1 in [START_REF] Huang | Convergence in gradient-like systems which are asymptotically autonomous and analytic[END_REF] in a discrete setting (see also Lemma 7.1 in [START_REF] Feireisl | Convergence for semilinear degenerate parabolic equations in several space dimensions[END_REF]). Lemma 5.5. Let 0 < θ < 1/2. Assume that (a n ) n≥0 is a sequence of nonnegative real numbers such that ∞ n=0 a 2 n < ∞, and there are a constant C > 0 and an integer n 0 such that

∞ k=n a 2 k ≤ Ca 1/(1-θ)
n for all n ≥ n 0 .

(5.23)

Then ∞ n=0 a n < ∞.
Proof. First replacing a n by max{a n , 1} for 0 ≤ n < n 0 , and then taking C large enough to ensure C ≥ ∞ n=0 a 2 n , we observe that (5.23) becomes valid for all n ≥ 0. So we may assume n 0 = 0. Set now

ρ n := ∞ k=n a 2 k and σ n = n k=0 a k for n ≥ 0.
Given any n ≥ 0, we first raise inequality (5.23) to the power 1 -θ > 0:

ρ 1-θ n ≤ C 1-θ a n .
Next, we sum this relation and we obtain

n k=0 ρ 1-θ k ≤ C 1 n k=0 a k = C 1 σ n .
We now apply a discrete integration-by-parts on the left-hand side

n k=0 [(k + 2) -(k + 1)]ρ 1-θ k = (n + 2)ρ 1-θ n -ρ 1-θ 0 + n k=1 (k + 1)(ρ 1-θ k-1 -ρ 1-θ k ).
Next, we notice that

ρ 1-θ k-1 -ρ 1-θ k = ρ k-1 ρ k (1 -θ)s -θ ds ≥ (1 -θ)a 2 k-1 ρ -θ k-1 , since ρ k-1 = ρ k + a 2 k-1 . This gives (n + 2)ρ 1-θ n -ρ 1-θ 0 + (1 -θ) n k=1 (k + 1)a 2 k-1 ρ -θ k-1 ≤ C 1 σ n .
It follows that, for every n ≥ 0,

(n + 1)ρ 1-θ n ≤ C 2 (1 + σ n ),
and

n k=1 (k + 1)a 2 k-1 ρ -θ k-1 ≤ C 2 (1 + σ n ),
where C 2 > 0 is a constant independent of n. Since the sequence (σ n ) is nondecreasing, the former estimate yields

ρ k-1 ≤ C 3 (1 + σ n ) 1/(1-θ) k -1/(1-θ) , 1 ≤ k ≤ n < ∞,
which we insert into the latter one, thus arriving at

n k=1 k 1+θ/(1-θ) a 2 k-1 ≤ C 4 (1 + σ n ) 1+θ/(1-θ) .
The constants C 3 , C 4 are independent of n ≥ 0. As a consequence, using the Cauchy-Schwarz inequality, we obtain for n ≥ 1,

σ n-1 = n k=1 a k-1 ≤ n k=1 k 1/(1-θ) a 2 k-1 1/2 n k=1 k -1/(1-θ) 1/2 ≤ C 5 (1 + σ n ) 1/2(1-θ) .
We conclude that the sequence (σ n ) n must be bounded, since 2(1 -θ) > 1. Indeed, assume by contradiction that (σ n ) n is unbounded, and set

r n = 1 + σ n ≥ 1. Then (r n ) is nondecreasing, r n → ∞ and r n-1 ≤ C 6 r 1/2(1-θ)
n , so that r n-1 /r n → 0, and we deduce that r n → 0. Contradiction. Thus (σ n ) is bounded, i.e.

∞ n=0 a n < ∞, as claimed. Remark 5.6. In the proof of convergence to equilibrium, we have used the fact that all norms are equivalent in V h . An interesting open question would be to prove a similar result for the time semi-discrete version of our problem.

Remark 5.7. By arguing as in the continuous case (see [START_REF] Grasselli | Well-posedness and longtime behavior for the modified phase-field crystal equation[END_REF]), using the energy estimate (cf. Lemma 3.1) and the Lojasiewicz inequality 5.17, it is possible to prove that any solution (u h , v h ) of the space semi-discrete scheme (3.1) converges to a single equilibrium, provided assumptions (H1) and (H4) hold.

Numerical results

We present some numerical results in one space dimension (obtained with the Scilab software 1 and in two space dimensions (obtained with the Freefem++ software [START_REF] Hecht | New development in FreeFem++[END_REF]). In every case, the nonlinearity f is given by (1.3) for some parameter ε, and we set f + = f , f -= 0 in assumption (H3). The space V h is the space of piecewise linear (P 1 ) finite elements.

6.1. Simulations in one space dimension. We first choose an interval Ω = (0, L) with L = 4π. In Table 1, we compute the error in the C 0 ([0, T ]; L 2 (Ω))-norm (which appears in Theorem 4.8). The parameters are ε = 0.5, β = 0.5 and T = 2. We use a uniform grid with space stepsize h = L/M and time stepsize τ = T /N . The initial value (u 0 h , v 0 h ) is the P 1 -interpolate of u 0 (x) = cos(x) + 0.3 cos(3x), v 0 (x) = 0.1. For the error of the time discretization, h = L/160 is fixed. Since the exact solution u h of the space semi-discrete scheme (3.1)-(3.2) is unknown, we use instead the solution on a fine grid with stepsize τ sol = T /5120. evaluated on the fine grid t k = kτ sol (k = 0, 1, . . . , 5120), and the ratio err h (τ )/err h (τ /2)

1 Scilab is freely available at http://www.scilab.org/ of consecutive errors. The computed ratio is close to 4, which means that the convergence error for the time discretization is O(τ 2 ), as expected.

For the error of the space discretization, τ = T /160 is fixed. We use again the solution u τ h sol on a fine grid with stepsize h = L/2560 for the comparison. Figure 1 shows the plot of the pseudo-energy E h (u n h , v n h ) (see (3.22)) versus the time t (in solid line). The domain is Ω = (0, L) with L = 4π. The parameters are ε = 0.5, β = 5, h = L/320 and dt = 0.005. The initial condition is the P 1interpolate of u 0 (x) = 0.1/(1 + 0.7 cos(x)), v 0 = 0.034. The black color corresponds to the second-order scheme (4.2); the blue color corresponds to a first-order scheme obtained by applying to the space semi-discrete scheme (3.1) the time discretization proposed by Wang and Wise [START_REF] Wang | Global smooth solutions of the three-dimensional modified phase field crystal equation[END_REF][START_REF] Wang | An energy stable and convergent finite-difference scheme for the modified phase field crystal equation[END_REF]. Both schemes are unconditionally stable.

The left figure shows the pseudo-energy on the interval [0, 40]. If we had v 0 = 0, then by (4.16), the pseudo-energy would be nonincreasing in both cases. Here, the pseudo-energy exhibits oscillations due to the fact that v 0 = 0. In both cases, the evolution is driven to a stationary state, as predicted by the theory (see Theorem 5.1). We notice that the first-order scheme has a smoothing effect which creates more dissipation, especially at the beginning of the evolution. This is seen in the right figure which shows the energy E h (u n h ) on the interval [0, 2.5] (in dashed and dashdot), in addition to the pseudo-energy E h (u h , v h ) (in solid line). The difference E h (u n h , v n h )-E h (u n h ) = (β/2)| vn h | 2 -1,h can be interpreted as a "kinetic energy". 6.2. Simulations in two space dimensions. The domain Ω is the square (0, 6π)× (0, 6π). It is decomposed in 50×50 squares, and each square is divided along the lower left/upper right diagonal, resulting in a uniform triangulation of ω. The parameters are β = 0.1, ε = 2 and the time step is τ = 0.25. The initial condition is the P 1interpolate of u 0 (x, y) = 0.2 + 0.2 cos(x) cos(y) and v 0 = 0. Figures 2 and3 show the evolution from stripes to a triangular distribution of drops. Numerical tests up to time t = 1250 indicate that the triangular distribution of drops is the steady state for this simulation.

For the continuous problem (1.1), using the translation invariance, from a triangular distribution of drops we easily build a two dimensional continuum of steady states. For the fully discrete scheme (4.2), the translation invariance is broken by the space discretization, but we expect a large number of steady states. This simulation illustrate the convergence to equilibrium result (Theorem 5.1) in a situation where the steady state is not unique.

  the inner product as well as the norm on L 2 (Ω) are simply indicated by (•, •) and • , respectively. For sake of simplicity, we assume that Ω 1dx = |Ω| = 1. The mean value of any function u ∈ L 2 (Ω) is denoted by u = Ω udx, and we set u = u -u . The dual space of H m per is denoted by H -m per , and it is equipped with the operator norm given by T -m = sup u m=1, u∈H m per |T (u)|. For an operator u ∈ H -m per , we denote u = u, 1 H -m per ,H m per and we set u = u -u .

  -≤ 0, and deg(F -) < deg(F ) (here, deg denotes the degree of the polynomial).

Theorem 4 . 3 (

 43 Existence for any τ ). For any

Lemma 4 . 4 .

 44 Let g ∈ C 3 ([0, 1]; R). Then the following identities hold

1 = u n+1 h

 1n+1 is uniquely determined (see (4.9)). It is sufficient to show that un+1 h ≃ U n+1 is uniquely determined by (4.10), for τ sufficiently small and independent of n. Then v n+1 h can be recovered by the first equation in(4.5).Assume that (4.10) has two solutions un+1 h ≃ U n+1 and un+1 h ≃ U n+1 . We subtract the two resulting systems (4.10), and we multiply by δ

F

  M,± (y) = F ± (M + y) and f M,± (y) = f ± (M + y),so that F M = F M,+ + F M,-and F (iv) M,+ ≥ 0, F (iv)M,-≤ 0. In particular, the function F M satisfies the decomposition (H3), and we havefM (r, s) = 1 2 (f M (r) + f M (s)) -1 12 (s -r) 2 (f ′′ M,+ (r) + f ′′ M,-(s)).

  k ) -u τ sol h (t k ) L 2 (0,L)

Figure 1 .

 1 Figure 1. (Pseudo-)Energy vs time

Figure 2 .Figure 3 .

 23 Figure 2. t = 5 (left) and t = 26.25 (right)

Table 1 .

 1 Convergence error for the time (left) and for the space (right) discretization

	N = T /τ	err h (τ ) ratio	M = L/h	err τ (h) ratio
	80	0.5018208	-	40	0.7770682	-
	160	0.1455507 3.45	80	0.2735932 2.84
	320	0.0368516 3.95	160	0.0706798 3.87
	640	0.0091325 4.04	320	0.0175677 4.02
	1280	0.0022163 4.12	640	0.0041882 4.20

Table 1 (

 1 left) shows the errorerr h (τ ) = max

Table (

 ( (t k ) -u τ h sol (t k ) L 2 (0,L)evaulated on the grid t k = kτ (k = 0, 1, . . . , 160) and the ratio err τ (h)/err τ (h/2) of consecutive errors. Again, the computed ratio is close to 4: the convergence error for the space discretization is O(h 2 ), as expected.

												1)
	(right) shows the error
	E	0.16 0.18 0.2 0.22 0.26 0.24	err τ (h) = max 0≤k≤160 20 25 30 35 40 2nd order 1st order h 0 u τ 5 10 15	E	0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.24 0.26 0.22	0	0.5	1	1.5	2	2nd order 2nd order 1st order 1st order	2.5
			t					t

Acknowledgments

The first author thanks the Laboratoire de Mathématiques et Applications for the kind hospitality and the Université de Poitiers which supported his stay in May 2014. The first author is member of the Gruppo Nazionale per l'Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA) of the Istituto Nazionale di Alta Matematica (INdAM). The second author thanks the Dipartimento di Matematica of the Politecnico di Milano for the invitation in February 2014.

as n → ∞. From (5.6) and Lemma 5.2, it follows that V n+1 -V n → 0 and that

Using the conservation law (5.5), we see that U ⋆ ∈ S M . Finally, the sequence (U n , V n ) is bounded, and we have seen that U n+1 -U n → 0, V n → 0 so the ω-limit set ω (U n , V n ) n≥0 is a nonempty compact and connected subset of

We notice that the functional E M,h is a polynomial of the variables (u 2 , . . . , u N h ) of total degree 2p + 2, so the following Lojasiewciz inequality holds: Lemma 5.4 ( Lojasiewicz' inequality [START_REF] Lojasiewicz | Ensembles semi-analytiques[END_REF]). Let U ⋆ ∈ R N h -1 be a critical point of E M,h . Then there exist constants θ ∈ (0, 1/2) and δ > 0 such that for any

(5.17)

Proof of Theorem 5.1. Let M = u 0 1 + βv 0 1 as previously. By Lemma 5.4, for every U ∞ ∈ ω (U n ) n , there exist some δ > 0 and θ ∈ (0, 1/2) that may depend on U ∞ such that the inequality (5.17) holds for all U in

..,m , where the constants δ i , θ i corresponding to U ∞ i in Lemma 5.4 are indexed by i. From the definition of ω ( U n ) n , we know that there exists a sufficiently large n (5.18) where E ∞ = W ∞ /2 is the constant value of E M,h on ω ( U n ) n .

Let us now set .19)


Machine Learning (ML) approaches are increasingly used to model data coming from sensor networks. Typical ML implementations are cpu intensive and are often running server-side. However, IoT devices provide increasing cpu capabilities and some classes of ML algorithms are compatible with distribution and downward scalability. In this demonstration we explore the possibility of distributing ML tasks to IoT devices in the sensor network. We demonstrate a concrete scenario of appliance recognition where a smart plug provides electrical measures that are distributed to WiFi nodes running the ML algorithms. Each node estimates class-conditional probabilities that are then merged for recognizing the appliance category. Finally, our architectures relies on Web technologies for complying with Web-of-Things paradigms.

I. INTRODUCTION

Nowadays, many objects are following the Internet-of-Things (IoT) paradigm by relying on well known IP technologies. In spite of increasing capabilities, those devices are still performing simple tasks such as sense and react. We also observe the emergence of data-driven ML algorithms that are able to leverage on the large amount of sensor data. ML algorithms are usually executed on powerful computers often residing server-side. In this demonstration, we explore the feasibility of deporting the execution of those algorithms directly on things. For this purpose, we use several OpenPicus FlyportPRO WiFi [START_REF]Openpicus flyportpro wifi[END_REF] modules which will act as distributed computing power. Our implementation is also relying on Web technologies as application layer according to Web-of-Things (WoT) paradigms [START_REF] Guinard | A Web of Things Application Architecture -Integrating the Real World into the Web[END_REF].

II. MACHINE LEARNING APPROACH

The experimental setup demonstrates the feasibility and effectiveness of our approach by its deployment in the field of the electrical appliance recognition. The electrical consumption is measured with smart plugs placed between the appliance plug and the wall socket. Different features can be measured by the smart plug, as the active power, reactive power, RMS current, voltage, phase and frequency. The aim of the electrical appliance recognition task consists in analysing a sequence of electrical appliance measures, called appliance signatures, and provide information about the appliance generating the signal. Potentially, several information can be inferred, as the brand and model of the appliance, its category (i.e. coffee machine) or its actual state (i.e. stand-by). In this work, we demonstrate the recognition of the appliance category and state using an implementation inspired from [START_REF] Ridi | Appliance and State Recognition using Hidden Markov Models[END_REF].

For training the models we use the ACS-F2 database [START_REF] Ridi | ACS-F2 -A New Database of Appliance Consumption Analysis[END_REF], containing 450 electrical signatures recorded from 225 appliances of different brands and / or models. A signature consists of one hour recording with a sampling frequency of 10 -1 Hz. The appliances are uniformly spread into 15 categories. The signatures in the ACS-F2 database are represented in a six dimensional space as described above. Given the small quantity of memory available on the OpenPicus, we decide to reduce as much as possible the feature space, removing the non relevant frequency and voltage features, as well as the phase angle which is redundant with the active and reactive power information [START_REF] Ridi | Appliance and State Recognition using Hidden Markov Models[END_REF]. We selected Hidden Markov Models (HMMs) to represent the state-based nature of the signals. In our setting, we assume that the HMMs are trained offline due to the large quantity of data requested by ML training procedures. The training of the HMMs is performed here using a classical Expectation Maximisation procedure that iteratively estimates the parameters of the HMMs, i.e. the transition probabilities between states and the emission probability estimators, in our case Gaussian Mixtures Models. We consider here two states models: on and off (or stand-by). After training, the models parameters are transferred to the OpenPicus that are then in charge of processing real-time data. The choice of generative models such as HMMs is also motivated by our need to spread the computation load in the network of OpenPicus. This is naturally done considering that we can train one model M i per category and spread the computation of class conditional probabilities p(x|M i ) (likelihoods) on each nodes. The decision about the category label is finally taken using the simple Bayes rule, i.e. electing the highest likelihood if all categories have equal priors. In a previous work the system reported increasing accuracy using larger number of Gaussian per state, typically up to 30 [START_REF] Ridi | ACS-F2 -A New Database of Appliance Consumption Analysis[END_REF]. Given the memory restriction of the devices we are here limiting the number of Gaussians to 8, trading off slightly the accuracy rate.

III. GENERAL ARCHITECTURE

As our architecture relies on the WoT paradigm, every device embeds a RESTful API accessible through CoAP (i.e a lightweight version of HTTP) [START_REF] Shelby | Constrained application protocol (coap)[END_REF]. In order to sense nonphysical measures like an appliance category, we introduce the concept of virtual sensors. Those type of sensors are associated to complex data-driven process composed of several models called virtual classes. The Virtual sensor exploits the computational power of nearly located smart things and deploys models on them, forming a set of virtual classes. Virtual classes are agents able to compute class conditional probabilities (likelihoods) as explained above. They are periodically executing the algorithm using the data coming from the sensors. From a client's point of view, a virtual sensor acts like any other type of sensing device offering a Web resource representing the actual state. With a single GET request, a client can retrieve the actual value inferred from the Virtual classes.

Offline training algorithms will deploy the runtime algorithms by creating new virtual sensors resources. This step is achieved by providing a JSON document containing all the parameters of the machine learning models. The newly created virtual sensor then handles all the distribution of the models on the OpenPicus modules, creating virtual classes. Due to the memory limitation, parameters are exchanged in float binary format (IEEE 754) in order to avoid the JSON verbosity. Virtual classes will further register at physical sensors to be notified with new values. Our implementation offers two access modes for clients (selectable when creating the virtual sensor). In the end-toend mode, clients are requesting the virtual sensor for the actual appliance category. In order to determine the appliance category, each virtual class will be requested sequentially for its current likelihood. The virtual sensor then compares the likelihood values and, according to the Bayes rule, the model having the highest likelihood is designated as winner (assuming equal priors). The winning appliance category and its actual state are then returned to the client. The sync-based mode is fully event-driven and is especially tailored for realtime systems. Virtual classes will notify the Virtual sensor each time a new likelihood is computed. The decision making is performed as soon as a likelihood changes and is cached on the virtual client. This allows to avoid the sequential requesting of virtual classes and improves the scalability as well as the round-trip time. Clients can optionally observe a virtual sensor that will send notifications as soon as the appliance category changes.

IV. DEMONSTRATION

Our demonstration is composed of a smart plug Plogg, a set of OpenPicus FlyportPRO WiFi nodes, a Raspberry Pi and a client tablet device. The OpenPicus nodes execute a model (virtual class resource), corresponding to a specific appliance category. The Raspberry Pi acts as virtual sensor for making the decision about the category of the appliance plugged in the smart plug. The virtual classes will regularly be notified about new values coming from the smart plug. We use a mobile device as client to feedback the user about information on the identified appliance. A web page using Sync-based mode is used and notified by the Virtual sensor in real-time. Figure1 illustrates the different entities and their data exchange. Offline experiments of the demonstration settings show correct identification rates of about 90%.

V. CONCLUSION

We have presented an architecture for executing the runtime part of a machine-learning process on smart things. An experimental setup demonstrates the feasibility of this architecture able to execute higher-level tasks on a network of IoT things. We successfully deployed complex trained models performing a classification task for appliance recognition. The extensive use of Web technologies among with the concepts of virtual sensor and virtual class allows a seamless integration in the world of IoT and WoT. Future works will include the recognition of user activities in houses, as well as the reconsidering the training part according to the WoT precepts.

VI. ACKNOWLEDGMENT

Fig. 1 .

 1 Fig. 1. Data exchange between the entities.

This work is supported by the grant Smart Living Green-Mod from the Hasler Foundation in Switzerland and by the HES-SO.