Erwan Faou

AND FR ÉD Romain Horsin

Éric Rousset

ON NUMERICAL LANDAU DAMPING FOR SPLITTING METHODS APPLIED TO THE VLASOV-HMF MODEL

Keywords: 1991 Mathematics Subject Classification. 35Q83, 35P25 Geometric Numerical Integration, Splitting methods, Vlasov equations, Landau damping, HMF model

We consider time discretizations of the Vlasov-HMF (Hamiltonian Mean-Field) equation based on splitting methods between the linear and non-linear parts. We consider solutions starting in a small Sobolev neighborhood of a spatially homogeneous state satisfying a linearized stability criterion (Penrose criterion). We prove that the numerical solutions exhibit a scattering behavior to a modified state, which implies a nonlinear Landau damping effect with polynomial rate of damping. Moreover, we prove that the modified state is close to the continuous one and provide error estimates with respect to the time stepsize.

INTRODUCTION

In this paper we consider time discretizations of the Vlasov-HMF model. This model has received much interest in the physics and mathematics litterature for many reasons: It is a simple ideal toy model that keeps several features of the long range interactions, it is a simplification of physical systems like charged or gravitational sheet models and it is rather easy to make numerical simulations on it. We refer for example to [START_REF] Barré | The Vlasov equation and the Hamiltonian Mean-Field model[END_REF], [START_REF] Yamaguchi | Stability criteria of the Vlasov equation and quasi stationary states of the HMF model[END_REF], [START_REF] Barré | Algebraic damping in the one-dimensional Vlasov equation[END_REF], [START_REF] Caglioti | Long time estimates in the mean field limit[END_REF], [START_REF] Caglioti | Quasi-stationary states for particle systems in the mean-field limit[END_REF] for more details. This model also has strong analogy with the Kuramoto model of coupled oscillators in its continuous limit [START_REF] Fernandez | Landau damping in the Kuramoto model[END_REF], [START_REF] Benedetto | Exponential dephasing of oscillators in the Kinetic Kuramoto Model[END_REF], [START_REF] Dietert | Stability and bifurcation for the Kuramoto model[END_REF]. A long time analysis of the Vlasov-HMF model around homogenous stationary states has been recently performed in [START_REF] Faou | Landau Damping In Sobolev Spaces For The Vlasov-HMF Model[END_REF] where a Landau damping result is proved in Sobolev regularity. The purpose of the present paper is in essence to show that this result persists through time discretization by splitting methods.

The Vlasov-HMF model reads

(1.1) ∂ t f (t, x, v) + v∂ x f (t, x, v) = ∂ x R×T P (x -y)f (t, y, u)dudy ∂ v f (t, x, v),
where (x, v) ∈ T × R and the kernel P (x) is given by P (x) = cos(x). We consider initial data under the form f 0 (x, v) = η(v) + εr 0 (x, v) where ε is a small parameter and r 0 is of size one (in a suitable functional space). This means that we study small perturbations of a stationary solution η(v). Writing the exact solution as

f (t, x, v) = η(v) + εr(t, x, v),
and setting

(1.2) g(t, x, v) = r(t, x + tv, v), the main result given in [START_REF] Faou | Landau Damping In Sobolev Spaces For The Vlasov-HMF Model[END_REF] is that if ε is small enough, g(t, x, v) converges towards some g ∞ (x, v) when t goes to ∞ in Sobolev regularity. This results implies a Landau damping phenomenon for the solution.

In this paper, we consider the time discretization of (1.1) by splitting methods based on the decomposition of the equation between the free part

(1.3) ∂ t f (t, x, v) + v∂ x f (t, x, v) = 0, f (0, x, v) = f 0 (x, v),
whose solution is given explicitely by ϕ t T (f 0 )(x, v) := f 0 (x -tv, v), and the potential part

(1.4) ∂ t f (t, x, v) = ∂ x R×T P (x -y)f (t, y, u)dudy ∂ v f (t, x, v), f (0, x, v) = f 0 (x, v),
whose solution is explicitely given by

ϕ t P (f 0 ) = f 0 (x, v + tE(f 0 , x)),
where E(f, x) = ∂ x R×T P (x -y)f (y, u)dudy is indeed kept constant during the evolution of (1.4).

The Lie splittings we consider are given by the formulas (1.5)

f n+1 = ϕ h P • ϕ h T (f n ), or f n+1 = ϕ h T • ϕ h P (f n ),
where h > 0 is the time step. The functions f n (x, v) defined above are a priori order one approximations of f (t, x, v) at time t = nh.

We also consider the Strang splitting (1.6)

f n+1 = ϕ h/2 T • ϕ h P • ϕ h/2
T (f n ) that should provide an order two approximation f n (x, v) of f (t, x, v) at time t = nh (the same being expected for the symmetric splitting where the roles of T and P are swapped).

We can then define the sequence of function r n (x, v) by the formula (1.7)

f n (x, v) = η(v) + εr n (x, v), and the functions (1.8)

g n (x, v) = r n (x + nhv, v)
which have to be thought as approximations of g(t, x, v) at time t = nh.

The main result of our paper is that if ε and h are small enough, g n (x, v) converges towards a limit function g ∞ h (x, v) when the n goes to ∞. Moreover, this solution is close to the exact limit function g ∞ (x, v) with an error estimate that scales in h for the Lie splitting, and in h 2 for the Strang splitting. Note that our results also imply convergence results in time which are uniform for positive times for g n (x, v) and give explicit convergence bounds for f n (x, v) in H s (Sobolev space, see (2.1) below) that depend on the final time T in a polynomial way.

The main idea of our proof can be compared with the classical backward error analysis methods widely used in Geometric Numerical Integration, see for instance [START_REF] Hairer | Geometric Numerical Integration[END_REF], [START_REF] Leimkuhler | Simulating Hamiltonian dynamics[END_REF]: we express the numerical solution as the exact solution of a continuous Vlasov type equation with time dependent kernel (with a poor regularity in time). Usually for Hamiltonian systems, the analysis has to be refined to make this equation independent of the time, implying the existence of a modified energy that is preserved by the numerical scheme. This "time averaging" introduces in general a remaining error term which is exponentially small (with respect to the time step) for finite dimensional systems (see [START_REF] Benettin | On the Hamiltonian interpolation of near to the identity symplectic mappings with application to symplectic integration algorithms[END_REF], [START_REF] Hairer | Geometric Numerical Integration[END_REF], [START_REF] Leimkuhler | Simulating Hamiltonian dynamics[END_REF], [START_REF] Reich | Backward error analysis for numerical integrators[END_REF]) or requires the use of a CFL (Courant-Friedrichs-Lewy) condition for semilinear Hamiltonian equations to be controlled (see [START_REF] Faou | Hamiltonian interpolation of splitting approximations for nonlinear PDE's[END_REF], [START_REF] Faou | Geometric numerical integration and Schrödinger equations[END_REF]).

Here the situation is completely different. The long time behavior of the solution is essentially controlled by a time convergent integral, which is a consequence of the dispersive effect of the free flow and ensures the existence of the continuous limit function g ∞ (x, v) (The Landau damping effect, see [START_REF] Mouhot | On Landau damping[END_REF], [START_REF] Bedrossian | Landau damping: paraproducts and Gevrey regularity[END_REF], [START_REF] Faou | Landau Damping In Sobolev Spaces For The Vlasov-HMF Model[END_REF]). As we will observe in the next section, the effect of the splitting approximation is essentially to discretize this convergent integral. As the integrand converges algebraically when the time goes to infinity, the numerical solution also yields a convergent time integral, even if the time appears in a discontinuous way in the evolution equation.

The proof of the uniform convergence estimates is based on a similar argument, but requires slightly more regularity for the functions than for the continuous case.

LANDAU DAMPING FOR THE VLASOV-HMF MODEL, MAIN RESULT

Before stating our main result, we first recall the scattering result derived in [START_REF] Faou | Landau Damping In Sobolev Spaces For The Vlasov-HMF Model[END_REF] (see also [START_REF] Mouhot | On Landau damping[END_REF], [START_REF] Bedrossian | Landau damping: paraproducts and Gevrey regularity[END_REF] for similar result with analytic or Gevrey regularity that are valid for much more singular interaction potentials).

We work in the following weighted Sobolev spaces: for a given ν > 1/2, we set

(2.1) f 2 H s ν = |p|+|q| s T×R (1 + |v| 2 ) ν |∂ p x ∂ q v f | 2 dxdv,
and we shall denote by H s ν the corresponding function space. We shall denote by • or F the Fourier transform on T × R given by (2.2)

fk (ξ) = 1 2π T×R f (x, v)e -ikx-iξv dxdv.
We shall need a stability property of the reference state η in order to control the linear part of the Vlasov equation (3.1). Let us denote by η = η(v) the spatially homogeneous stationary state and let us define the functions

(2.3) K(n, t) = -np n nt η0 (nt), K 1 (n, t) = -np n nt η0 (nt)1 t≥0 , t ∈ R, n ∈ Z,
where (p k ) k∈Z are the Fourier coefficients of the kernel P (x). We shall denote by K1 (n, τ ) = R e -iτ t K 1 (n, t) dt the Fourier transform of K 1 (n, •). We shall assume that η satisfies the following condition

(H) η(v) ∈ H 5 3 and ∃ κ > 0, inf Im τ ≤0 |1 -K1 (n, τ )| ≥ κ, n = ±1.
Note that thanks to the localization property of η in the first part of the assumption, the Fourier transform of K can be indeed continued in the half plane Im τ ≤ 0. Here, the assumption is particularly simple due to the fact that for our kernel, there are only two non-zero Fourier modes.

This assumption is very similar to the one used in [START_REF] Mouhot | On Landau damping[END_REF], [START_REF] Bedrossian | Landau damping: paraproducts and Gevrey regularity[END_REF] and can be related to the standard statement of the Penrose criterion. In particular it is verified for the states η(v) = ρ(|v|) with ρ non-increasing which are also known to be Lyapounov stable for the nonlinear equation (see [START_REF] Marchioro | A note on the nonlinear stability of a spatially symmetric Vlasov-Poisson flow[END_REF]). We also use the notation x = (1 + |x| 2 ) 1/2 for x ∈ R. In [START_REF] Faou | Landau Damping In Sobolev Spaces For The Vlasov-HMF Model[END_REF], the following result is proved:

Theorem 2.1. Let us fix s ≥ 7, ν > 1/2 and assume that η ∈ H s+4 ν satisfies the assumption (H). Assume that g(0, x, v) is in H s ν . Then there exists ε 0 > 0 and a constant C > 0 such that for every

ε ∈ (0, ε 0 ] there exists g ∞ (x, v) ∈ H s-4 ν
such that for all r ≤ s -4 and r ≥ 1,

(2.4) ∀ t ≥ 0, g(t, x, v) -g ∞ (x, v) H r ν ≤ C t s-r-3 .
In this paper, we prove the following semi-discrete version of the previous result: Theorem 2.2. Let us fix s ≥ 7, ν > 1/2 and assume that η ∈ H s+4 ν satisfies the assumption (H). Assume that g(0, x, v) is in H s ν . For a time step h, let g n (x, v), n ≥ 0, be the sequence of functions defined by the formula (1.8) from iterations of the splitting methods (1.5) (Lie), or (1.6) (Strang), with g 0 (x, v) = g(0, x, v). Then there exists ε 0 > 0, h 0 > 0 and a constant C > 0 such that for every ε ∈ (0, ε 0 ] and every h ∈ (0, h 0 ], there exists g

∞ h (x, v) ∈ H s-4 ν
such that for all r ≤ s -4 and r ≥ 1,

(2.5) ∀ n ≥ 0, g n (x, v) -g ∞ h (x, v) H r ν ≤ C nh s-r-3 .
If moreover ν > 3/2 and s ≥ 8, we have for the Lie splitting methods (1.5) the estimate

(2.6) g n (x, v) -g(nh, x, v) H s-6 ν-1 ≤ Ch ∀n ∈ AE,
where g(t, x, v) is the solution (1.2) associated with the continuous equation with the same initial value.

In the case of the Strang splitting method (1.6), we have if ν > 5/2 and s ≥ 9 the second order estimate

(2.7) g n (x, v) -g(nh, x, v) H s-7 ν-2 ≤ Ch 2 ∀n ∈ AE.
Let us make the following comments:

a) The estimates (2.6) and (2.7) exhibit a convergence rates in time of order 1 and 2 respectively for the numerical solutions. These estimates hold uniformly in time. Note however that these results do not imply convergence results uniform in time for the functions f n (x, v) to f (nh, x, v) given by the splitting methods (1.5) and (1.6). It is easy to check, using the formula f n (x, v) = η(v) + εg n (x -nhv, v) that we have an estimate of the form

∀ n ≥ 0, f n (x, v) -f (nh, x, v) H s-6 ν-1 ≤ Cε nh s-6 h.
for the Lie splitting methods (1.5). In the case of the Strang splitting (1.6), we have from the same arguments:

∀ n ≥ 0, f n (x, v) -f (nh, x, v) H s-7 ν-2 ≤ Cε nh s-7 h 2 .
Hence we obtain convergence results which are global in time only if we measure the error in L 2 . If we measure the error in H σ , σ > 0, then for a fixed time horizon nh ≤ T , the error grows like T σ . This is however better than the rough e T , estimate that is usually obtained through Gronwall type arguments (Note that convergence results can be found in [START_REF] Einkemmer | Convergence analysis of Strang splitting for Vlasov-type equations[END_REF] for the case of compactly supported data, and in [START_REF] Casas | High-order Hamiltonian splitting for Vlasov-Poisson equations[END_REF] for the Vlasov-Poisson case). Let us also mention as an easy consequence of (2.4), (2.5), (2.6) and (2.7), that the following estimates hold for the limit state of the equation: For the Lie splitting

g ∞ (x, v) -g ∞ h (x, v) H s-6 ν-1 ≤ Ch,
and for the Strang splitting

g ∞ (x, v) -g ∞ h (x, v) H s-7 ν-2 ≤ Ch 2 .
b) The long time behavior of the exact solution (1.2) is essentially controlled by a time convergent integral (see [START_REF] Faou | Landau Damping In Sobolev Spaces For The Vlasov-HMF Model[END_REF]). We shall see (Proposition 3.1 below) that the splitting method provides a discretization of this integral, but essentially without changing the decay in time of the integrand. Thus the numerical solution also yields a convergent time integral, even if the time appears in a discontinous way, giving us the long time behavior. Moreover, this discretization is performed by rectangle methods in the case of Lie splittings, and by the midpoint rule in the case of Strang splitting. Estimates (2.6) and (2.7) reflect the respective accuracies of these two methods. The second order estimate requires a more refined analysis than the first order, for it is obtained by tracking the cancellations provided by the midpoint rules. We mention that this result remains true for Strang splitting of the form (1.6) where the role of T and P are exchanged but the complete proof is given for (1.6) only (the time integration rule being the trapezoidal rule and the arguments identical for both cases). Finally, let us mention that the proofs of the convergence results (for Lie or Strang) widely use the long time behavior of both the exact and discrete solutions, in particular uniform bounds on their regularity. This can be understood as stability results for the numerical schemes. The convergence results are essentially the combination of these stability results, and the accuracy of the discretization of the integral. c) Our results hold only for time discretization of the equation. Fully discrete scheme including for example an interpolation procedure at each step (semi-Lagrangian methods) traditionally exhibit recurrence phenomena due the discretization in the v variable. Indeed, the Landau damping effect reflects essentially the fact that the solution of the free Vlasov equation is a superposition of travelling wave in the Fourier variable ξ. At the discrete level, the ξ variable is only discretized by a finite number of points which causes numerical interactions of these travelling waves preventing the mixing effect to occur for very long time. Typically, the previous result is hence valid a priori for a time of order O(1/δv) only, if δv is the size of the mesh variable in v. Solutions exist to remedy these difficulties, for examples by putting absorbers in the Fourier spaces, see for instance [START_REF] Einkemmer | A strategy to suppress recurrence in grid-based Vlasov solvers[END_REF]. The analysis of these space discretization effects will be the subject of further studies.

Let us finally explain how the previous scattering results imply Landau damping effects for the solution f (t, x, v). Let us recall the following elementary Lemma: Lemma 2.3. For every α, β, γ, s ∈ N with α + β = s, and γ < ν -1 2 . we have the following inequality:

(2.8) ∀ k ∈ Z, ∀ ξ ∈ R, |∂ γ ξ fk (ξ)| 2 s/2 C(ν) k -α ξ -β f H s ν ,
where C(ν) depends only on ν > 1/2.

Proof. We have by using the Cauchy-Schwarz inequality that

k α ξ β ∂ γ ξ fk (ξ) = 1 2π T×R ∂ α x ∂ β v (v γ f (x, v))e -ikx e -ivξ dxdv C f H s ν R (1 + |v| 2 ) γ-ν dv 1/2
.

The previous inequality with α = β = 0 yields the result when k = 0 or |ξ| 1 and we conclude by using x 2 α/2 |x| 2 for |x| > 1 and the fact that ν -γ > 1/2.

As a consequence of Theorem 2.2, we have the nonlinear Landau damping effect for the semidiscrete solution: The functions g n (x, v) being bounded in

H s-4 ν , f n (x, v) = η(v) + εr n (x, v) = η(v) + εg n (x -nhv, v) satisfy ∀ k ∈ Z * , ∀ξ ∈ R, ∀ α + β = s -4, | f n k (ξ)| = ε|ĝ n k (ξ + knh)| ≤ Cε ξ + knh α k β ,
the last estimate being a consequence of the embedding Lemma 2.3. This yields that for every k = 0, f n k (ξ) tends to zero when nh → ∞ with a polynomial rate, but with a speed depending on k.

Moreover, by setting

η ∞ h (v) := η(v) + ε 2π T g ∞ h (x, v
)dx, we have by the previous Lemma 2.3 that for r ≤ s -4,

∀ ξ ∈ R, | f n 0 (ξ) -η∞ h (ξ)| ≤ C ξ r nh s-r-3 .
In other words,

f n (x, v) converges weakly towards η ∞ h (v). Moreover, this weak limit η ∞ h (v) is O(h) for Lie splittings (or O(h 2 ) for Strang splitting) close to the exact limit η ∞ (v) := η(v) + ε 2π T g ∞ (x, v)dx,
which exists by Theorem 2.1.

BACKWARD ERROR ANALYSIS

The unknown g(t, x, v) defined in (1.2) is solution of the equation

(3.1) ∂ t g = {φ(t, g(t)), η} + ε{φ(t, g(t)), g}.
where

(3.2) φ(t, g)(x, v) = R×T (cos(x -y + t(v -u)))g(y, u)dudy,
and {f, g} = ∂ x f ∂ v g -∂ v f ∂ x g is the usual microcanonical Poisson bracket.
In Fourier space, we have the expression:

φ(t, g) = 1 2 k∈{±1}
e ikx e iktv ĝk (kt).

In the evolution of the solution g(t, x, v) of (3.1), an important role is played by the quantity

(3.3) ζ k (t) = ĝk (t, kt), k ∈ {±1}, such that φ(t, g(t)) = 1 2 k∈{±1} e ikx e iktv ζ k (t). Note that for k = 0, ζ k (t) is the Fourier coefficient in x of the density ρ(t, x) = R f (t, x, v) dv.
The following result shows that semi-discrete solution g n (x, v) also satisfies an equation of the form (3.1), but with a discontinuous dependence with respect to the time: Proposition 3.1. For h > 0, the solution g n (x, v) given by the splitting method (1.5) (Lie) or (1.6) (Strang), and formula (1.8) coincides at times t = nh with the solution g(t, x, v) of the equation

(3.4) ∂ t g = {Φ h (t, g(t)), η} + ε{Φ h (t, g(t)), g(t)}, where Φ h (t, g(t)) = φ(s h (t), g(t)
) with the definition of φ given in (3.2).

In the case of Lie splittings (1.5) , s h (t) are given by the formulas

(3.5) s h (t) = t h h + h, and s h (t) = t h h,
respectively. In the case of the Strang splitting method (1.6), s h (t) is given by the formula

(3.6) s h (t) = t h h + h 2 .
Proof. We prove the result in the case of Strang splitting (1.6)-(1.8), the proof being analogous for Lie splittings. By definition, the function f n (x, v) satisfies the recurrence relation (1.6). Hence, we have (using the linearity of ϕ t T and the fact that

ϕ t T (η) = η for all t ∈ R) η + εg n = ϕ -nh T (f n ) = ϕ -nh T • ϕ h/2 T • ϕ h P • ϕ h/2 T (f n-1 ) = ϕ -nh+h/2 T • ϕ h P • ϕ h/2+(n-1)h T ϕ -(n-1)h T (f n-1 ) = ϕ -nh+h/2 T • ϕ h P • ϕ nh-h/2 T (η + εg n-1
).

Now we verify that for

t ∈ [0, h], the application t → ϕ -nh+h/2 T • ϕ t P • ϕ nh-h/2 T (η + εg n-1
) is the solution of the equation ∂ t g = {φ(s h (t), g), g}.

with inital data g(0) = η + εg n-1 . Using the fact that η is a stationary state of the equation, we easily get the result.

For notational convenience, we will often write in the following s(t) instead of s h (t). As in (3.3), we define

(3.7) Z k (t) = ĝk (t, ks(t)), k ∈ {±1}, such that Φ h (t, g(t)) = φ(s(t), g(t)) = 1 2 k∈{±1} e ikx e iks(t)v Z k (t).
Of course, we expect the Z k (t) to be approximations of the terms ζ k (t) defined in (3.3).

Lemma 3.2. Let h 0 > 0 be given. There exist two constants c and C > 0 such that for all h ∈ (0, h 0 ] and all t > 0,

c t ≤ s h (t) ≤ C t
and for all t and σ, c t

± σ ≤ s h (t) ± s h (σ) ≤ C t ± σ .
Proof. For t ∈ R, we can write t = nh + µ with µ ∈ [0, h). In the case of Strang splitting (3.6), we thus have

s(t) = nh + h/2 = t + h/2 -µ. Hence we have t -s(t) ∈ [-h/2, h/2)
which clearly implies the first inequality. The second is proved using the fact that with similar calculations t ± σ = s(t) ± s(σ) + O(h). The proof is analogous for Lie splittings (3.5).

As in [START_REF] Faou | Landau Damping In Sobolev Spaces For The Vlasov-HMF Model[END_REF] we introduce the weighted norms:

(3.8) N T,s,ν (g) = sup t∈[0,T ] g(t) H s ν t 3 , M T,γ (Z) = sup t∈[0,T ] sup k∈{±1} t γ |Z k (t)| and (3.9) Q T,s,ν (g) = N T,s,ν (g) + M T,s-1 (Z) + sup [0,T ] g(t) H s-4 ν .
We shall prove the following result:

Theorem 3.3. Let us fix s ≥ 7, ν > 1/2 and R 0 > 0 such that Q 0,s,ν (g) ≤ R 0 , and assume that η ∈ H s+4 ν
satisfies the assumption (H). Then there exist R > 0, h 0 > 0 and ε 0 > 0 such that for every ε ∈ (0, ε 0 ], h ∈ (0, h 0 ] and for every T ≥ 0, the solution of (3.4) satisfies the estimate

Q T,s,ν (g) ≤ R.
This result is a semi-discrete version of the main Theorem in [START_REF] Faou | Landau Damping In Sobolev Spaces For The Vlasov-HMF Model[END_REF] where the same norms are used to control the solution g(t) of the equation (3.1). Let us also mention that it holds for any of the three formulas (3.6)-(3.5) defining s(t), the only property being used is the fact that s(t) satisfies Lemma 3.2.

ESTIMATES

We fix now s ≥ 7 and R 0 as in the previous Theorem. In the following a priori estimates, C stands for a number which may change from line to line and which is independent of R 0 , R, h, ε and T . 4.1. Estimate of M T,s-1 (Z). Towards the proof of Theorem 3.3, we shall first estimate Z k (t), k = ±1. Proposition 4.1. Assuming that η ∈ H s+2 ν verifies the assumption (H), then there exist C > 0 and h 0 > 0 such that for every T > 0 and h ∈ (0, h 0 ], every solution of (3.4)

such that Q T,s,ν (g) ≤ R enjoys the estimate (4.1) M T,s-1 (Z) ≤ C R 0 + εR 2 .
Proof. The main ingredient of the proof of the previous result is to write the equation (3.4) in Fourier:

(4.2) ĝn (t, ξ) = ĝn (0, ξ) + t 0 p n Z n (σ)η 0 (ξ -ns(σ))(n 2 s(σ) -nξ)dσ + ε k∈Z p k t 0 Z k (σ)ĝ n-k (σ, ξ -ks(σ))(nks(σ) -kξ)dσ,
for all (n, ξ) ∈ Z × R, with p k = 1 2 for k ∈ {±1} and p k = 0 for k = ±1, and where the Z k (t) are defined by (3.7). Setting ξ = ns(t) in (4.2), the equation satisfied by (Z n (t)) n=±1 can be written under the almost closed form

(4.3) Z n (t) = ĝn (0, ns(t)) - t 0 p n Z n (σ)η 0 (n(s(t) -s(σ)))n 2 (s(t) -s(σ))dσ -ε k∈{±1} p k t 0 Z k (σ)ĝ n-k (σ, ns(t) -ks(σ))kn(s(t) -s(σ))dσ.
Remark 4.2. Note that, for every m ∈ AE, we have for t ∈ [mh, (m + 1)h] the formula

Z n (t) -ĝn (mh, ns(t)) = - t mh p n Z n (σ)η 0 (n(s(t) -s(σ)))n 2 (s(t) -s(σ))dσ -ε k∈{±1} p k t mh Z k (σ)ĝ n-k (σ, ns(t) -ks(σ))kn(s(t) -s(σ))dσ.
As s(t)s(σ) = 0 for almost every σ ∈ [mh, t], we notice that the function t → ĝn (t, ns(t)) = Z n (t) is constant on the small intervalls [mh, (m + 1)h]. This is due to the fact that the electric field is constant during the evolution of (1.4).

To study the equation (4.3), we shall first consider the corresponding linear equation, that is to say that we shall first see

(4.4) F n (t) := ĝn (0, ns(t)) -ε k∈{±1} p k t 0 Z k (σ)ĝ n-k (σ, ns(t) -ks(σ))kn(s(t) -s(σ))dσ
as a given source term and we shall study the linear integral equation

(4.5) Z n (t) = t 0 K(n, s(t) -s(σ))Z n (σ) dσ + F n (t) n = ±1,
where the kernel K(n, t) has been introduced in (2.3). We rewrite this equation as

(4.6) Z n (t) = t 0 K(n, t -σ)Z n (σ) dσ + F n (t) + G n (t) n = ±1, where (4.7) G n (t) = t 0 K(n, s(t) -s(σ)) -K(n, t -σ) Z n (σ) dσ
is an error term due to the time discretization. To study the linear equation (4.5)-(4.6), we use the result given by the following Lemma. The proof of this Lemma can be found in [START_REF] Faou | Landau Damping In Sobolev Spaces For The Vlasov-HMF Model[END_REF]. For completeness, we recall it in Appendix of the present paper.

Lemma 4.3. Let γ ≥ 0, and assume that η ∈ H γ+3 ν satisfies (H). Then, there exists C > 0 such for every T ≥ 0, we have

M T,γ (Z) ≤ CM T,γ (F + G).
From this Lemma and (2.8), we first get that

(4.8) M T,s-1 (Z) ≤ C g(0) H s ν + M T,s-1 (G) + εM T,s-1 (F 1 ) + εM T,s-1 (F 2 ) where F 1 n (t) = n 2 p -n t 0 Z -n (σ)ĝ 2n (σ, n(s(t) + s(σ)))(s(t) -s(σ)) dσ, n = ±1, F 2 n (t) = -n 2 p n t 0 Z n (σ)ĝ 0 (σ, n(s(t) -s(σ)))(s(t) -s(σ)) dσ, n = ±1.
Let us estimate F 1 n . By using again (2.8) and the definition (3.8) of N σ,s,ν , we get using Proposition 3.1 that

|F 1 n (t)| ≤ C t 0 (s(t) -s(σ)) σ 3 M σ,s-1 (Z)N σ,s,ν (g) σ s-1 s(t) + s(σ) s dσ ≤ C R 2 t s-1 +∞ 0 1 σ s-4 dσ ≤ C R 2 t s-1
provided s ≥ 6. This yields that for all T ≥ 0

M T,s-1 (F 1 ) ≤ CR 2 .
To estimate F 2 n , we split the integral into two parts: we write

F 2 n (t) = I 1 n (t) + I 2 n (t) with I 1 n (t) = -n 2 p n t 2 0 Z n (σ)ĝ 0 (σ, n(s(t) -s(σ)))(s(t) -s(σ)) dσ, n = ±1, I 2 n (t) = -n 2 p n t t 2 Z n (σ)ĝ 0 (σ, n(s(t) -s(σ)))(s(t) -s(σ)) dσ, n = ±1.
For I 1 n , we proceed as previously,

|I 1 n (t)| ≤ CR 2 t 2 0 σ 3 (s(t) -s(σ)) σ s-1 s(t) -s(σ) s dσ ≤ CR 2 t s-1 +∞ 0 1 σ s-4 dσ
and hence since s ≥ 6, we have M T,s-1 (I 1 ) ≤ CR 2 . To estimate I 2 n , we shall rather use the last factor in the definition of Q T,s,ν in (3.9). By using again (2.8), we write

|I 2 n (t)| ≤ C t t 2 M σ,s-1 (Z) σ s-1 g(σ) H s-4 ν s(t) -s(σ) s-5 dσ ≤ CR 2 t s-1 +∞ 0 1 σ s-5 dσ ≤ CR 2 t s-1
and hence since s ≥ 7, we find again

M T,s-1 (I 2 ) ≤ CR 2 .
Finally, we have

t s-1 |G n (t)| ≤ M T,s-1 (Z) t 0 t s-1 σ s-1 s(t)-s(σ) t-σ |∂ t K(n, θ)|dθ dσ ≤ CM T,s-1 (Z) t 0 s(t)-s(σ) t-σ t s-1 θ s+2 σ s-1 dθ dσ,
where we used the fact that η ∈ H s+2 ν and Lemma 2.3. Now, since

s(t)-s(σ) t-σ dθ θ s+2 ≤ Ch t -σ s+2 , we get M T,s-1 (G) ≤ ChM T,s-1 (Z) t 0 t s-1 t -σ s+2 σ s-1 dσ. As the integral t 0 t s-1 t -σ s+2 σ s-1 dσ ≤ C t/2 0 1 t 3 σ s-1 dσ + C t t/2 1 t -σ s+2 dσ.
is uniformly bounded in time, we conclude that (4.9)

M T,s-1 (G) ≤ ChM T,s-1 (Z).
By combining the last estimates and (4.8), we thus obtain (4.1).

M T,s-1 (Z) ≤ C R 0 + hM T,s-1 (Z) + εR 2 .
By taking h ≤ h 0 small enough, this ends the proof of Proposition 4.1.

4.2.

Estimate of N T,s,ν (g).

Proposition 4.4. Assuming that η ∈ H s+2 ν verifies the assumption (H), then there exists C > 0 and h 0 > 0 such that for every T > 0 and h ∈ (0, h 0 ], every solution of (3.4)

such that Q T,s,ν (g) ≤ R enjoys the estimate N T,s,ν (g) ≤ C(R 0 + εR 2 )(1 + εR)e CεR .
Proof.

2 L 2 = 2ε Dg(t), D(L s(t) [g(t)]g(t)) L 2 + 2 Dg(t), D(L s(t) [g(t)](η)) L 2 = 2ε Dg(t), L s(t) [g(t)]Dg(t) L 2 + 2ε Dg(t), [D, L s(t) [g(t)]]g(t) L 2 +2 Dg(t), D(L s(t) [g(t)](η)) L 2 ,
where [D, L s(t) ] denotes the commutator between the two operators D and L s(t) . The first term in the previous equality vanishes since L s(t) [g] is the transport operator associated with a divergence free Hamiltonian vector field. Consequently, we get that (4.10)

d dt Dg(t) 2 L 2 2ε Dg(t) L 2 [D, L s(t) [g(t)]]g(t) L 2 + 2 Dg(t) L 2 D(L s(t) [g(σ)](η)) L 2 .
To get the estimates of Proposition 4.4, we shall use the previous estimates with the operator D = D m,p,q defined as the Fourier multiplier by k p ξ q ∂ m ξ for (m, p, q) ∈ N 3d such that p + q s, m ν and the definition (2.1) of the H s ν norm. To evaluate the right hand-side of (4.10), we shall use the following Lemma, whose proof is given in Appendix (see also [START_REF] Faou | Landau Damping In Sobolev Spaces For The Vlasov-HMF Model[END_REF]). Lemma 4.5. For p + q γ and r ν, and functions h(x, v) and g(x, v), we have the estimates

D r,p,q , L σ [g] h L 2 ≤ C m σ,γ+1 (ζ) h H 1 ν + m σ,2 (ζ) h H γ ν , (4.11) D r,p,q L σ [g] h L 2 ≤ C m σ,γ+1 (ζ) h H 1 ν + m σ,2 (ζ) h H γ+1 ν , (4.12)
for all σ, where the sequence ζ k is defined by ζ k = ĝk (kσ), k ∈ {±1}, and where

m σ,γ (ζ) = σ γ sup k∈{±1} |ζ k | ,
with a constant C depending only on γ, and in particular, not depending on σ.

Let us finish the proof of Proposition 4.4. By using the previous lemma with γ = s, σ = s(t), g = g(t) and h = g(t) or h = η, we obtain from (4.10) that

d dt g(t) 2 H s ν ≤ C t 2 m t,s-1 (Z(t)) η H 1 ν + ε g(t) H 1 ν g(t) H s ν + C t s-3 m t,s-1 (Z(t)) η H s+1 ν g(t) H s ν + Cε t s-3 m t,s-1 (Z(t)) g(t) 2 H s ν .
This yields, using the fact that M t,γ (Z) = sup σ∈[0,t] m σ,γ (Z(σ)),

g(t) H s ν ≤ g(0) H s ν + C t 3 M t,s-1 (Z) η H s+1 ν + εR + CεR t 0 1 σ s-3 g(σ) H s ν dσ for t ∈ [0, T ].
From the Gronwall inequality, we thus obtain

g(t) H s ν ≤ g(0) H s ν + C t 3 M t,s-1 (Z) η H s+1 ν + εR e CεR +∞ 0 dσ σ s-3 .
By using Proposition 4.1, this yields

N T,s,ν (g) ≤ R 0 + C(R 0 + εR 2 )(1 + εR) e CεR .
This ends the proof of Proposition 4.4. verifies the assumption (H), then there exists C > 0 and h 0 > 0 such that for every T > 0 and h ∈ (0, h 0 ], every solution of (3.4)

such that Q T,s,ν (g) ≤ R enjoys the estimate g(t) H s-4 ν ≤ C R 0 + εR 2 )e CεR , ∀t ∈ [0, T ].
Proof. We use again (4.10) with D = D m,p,q but now with p + q ≤ s -4. By using Lemma 4.5, we find

(4.13) d dt g(t) 2 H s-4 ν ≤ Cm t,s-3 (Z(t)) η H s-3 ν g(t) H s-4 ν + ε g(t) 2 H s-4 ν .
This yields

g(t) H s-4 ν ≤ g(0) H s-4 ν + C η H s-3 ν M t,s-1 (Z) t 0 1 σ 2 dσ + CεM t,s-1 (Z) t 0 1 σ 2 g(σ) H s-4 ν dσ.
By using Proposition 4.1, we thus get

g(t) H s-4 ν ≤ C R 0 + εR 2 ) + CεR t 0 1 σ 2 g(σ) H s-4 ν dσ.
From the Gronwall inequality, we finally find

g(t) H s-4 ν ≤ C R 0 + εR 2 )e CεR .
This ends the proof of Proposition 4.6.

Q T,s,ν (g) ≤ C(R 0 + εR 2 )(1 + εR)e CεR
assuming that Q T,s,ν (g) ≤ R. Consequently, let us choose R such that R > CR 0 , then for ε sufficiently small we have R > C(R 0 + εR 2 )(1 + εR)e CεR and hence by usual continuation argument, we obtain that the estimate Q T,s,ν (g) ≤ R is valid for all times.

To prove (2.5), let us define g ∞ h (x, v) by

g ∞ h (x, v) = g(0, x, v) + +∞ 0 L s(σ) [g(σ)](η + ǫg(σ))dσ,
To justify the convergence of the integral, we note that thanks to (4.12), we have for all σ

L s(σ) [g(σ)](η + ǫg(σ)) H s-4 ν ≤ C m σ,s-3 (Z(σ)) η + ǫg(σ) H 1 ν +m σ,2 (Z(σ)) η + ǫg(σ) H s-3 ν ,
where Z k (t) = ĝk (t, ks(t)), k = ±1, and

m σ,γ (Z(σ)) = σ γ sup k=±1 |Z k (σ)| .
By interpolation, we have

η + ǫg(σ) H s-3 ν ≤ C η + ǫg(σ) 3 4 H s-4 ν η + ǫg(σ) 1 4
H s ν , and thus, using the bound provided by Theorem 3.3, we have

η + ǫg(σ) H s-3 ν ≤ C(R) σ 3 4 .
Using again Theorem 3.3, we thus find that

(4.14) L s(σ) [g(σ)](η + ǫg(σ)) H s-4 ν ≤ C(R) 1 σ 2 + 1 σ s-3-3 4
, and g ∞ h (x, v) is then well defined, and belongs to H s-4 ν . Since we have for all t

g(t) -g ∞ h (x, v) = +∞ t L s(σ) [g(σ)](η + ǫg(σ))dσ,
we find by using again (4.14) that

g(t) -g ∞ h H s-4 ≤ C(R) +∞ t 1 σ 2 + 1 σ s-3-3 4 dσ ≤ C(R) t .
In a similar way, by using again (4.12), we have for r ≤ s -4 and r ≥ 1,

g(t) -g ∞ h H r ν ≤ C(R) +∞ t 1 σ s-r-2 + 1 σ s-3 dσ ≤ C(R) 1 t s-r-3 + 1 t s-4 ≤ C(R) t s-r-3 ,
which gives the result using the fact that g(nh) = g n the solution given by the numerical scheme.

PROOF OF THE CONVERGENCE ESTIMATE (2.6)

We shall now prove the convergence estimate (2.6). Note that in view of the previous result and of the analysis in [START_REF] Faou | Landau Damping In Sobolev Spaces For The Vlasov-HMF Model[END_REF], the functions g(t, x, v) and g(t, x, v) satisfy the same estimates. In particular, we can assume that Q T,s,ν (g), Q T,s,ν-1 (g), Q T,s,ν (g) and Q T,s,ν-1 (g) are both uniformly bounded by the same constant R, provided that ν > 3/2.

By using the equations (3.1) and (3.4), we get that δ = gg solves the equation

(5.1) ∂ t δ(t) = {φ(t, δ(t)), η} + ε{φ(t, δ(t)), g(t)} + ε{φ(t, g(t)), δ(t)} -ε{φ(t, δ(t)), δ(t)} + R with (5.2) R(t, x, v) = {φ(t, g(t)) -φ(s(t), g(t)), η} + ε{φ(t, g(t)) -φ(s(t), g(t)), g(t)}
and with zero initial data. It will be useful to use the expression of R in Fourier which reads

(5.3) Rn (t, ξ) = np n ĝn (t, nt)η 0 (ξ -nt)(nt -ξ) -ĝn (t, ns(t))η 0 (ξ -ns(t))(ns(t) -ξ) + ε k=±1 kp k ĝk (t, kt)ĝ n-k (t, ξ -kt)(nt -ξ) -ĝk (t, ks(t))ĝ n-k (t, ξ -ks(t))(ns(t) -ξ) .
Let T be a positive real number. By using the weighted norms defined in (3.8), we now consider for s ≥ 8 and ν > 3/2, the quantity

Q T,s-2,ν-1 (δ) = M T,s-3 (d) + N T,s-2,ν-1 (δ) + sup t∈[0,T ] δ(t) H s-6 ν-1
with N T,s,γ and M T,γ defined in (3.8), and where we set

d k (t) = ĝk (t, kt) -ĝk (t, kt).
We shall prove the following result: Proposition 5.1. For s ≥ 8, ν > 3/2, assume that η ∈ H s+4 ν satisfies the assumption (H). Then there exists R 1 > 0, h 0 > 0 and ǫ 0 > 0 such that for every h ∈ (0, h 0 ], every ε ∈ (0, ε 0 ] and every T ≥ 0, the solution of (5.1) satisfies the estimate

Q T,s-2,ν-1 (δ) ≤ R 1 h.
Proposition 5.1 clearly implies the convergence estimate (2.6). It actually proves that the interpolation g(t) of the sequence of functions g n (x, v) given by the splitting methods (Strang or Lie) is always at least an approximation of order one of the exact solution g(t) at all times. The proof will use the same steps as in the proof of Theorem 3.3. Remark 5.2. We mentioned that Z k (t) is expected to be an approximation of ζ k (t). By Taylor expanding in ξ, it is indeed clear that Proposition 5.1, together with the uniform bound on Q T,s,ν (g) and Q T,s,ν (g), implies that sup

t∈Ê k=±1 |ζ k (t) -Z k (t)| = O(h).
In fact, we could have proved without any major difference Proposition 5.1 using the quantity

dk (t) = ζ k (t) -Z k (t) instead of d k (t)
. The use of d k (t) will however be crucial to prove the second order estimate, since it is clear that the quantity dk (0) does not scale in h 2 .

Let us now begin the proof of Proposition 5.1.

Estimate of M T,s-3 (d).

We shall first estimate d k (t), k = ±1.

Proposition 5.3. Assuming that η ∈ H s+4 ν satisfies the assumption (H), there exists C > 0, ε 0 > 0 and h 0 > 0 such that for every h ∈ (0, h 0 ], every ε ∈ (0, ε 0 ] and every T > 0, every solution of (5.1)

such that Q T,s-2,ν-1 (δ) ≤ R 1 h enjoys the estimate M T,s-3 (d) ≤ C(R)h(1 + (ε + εh)R 2 1 )
, where C(R) is a number that depends only on R (one can take

C(R) = (1 + R + R 2 )C).
Proof. From (5.1), we obtain by taking the Fourier transform, integrating in time and setting ξ = nt that (5.4)

d n (t) = t 0 K(n, t -σ)d n (σ) dσ + G n (t) + H n (t)
where

H n (t) = t 0 Rn (σ, nt) dσ and G n (t) = ε k=±1 t 0 kp k d k (σ)ĝ n-k (σ, nt -kσ) + ζ k (σ) δn-k (σ, nt -kσ) -d k (σ) δn-k (σ, nt -kσ) n(σ -t) dσ.
The kernel K(k, t) is still defined by (2.3). By using Lemma 4.3, we find the estimate

(5.5) M T,s-3 (d) ≤ C (M T,s-3 (G) + M T,s-3 (H)) .
To estimate M T,s-3 (G), we proceed as in the proof of Theorem 3.3. We first split

G = G 1 + G 2
where G 1 corresponds to the term k = -n in the sum and G 2 corresponds to k = n. For G 1 , we obtain

|G 1 n | ≤ Cε t 0 (t -σ) 1 σ s-3 σ 3 t + σ s M σ,s-3 (d)N σ,s,ν (g) + 1 σ s-1 σ 3 t + σ s-2 M σ,s-1 (ζ)N σ,s-2,ν-1 (δ) + 1 σ s-3 σ 3 t + σ s-2 M σ,s-3 (d)N σ,s-2,ν-1 (δ) dσ,
and hence that

|G 1 n | ≤ Cε(RR 1 h + h 2 R 2 1 ) 1 t s-2 since s ≥ 8. To estimate G 2 n , we can again write G 2 n = J 1 n + J 2 n with J 1 n = t 2 0 j n dσ, J 2 n = t t 2 j n dσ with j n = εnp n d n (σ)ĝ 0 (σ, nt -nσ) + ζ n (σ) δ0 (σ, nt -nσ) -d n (σ) δ0 (σ, nt -nσ) n(σ -t).
As in the proof of Proposition 4.1, we can prove by using the same estimates as above that

|J 1 n | ≤ Cε(RR 1 h + R 2 1 h 2 ) 1 t s-2 .
To estimate J 2 n , we also proceed as in the proof of Proposition 4.1 and write

|J 2 n | ≤ Cε t t 2 (t-σ) 1 σ s-3 1 t -σ s-4 M σ,s-3 (d) g H s-4 ν + 1 σ s-1 1 t -σ s-6 M σ,s-1 (ζ) δ H s-6 ν-1 + 1 σ s-3 1 t -σ s-6 M σ,s-3 (d) δ H s-6 ν-1 dσ.
This yields, since s ≥ 8,

|J 2 n | ≤ Cε(RR 1 h + R 2 1 h 2 ) 1 t s-3 . We have thus proven that M T,s-3 (G) ≤ Cε(RR 1 h + R 2 1 h 2 ). It remains to estimate M T,s-3 (H). We shall prove that M T,s-3 (H) ≤ C(R)h.
At first, we can write (5.6)

H n = H L n + εH P n with H L n = t 0 np n ĝn (σ, nσ)η 0 (nt -nσ)(nσ -nt) -ĝn (σ, ns(σ))η 0 (nt -ns(σ))(ns(σ) -nt) dσ and
H P n = t 0 k=±1 kp k ĝk (σ, kσ)ĝ n-k (σ, nt -kσ)(nσ -nt) -ĝk (σ, ks(σ))ĝ n-k (σ, nt -ks(σ))(ns(σ) -nt) dσ.
We shall focus on the estimate of H P n , the estimate of H L n is easier to obtain since η can be assumed as smooth as we need (here η ∈ H s+4 ν suffices). By Taylor expanding up to the first order, we have

H P n ≤ h k=±1 p k t 0 sup α+β≤1 sup ξ∈[ks(σ);kσ] ∂ α ξ ĝk (σ, ξ)∂ β ξ ĝn-k (σ, nt -ξ) |nσ -nt| dσ.
As previously, we distinguish the cases k = n and k = -n.

When k = -n, we have by Lemma 2.3

t 0 sup α+β≤1 sup ξ∈[-ns(σ);-nσ] ∂ α ξ ĝ-n (σ, ξ)∂ β ξ ĝ2n (σ, nt -ξ) |σ -t| dσ ≤ C t 0 σ 6 N T,s,ν (g)N T,s,ν (g) σ s t + s s-1 dσ ≤ CR 2 t s-1 , since s ≥ 8.
When k = n, we split as previously the integral into

t 0 sup α+β≤1 sup ξ∈[ns(σ);nσ] ∂ α ξ ĝn (σ, ξ)∂ β ξ ĝ0 (σ, nt -ξ) |nσ -nt| dσ = t/2 0 sup α+β≤1 sup ξ∈[ns(σ);nσ] ∂ α ξ ĝn (σ, ξ)∂ β ξ ĝ0 (σ, nt -ξ) |nσ -nt| dσ + t t/2 sup α+β≤1 sup ξ∈[ns(σ);nσ] ∂ α ξ ĝn (σ, ξ)∂ β ξ ĝ0 (σ, nt -ξ) |nσ -nt| dσ.
Using the same estimates as above, we have for the first term

∂ α ξ ĝn (σ, ξ)∂ β ξ ĝ0 (σ, nt -ξ) |nσ -nt| dσ ≤ CR 2 t s-1 .
For the second term, we rather use the quantity g H s-4 ν instead of N T,s,ν (g). We obtain the estimate

∂ α ξ ĝn (σ, ξ)∂ β ξ ĝ0 (σ, nt -ξ) |nσ -nt| dσ ≤ C t t/2 σ 3 N T,s,ν (g) g(σ) H s-4 ν σ s t -s s-5 dσ ≤ CR 2 t s-3 , since s ≥ 8. Hence M T,s-3 (H P ) ≤ ChR 2 ,
and by the same arguments

M T,s-3 (H L ) ≤ ChR 2 . Therefore M T,s-3 (H) ≤ ChR 2 ,
which concludes the proof of proposition 5.3.

Estimate of δ H s-6

ν-1 .

Proposition 5.4. Assuming that η ∈ H s+4 satisfies the assumption (H), there exists C > 0, ε 0 > 0 and h 0 > 0 such that for every h ∈ (0, h 0 ], every ε ∈ (0, ε 0 ] and every T > 0, every solution of (5.1)

such that Q T,s-2,ν-1 (δ) ≤ R 1 h enjoys the estimate sup t∈[0,T ] δ(t) H s-6 ν-1 ≤ C(R)h(1 + (ε + εh)R 2 1 )(1 + hR 1 )
where C(R) is a number that depends only on R.

Proof. By using the notation L t [g]h = {φ(t, g), h}, we can rewrite (5.1) as

(5.7) ∂ t δ(t) = L t [δ(t)](η -εδ) + εL t [δ(t)]g + εL t [g(t)]δ + R.
Let D be the linear operator defined as the Fourier multiplier by k p ξ q ∂ m ξ for (m, p, q) ∈ AE 3d such that p + q ≤ s -6, and m ≤ ν -1. From an energy estimate as in proof of Proposition 4.4, we obtain that

(5.8) Dδ(t) 2 L 2 = t 0 -ε [D, L σ [δ(σ)]]δ(σ), Dδ(σ) L 2 + ε [D, L σ [g(σ)]]δ(σ), Dδ(σ) L 2 + D (L σ [δ(σ)](η + εg(σ))) , Dδ(σ) L 2 + ε DR, Dδ(σ) L 2 dσ.
By using Lemma 4.5, we thus obtain that

sup [0,T ] δ(t) H s-6 ν-1 ≤ C T 0 ε(m σ,s-5 (d(σ)) δ(σ) H 1 ν-1 + m σ,2 (d(σ)) δ(σ) H s-6 ν-1 ) + ε(m σ,s-5 (ζ(σ)) δ(σ) H 1 ν-1 + m σ,2 (ζ(σ)) δ(σ) H s-6 ν-1 ) + m σ,s-5 (d(σ)) η + εg(σ) H 1 ν-1 + m σ,2 (d(σ)) η + εg(σ) H s-5 ν-1 + DR L 2 dσ.
Next, we can use the fact that Q t,s,ν (g) ≤ R and Proposition 5.3 to obtain that

sup t∈[0,T ] δ(t) H s-6 ν-1 ≤ C(R)h(1 + (ε + εh)R 2 1 )(1 + hR 1 ) T 0 1 σ 2 dσ + C T 0 R(σ) H s-6 ν-1 dσ ≤ C(R)h(1 + (ε + εh)R 2 1 )(1 + hR 1 ) + C T 0 R(σ) H s-6
ν-1 dσ. It remains to estimate the last integral to conclude. By using the expression (5.3), we get as in the proof of Lemma 4.5 that

R(σ) H s-6 ν-1 ≤ Ch m (1) σ,s-5 (g(σ))( η H 1 ν + g(σ) H 1 ν ) + m (1) σ,2 (g(σ))( η H s-5 ν + g(σ) H s-5 ν )
where we have set

m (1) σ,γ (h) = σ γ sup n=±1 sup ξ∈[nσ,ns(σ)] (| ĥn (ξ)| + |∂ ξ ĥn (ξ)|).
Note that thanks to lemma 2.3, we have that (5.9) m

σ,s-5 (g(σ)) ≤ C σ s-5+3 σ s g(σ) H s ν ≤ C σ 2 R since Q t,s,ν (g) ≤ R. This yields (5.10) t 0 R(σ) H s-6 ν-1 ≤ CRh t 0 1 σ 2 dσ ≤ CRh. (1)
Consequently, we get that

sup t∈[0,T ] δ(t) H s-6 ν-1 ≤ C(R)h(1 + (ε + εh)R 2 1 )(1 + hR 1 ).
This ends the proof of Proposition 5.4.

5.3.

Estimate of N T,s-2,ν-1 (δ).

Proposition 5.5. Assuming that η ∈ H s+4 satisfies the assumption (H), there exists C > 0, ε 0 > 0 and h 0 > 0 such that for every h ∈ (0, h 0 ], every ε ∈ (0, ε 0 ] and every T > 0, every solution of (5.1) such that Q T,s-2,ν-1 (δ) ≤ hR 1 enjoys the estimate

N T,s-2,ν-1 (δ) ≤ C(R)h(1 + (ε + εh)R 2 1 )(1 + hR 1 )
where C(R) is a number that depends only on R.

Proof. We proceed as in the previous proof. Let D be the linear operator defined as the Fourier multiplier by k p ξ q ∂ m ξ for (m, p, q) ∈ AE 3d such that p + q ≤ s -2, and m ≤ ν -1. Using (5.8) with this operator D, and Lemma 4.5, we obtain that

δ(t) 2 H s-2 ν-1 ≤ C t 0 ε(m σ,s-1 (d(σ)) δ(σ) H 1 ν-1 δ(σ) H s-2 ν-1 + m σ,2 (d(σ)) δ(σ) 2 H s-2 ν-1
)

+ ε(m σ,s-1 (ζ(σ)) δ(σ) H 1 ν-1 δ(σ) H s-2 ν-1 + m σ,2 (ζ(σ)) δ(σ) 2 H s-2 ν-1 ) + m σ,s-1 (d(σ)) η + εg(σ) H 1 ν-1 δ(σ) H s-2 ν-1 + m σ,2 (d(σ)) η + εg(σ) H s-1 ν-1 δ(σ) H s-2 ν-1 + DR L 2 δ(σ) H s-2 ν-1
dσ.

Using now Q t,s,ν (g) ≤ R and Proposition 5.3, we obtain that

δ(t) 2 H s-2 ν-1 ≤ C t 6 h 2 C(R) 2 (1 + (ε + εh)R 2 1 ) 2 (1 + hR 1 ) 2 + t 3 hC(R)(1 + (ε + εh)R 2 1 ) t 0 R(σ) H s-2
ν-1 dσ. It remains to estimate the last integral to conclude. By using the expression (5.3), we get as in the proof of Lemma 4.5 that

R(σ) H s-2 ν-1 ≤ Ch m (1) σ,s-1 (g(σ))( η H 1 ν + g(σ) H 1 ν ) + m (1) σ,2 (g(σ))( η H s-1 ν + g(σ) H s-1 ν ) with m (1) σ,γ (h) = σ γ sup n=±1 sup ξ∈[nσ,ns(σ)] (| ĥn (ξ)| + |∂ ξ ĥn (ξ)|).
As in the previous proof we have that (5.11) m

(1)

σ,s-1 (g(σ)) ≤ C σ s-1+3 σ s g(σ) H s ν ≤ C σ 2 R since Q t,s,ν (g) ≤ R. This yields (5.12) t 0 R(σ) H s-2 ν-1 ≤ C(R)h t 3
. Collecting all the above estimates, we obtain

sup t∈[0,T ] t -3 δ(t) H s-2 ν-1 ≤ C(R)h(1 + (ε + εh)R 2 1 )(1 + hR 1 ).
This ends the proof of Proposition 5.5.

5.4.

Proof of Proposition 5.1. From Propositions 5.3, 5.4 and 5.5, we have the estimate

Q T,s-2,ν-1 (δ) ≤ C(R)h(1 + (ε + εh)R 2 1 )(1 + hR 1 ), under the assumption Q T,s-2,ν-1 (δ) ≤ hR 1 . Choosing R 1 > C(R), we have C(R)h(1 + (ε + εh)R 2 1 )(1 + hR 1 ) < R 1 h, if ε
and h are small enough. Hence, by an usual continuation argument, the estimate Q T,s-2,ν-1 (δ) ≤ R 1 h is valid for all T > 0, thus proving Proposition 5.1 and with it the convergence estimate 2.6.

PROOF OF THE CONVERGENCE ESTIMATE (2.7) FOR STRANG SPLITTING

From now on, we only consider the case of Strang splitting (1.6), and thus s(t) is given by formula (3.6). Proposition 5.1 implies that, up to the loss of two derivatives, g(t, x, v) is an approximation of order one (with respect to h) of the exact solution of the Vlasov-HMF equation g(t, x, v). Getting the rate of order 2 brings technical complications in order to take advantage of the cancellations provided by the midpoint rule. In view of theorem 3.3 and the main result of [START_REF] Faou | Landau Damping In Sobolev Spaces For The Vlasov-HMF Model[END_REF], we can assume that for all α ∈ {0, 1, 2}, Q T,s,ν-α (g) and Q T,s,ν-α (g) are bounded by the same constant R > 0, provided that ν > 5/2. To prove the result, we proceed as before and start from the equation (5.1) on the error term δ(t). From now on, using the weighted norms defined in (3.8), we consider the quantity:

Q T,s-3,v-2 (δ) = M T,s-4 (d) + N T,s-3,v-2 (δ) + sup t∈[0,T ] δ(t) H s-7 ν-2 .
The convergence result (2.7) will be a consequence of the following proposition: Proposition 6.1. Let us fix s ≥ 9 and ν > 5/2. Assuming that η ∈ H s+4 ν satisfies the assumption (H), there exists R 2 > 0, h 0 > 0 and ε 0 > 0 such that for every h ∈ (0, h 0 ], every ε ∈ (0, ε 0 ] and every T > 0, the solution of (5.1) satisfies the estimate

Q T,s-3,ν-2 (δ) ≤ R 2 h 2 .
The crucial point of the proof of proposition 6.1 is the cancellation provided by the midpoint rule. We can summarize it by the following easy lemma: Lemma 6.2. For t ∈ Ê and a fixed h > 0, let s(t) = s h (t) be given by formula (3.6). Then, for all n ∈ AE, (n+1)h nh (σs(σ))dσ = 0.

Estimate of M T,s-4 (d).

Proposition 6.3. Assuming that η ∈ H s+4 ν satisfies the assumption (H), there exists C > 0, ε 0 > 0 and h 0 such that for every h ∈ (0, h 0 ], every ε ∈ (0, ε 0 ] and every T > 0, every solution of (5.1) such that Q T,s-3,ν-2 (δ) ≤ R 2 h 2 enjoys the estimate

M T,s-4 (d) ≤ C(R)h 2 (1 + (ε + εh 2 )R 2
2 ) Proof. From Equation (5.4) and Lemma 4.3, we still get that M T,s-4 (d) ≤ C (M T,s-4 (G) + M T,s-3 (H)) .

By using the same arguments as in the proof of Proposition 5.3, we can easily obtain that

M T,s-4 (G) ≤ Cε(RR 2 h 2 + h 4 R 2 2
). It thus remains the estimate of M T,s-4 (H) that requires a refined analysis of the cancellations in the integral. By using again the decomposition (5.6), we shall focus on the term H P that is more difficult. By Taylor expanding up to second order, we find that

H P n = t 0 k=±1 k 2 p k (σ -s(σ)) ∂ ξ ĝk (s(σ), ks(σ))ĝ n-k (s(σ), nt -ks(σ)) + ĝk (s(σ), ks(σ))∂ ξ ĝn-k (s(σ), nt -ks(σ)) (ns(σ) -nt) dσ + K n where |K n | ≤ h 2 k=±1 p k t 0 sup α+β 1 ≤2, α =2 α+β 2 ≤2, α =2 sup τ, ξ∈[kσ,ks(σ)] |∂ α t ∂ β 1 ξ ĝk (τ, ξ)||∂ α t ∂ β 2 ξ ĝn-k (τ, nt -ξ)||nσ -nt| dσ.
To estimate the remainder, K n , we can again distinguish the cases k = n and k = -n. Since Q T,s,ν (g) ≤ R, by using the equation (3.4) and Lemma 4.5, we also have that

∂ t g(t) H s-1 ν ≤ C(R) t 2 .
Since ν > 5/2, we can use Lemma 2.3 to obtain, by similar arguments to the ones used in the proof of Proposition 5.3, that

|K n | ≤ h 2 C(R) t 0 1 σ s-4 1 t + σ s-3 (t -σ) dσ + t t/2 1 σ s-4 1 t -σ s-3 (t -σ) dσ and hence that |K n | ≤ 1 t s-4 C(R)h 2 .
To estimate the main term in H P n , assuming that Nh ≤ T ≤ (N + 1)h for some N, we can split the time integral into

t 0 = N -1 j=0 (j+1)h jh + t N h
and we observe that all the integrals (j+1)h jh vanish due to the symmetry of σs(σ) (see Lemma 6.2). We thus obtain that

|H P n | ≤ 1 t s-4 C(R) t N h |t -s(σ)| dσ + 1 t s-4 C(R)h 2 ≤ 1 t s-4 C(R)h 2 .
From the same argument (slightly easier since η does not depend on time and is smoother), we can also prove that

|H L n | ≤ 1 t s-4 C(R)h 2 .
Consequently, by collecting the previous estimates, we find that

M T,s-4 (G) ≤ C(R)h 2 (1 + (ε + εh 2 )R 2
2 ). This ends the proof. 6.2. Estimate of δ H s-7 ν-2 . Proposition 6.4. Assuming η ∈ H s+4 ν satisfies the assumption (H), there exists C > 0, ε 0 and h 0 such that for every h ∈ (0, h 0 ], every ε ∈ (0, ε 0 ] and every T > 0, every solution of (5.1)

such that Q T,s-3,ν-2 (δ) ≤ h 2 R 2 enjoys the estimate sup [0,T ] δ H s-7 ν-2 ≤ C(R)h 2 (1 + (ε + εh 2 )R 2 2 )(1 + h 2 R 2 ).
Proof. We can start as in the proof of Proposition 5.4. By using the energy identity (5.8) with D the Fourier multiplier by k p ξ q ∂ m ξ for (m, p, q) ∈ AE 3d such that p + q ≤ s -7, and m ≤ ν -2 and Lemma 4.5, we first obtain that

Dδ(t) 2 L 2 ≤ C t 0 Dδ(t) L 2 ε(m σ,s-6 (d(σ)) δ(σ) H 1 ν-2 + m σ,2 (d(σ)) δ(σ) H s-7 ν-2 ) + ε(m σ,s-6 (ζ(σ)) δ(σ) H 1 ν-2 + m σ,2 (ζ(σ)) δ(σ) H s-7 ν-2 ) + m σ,s-6 (d(σ)) η + εg(σ) H 1 ν-2 + m σ,2 (d(σ)) η + εg(σ) H s-6 ν-2 dσ + t 0 DR, Dδ(σ) L 2 dσ .

This yields

Dδ(t) 2 L 2 ≤ C(R)h 2 (1 + (ε + εh 2 )R 2 2 )(1 + h 2 R 2 ) t 0 Dδ(σ) L 2 σ 2 dσ + t 0 DR, Dδ(σ) L 2 dσ ≤ C(R)h 2 (1 + (ε + εh 2 )R 2 2 )(1 + h 2 R 2 ) sup σ∈[0,t] δ(σ) H s-7 ν-2 + t 0 DR(σ), Dδ(σ) L 2 dσ . (6.1)
The main difficulty is now to use the cancellation in the midpoint quadrature rule in order to estimate the last integral. Let us define

I(t) = t 0 DR(σ), Dδ(σ) L 2 dσ .
By Taylor expanding, integrating by parts once and using the estimate (5.10), we first write that

I(t) ≤ t 0 DR(σ), Dδ(s(σ)) L 2 dσ + h sup σ∈[0,T ] ∂ t δ(σ) H s-8 ν-2 sup t 0 R H s-6 ν-2 dσ ≤ t 0 DR(σ), Dδ(s(σ)) L 2 dσ + h 2 C(R) sup σ∈[0,T ] ∂ t δ(σ) H s-8 ν-2 .
By using the equation (5.7) and Lemma 4.5, we have that

∂ t δ(t) H s-8 ν-2 ≤ Cm s-7 (d(t)) η H s-7 ν-2 + ε m s-7 (d(t))( δ H s-7 ν-2 + g H s-7 ν-2 ) + m s-7 (ζ) δ H s-7 ν-2
and thus by using Proposition 6.3 we find

∂ t δ(t) H s-8 ν-2 ≤ C(R)h 2 (1 + (ε + εh 2 )R 2 2 ) + εh 4 R 2 2
. This yields

I(t) ≤ t 0 DR(σ), Dδ(s(σ)) L 2 dσ + C(R)h 4 (1 + (ε + εh 2 )R 2
2 ).

Next, by using the definition of R provided in (5.3), we observe that by Taylor expanding in time and in the ξ variable, we can write

(6.2) R(t, x, v) = R 1 (t, x, v) + R 2 (t, x, v)
where R and the remainder R 2 is such that

R 2 (t) H s-7 ν-2 ≤ Ch 2 m (2) t,s-6 (g) η H 1 ν + ε g H 1 ν + m (2) t,2 (g) η H s-6 ν + ε g H s-6 ν with m (2)
t,s-6 (g) ≤ C(R) t 2 . (2)
This yields

R 2 (t) H s-7 ν-2 ≤ C(R) t 2 , and hence (6.3) I(t) ≤ t 0 DR 1 , Dδ(s(σ)) L 2 dσ +C(R)h 4 (1+(ε+εh 2 )R 2 2 )+C(R)h 2 sup σ∈[0,t] δ(σ) H s-7 ν-2 .
Thanks to the definition of R 1 , assuming that Nh ≤ t ≤ (N + 1)h, we obtain that

t 0 DR 1 (σ), Dδ(s(σ)) L 2 dσ = t N h DR 1 (σ), Dδ(s(σ)) L 2 dσ ≤ C(R)h 2 sup σ∈[0,t] δ(σ) H s-7 ν-2 .
Here we have also used the cancellation argument of Lemma 6.2. Consequently, from (6.1), (6.3) and the above estimate, we obtain that

sup σ∈[0,T ] δ(σ) 2 H s-7 ν-2 ≤ C(R)h 4 (1 + (ε + εh 2 )R 2 2 )(1 + h 2 R 2 ) + C(R)h 2 (1 + (ε + εh 2 )R 2 2 )(1 + h 2 R 2 ) sup σ∈[0,T ] δ(σ) H s-7 ν-2 .
By using the Young inequality this yields

sup σ∈[0,T ] δ(σ) H s-7 ν-2 ≤ C(R)h 2 (1 + (ε + εh 2 )R 2 2 )(1 + h 2 R 2 ).
This ends the proof.

6.3. Estimate of N T,s-3,ν-2 (δ).

Proposition 6.5. Assuming η ∈ H s+4 satisfies the assumption (H), there exists C 0, ε 0 > 0, and h 0 such that for every h ∈ (0, h 0 ], every ε ∈ (0, ε 0 ] and every T > 0, every solution of (5.1)

such that Q T,s-3,ν-2 (δ) ≤ h 2 R 2 enjoys the estimate N T,s-3,ν-2 (δ) ≤ C(R)h 2 (1 + (ε + εh 2 )R 2 2 )(1 + h 2 R 2 ).
Proof. We can start once more as in the proof of Proposition 5.4. By using the energy identity (5.8) with D the Fourier multiplier by k p ξ q ∂ m ξ for (m, p, q) ∈ AE 3d such that p + q ≤ s -3, and m ≤ ν -2 and Lemma 4.5, we first obtain that

Dδ(t) 2 L 2 ≤ C t 0 Dδ(t) L 2 ε(m σ,s-2 (d(σ)) δ(σ) H 1 ν-2 + m σ,2 (d(σ)) δ(σ) H s-3 ν-2 ) + ε(m σ,s-2 (ζ(σ)) δ(σ) H 1 ν-2 + m σ,2 (ζ(σ)) δ(σ) H s-3 ν-2 ) + m σ,s-2 (d(σ)) η + εg(σ) H 1 ν-2 + m σ,2 (d(σ)) η + εg(σ) H s-3 ν-2 dσ + t 0 DR, Dδ(σ) L 2 dσ .

This yields

Dδ(t) 2 L 2 ≤ C(R) 2 h 2 (1 + (ε + εh 2 )R 2 2 )(1 + h 2 R 2 ) t 0 σ 2 Dδ(σ) L 2 dσ + t 0 DR, Dδ(σ) L 2 dσ ≤ C(R) 2 t 6 h 4 (1 + (ε + εh 2 )R 2 2 ) 2 (1 + h 2 R 2 ) 2 + t 0 DR(σ), Dδ(σ) L 2 dσ . (6.4)
As in the previous proof we consider

I(t) = t 0 DR(σ), Dδ(σ) L 2 dσ ,
and obtain as previously

I(t) ≤ t 0 DR(σ), Dδ(s(σ)) L 2 dσ + h sup σ∈[0,T ] ∂ t δ(σ) H s-4 ν-2 sup t 0 R H s-2 ν-2 dσ ≤ t 0 DR(σ), Dδ(s(σ)) L 2 dσ + h 2 C(R) t 3 sup σ∈[0,T ] ∂ t δ(σ) H s-4 ν-2 .
Here we have used estimate (5.12). By using the equation (5.7) and Lemma 4.5, we have that

∂ t δ(t) H s-4 ν-2 ≤ Cm s-3 (d(t)) η H s-3 ν-2 + ε m s-3 (d(t))( δ H s-3 ν-2 + g H s-3 ν-2 ) + m s-3 (ζ) δ H s-3 ν-2
and thus by using Proposition 6.3 we find

∂ t δ(t) H s-4 ν-2 ≤ C(R) t 3 h 2 (1 + (ε + εh 2 )R 2 2 ) + t 3 εh 4 R 2 2
. This yields

I(t) ≤ t 0 DR(σ), Dδ(s(σ)) L 2 dσ + t 6 C(R) 2 h 4 (1 + (ε + εh 2 )R 2 2 ) 2 .
As in the previous proof, we write

R(t, x, v) = R 1 (t, x, v) + R 2 (t, x, v)
where R and the remainder R 2 is such that

R 2 (t) H s-3 ν-2 ≤ Ch 2 m (2) t,s-2 (g) η H 1 ν + ε g H 1 ν + m (2) t,2 (g) η H s-2 ν + ε g H s-2 ν with m (2) t,γ (h) = t γ sup n=±1 sup ξ∈[nt,ns(t)] α+β≤2, α =2 |∂ α t ∂ β ξ ĥn (ξ)|. Since ν > 5/2, we have that m (2) t,s-2 (g) ≤ C t 2 g H s-4 ν + t 2 ∂ t g H s-4 ν .
and hence thanks to the equation (3.4) and Theorem 3.3, we obtain that m

t,s-2 (g)C(R) t 2 This yields R 2 (t) H s-3 ν-2 ≤ C(R) t 2 and hence (6.5)

I(t) ≤ t 0 DR 1 , Dδ(s(σ)) L 2 dσ + C(R) 2 t 6 h 4 (1 + (ε + εh 2 )R 2 2 ) 2 + C(R) t 3 h 2 sup [0,t] δ H s-3 ν-2 .
Note that by definition of K(t), the function K(t) is continuous in t = 0, but not C 1 . Using an integration by parts on the definition of the Fourier transform, we then get that (6.10)

|∂ α τ K(τ )| ≤ C τ 2 , α ≤ 2.
To get this, we have used that the function t η0 (t) vanishes at zero. By using this estimate on K, (H) and that F (τ ) ∈ H 1 τ (the Sobolev space in τ ) since F is compactly supported in time, we easily get that y defined via its Fourier transform by (6.8) belongs to H 1 τ . This implies that t y ∈ L 2 and thus that y ∈ L 1 t . These remarks justify the use of the Fourier transform and that the function y defined through its Fourier transform via (6.8) is a solution of (6.6). Moreover, thanks to (6.8) and (H), we get that ŷ can be continued as an holomorphic function in Im τ ≤ 0 and thanks to a Paley Wiener type argument, that y vanishes for t ≤ 0. We have thus obtained an L 1 solution of (6.6) that vanishes for t ≤ 0. By a Gronwall type argument, we easily get that there is a unique solution in this class of (6.6) and thus we have obtained the expression of the unique solution.

We can thus now focus on the proof of the estimate stated in Lemma 4.3. Note that a L 2 -based version of this estimate would be very easily obtained. The difficulty here is to get the uniform L ∞ in time estimate we want to prove.

We shall first prove the estimate for γ = 0. Let us take χ(τ ) ∈ [0, 1] a smooth compactly supported function that vanishes for |τ | ≥ 1 and which is equal to one for |τ | ≤ 1/2. We define χ R (τ ) = χ(τ /R) and χ R (∂ t ) the corresponding operator in t variable corresponding to the convolution with the inverse Fourier transform of χ R (τ ). Thanks to (6.10), we have that for R large To estimate the solution y of (6.6), we shall write that y = χ 2R (∂ t )y + (1 -χ 2R (∂ t ))y =: y l + y h .

By applying (1 -χ 2R (∂ t )) to (6.6), we get that

y h = K * y h + (1 -χ 2R (∂ t ))F = (1 -χ R (∂ t )K * y h + (1 -χ 2R (∂ t ))F
since (1 -χ R ) = 1 on the support of 1 -χ 2R . Therefore, we obtain thanks to (6.11) and the fact that χ 2R (∂ t ) is a convolution operator with a L 1 function, that

y h L ∞ ≤ 1 2 y h L ∞ + C F L ∞
and hence y h L ∞ ≤ 2C F L ∞ . For the low frequencies, we can use directly the form (6.6) of the equation: We can write that ŷl (τ ) = χ 2R (τ )

1 -K(τ ) χ R (τ ) F (τ ).

Since the denominator does not vanish thanks to (H), we obtain again that y l can be written as the convolution of an L 1 function -which is the inverse Fourier transform of χ 2R (τ )/(1 -K(τ )) -by the function χ R (∂ t )F which is a convolution of F by a smooth function. Thus we obtain by using again the Young inequality that y l L ∞ ≤ C F L ∞ . Since y L ∞ ≤ y l L ∞ + y h L ∞ , we get the desired estimate for γ = 0. To get the estimate for arbitrary γ, we can proceed by induction. We observe that ty(t) = K * (ty) + F 1 with F 1 = (tK) * y + tF . Using the result γ = 0, we obtain that ty L ∞ ≤ C F 1 L ∞ . Now since η 0 ∈ H γ+3 , for γ = 1, we obtain that tK ∈ L 1 and thus

F 1 L ∞ ≤ C tF L ∞ + y L ∞ ) ≤ C (1 + t)F L ∞ .
The higher order estimates follow easily in the same way.

Proof of Lemma 4.5. We give the proof of (4.11), the proof of the second estimate being slightly easier. In the Fourier side, we have for L σ [g](h) the expression

(F L σ [g]h) n (ξ) = k∈{±1} kp k Z k (σ) ĥn-k (σ, ξ -ks(σ))(ns(σ) -ξ).
Consequently, we obtain that F ([D r,p,q , L σ [g]]h) ) n (ξ) = k∈{±1} kp k Z k (σ) n p ξ q ∂ r ξ ĥn-k (σ, ξ -ks(σ))(ns(σ) -ξ) -(n -k) p (ξs(σ)) q ∂ r ξ ĥn-k (σ, ξ -ks(σ))(ns(σ) -ξ) . For k = ±1, we can thus expand the above expression into a finite sum of terms under the form I k n (σ, ξ) = kp k Z k (s(σ))k p 1 (n -k) p-p 1 +α ks(σ) q 1 +α (ξ -ks(σ)) q-q 1 +β ∂ r 1 ξ ĥn-k (σ, ξ -ks(σ)) where 0 ≤ p 1 ≤ p, 0 ≤ q 1 ≤ q, m -1 ≤ r 1 ≤ m, α + β = r 1 -m + 1, α, β ≥ 0.

Moreover, if r 1 = r, then we have p 1 + q 1 > 0.

We have to estimate n ξ | k∈±1 I k n (σ, ξ)| 2 dξ by isometry of the Fourier transform. We note that for a fixed k ∈ {±1} then for |n -k| + |ξ -ks(σ)| ≤ |k|s(σ), we have This ends the proof of the Lemma.

4. 4 .

 4 Proof of Theorem 3.3 and estimate (2.5). The proof of Theorem 3.3 follows from the a priori estimates in Propositions 4.1, 4.4 and 4.6 and a continuation argument. Indeed, by combining the estimates of these three propositions, we get that

t 2 |( 1 -≤ 1 2

 212 (1 -χ R (∂ t ))K(t)| ≤ C α≤2 ∂ α τ ((1 -χ R (τ )) K(τ )) L 1 (R) ≤ C χ R (∂ t ))K(t) L 1 (R) ≤ CR for R sufficiently large. This choice fixes R.

  |I k n (σ, ξ)| ≤ Cs(σ) p+q+1 |Z k (s(σ))||n -k||∂ r 1 ξ ĥn-k (σ, ξ -ks(σ))| whereas for |n -k| + |ξ -ks(σ)| ≥ |k|s(σ), we have |I k n (σ, ξ)| ≤ C σ 2 |Z k (s(σ))|(|n -k| + |ξ -ks(σ)|) γ |∂ r 1 ξ ĥn-k (σ, ξ -ks(σ))|.Consequently by taking the L 2 norm, we find that k∈±1I k n (σ, ξ) L 2 ≤ C m σ,γ+1 (Z) h(σ) H 1 + m σ,2 (Z) h(σ) H m .

This work is partially supported by the ERC starting grant GEOPARDI No. 279389.

Thanks to the definition of R 1 , assuming that Nh ≤ t ≤ (N + 1)h, we obtain that

Once more we have used the cancellation argument (see lemma 6.2). Consequently, from (6.4), (6.5) and the above estimate, we obtain that

By using the Young inequality this yields

This ends the proof.

6.4. Proof of proposition 6.1. Using Propositions 6.3, 6.4 and 6.5, we have proven that

We can thus obtain Proposition 6.1 with the same bootstrap argument as before by choosing R 2 > C(R) and then by taking ε and h sufficiently small. APPENDIX Proof of Lemma 4.3. Let us first note that the equation (4.6) only involves K(n, t) for positive t, hence K(n, t) can be replaced here by K 1 (n, t) (see (2.3)). Let us take T > 0, and let us set for the purpose of the proof K(t) = K 1 (n, t), F (t) = (F n (t) + G n (t))1 0≤t≤T . Since we only consider the cases n = ±1, we do not write down anymore explicitly the dependence in n. We consider the equation (6.6)

setting y(t) = 0 for t ≤ 0. Note that the solution of this equation coincides with ζ n (t) on [0, T ] since the modification of the source term for t ≥ T does not affect the past. By taking the Fourier transform in t (that we still denote by • ), we obtain

Under the assumption (H), the solution of (6.7) is given explicitely by the formula (6.8)

Let us observe that since (1 + v 2 )η 0 ∈ H 5 , we have by (2.8) that for α ≤ 2 and for t > 0 (6.9)


M Gallo

K Ege

M Rébillat
email: marc.rebillat@ensam.eu

J Antoni

A MULTI-SINE SWEEP METHOD FOR THE CHARACTERIZATION OF WEAK NON-LINEARITIES; PLANT NOISE AND VARIABILITY ESTIMATION

Weak non-linearities in vibrating structures can be characterized by a signal-model approach based on cascade of Hammerstein models. The experiment consists in exciting a device with a sine sweep at different levels, in order to assess the evolutions of non linearities on a wide frequency range. The method developed in [1], based on exponential sine sweep, is able to give an approximative identification of the Hammerstein models, but cannot make the distinction between nonlinear distortion and stationary plant noise. Therefore, this paper proposes improvements on the method that provide a more precise estimation of the Hammerstein models through the cancellation of the plant noise: it relies on the repetition of the signal on a certain amount of periods (multi-sine sweeps) and then on the consideration of the synchronous average out of the different periods from the resulting signal. Mathematical foundations and practical implementation of the method are discussed. The second main point of improvement concerning the study of the vibrating device is the use of the Bootstrap analysis. By considering some periods randomly chosen among the multisine sweep, one can study the variability of the experiments. The method becomes more robust.

INTRODUCTION

Many structures are generally assumed to undergo linear vibrations when they are driven at low levels, as the transverse motion w remains in a smaller range than the plate thickness, which is in practice not perfectly the case. Also, as non linearities appear -to a certain point -at a lower amplitude than the linear response, it seems fundamental to measure theses non linearities with precision. Therefore the aim is to estimate measurement noise, whether it come from electronics or instruments, and then extract it out of the resulting signal.

This paper is thus organized around on the improvements concerning the exponential sine sweep based method, mainly developed in the work of Rébillat et al. [START_REF] Rébillat | Identification of cascade of hammerstein models for the description of nonlinearities in vibrating devices[END_REF], in order to study non linear behavior of vibrating devices. The method at stake is described in Sec. 2. Then, the noise estimation is developed in Sec. [START_REF] Novak | Modeling of non linear audio systems using swept-sine signals: application to audio effects[END_REF]. Concerning the improvement of the method itself, a more precise estimation of the kernels of Hammerstein models is discussed in Sec. [START_REF] Gallman | Iterative method for identification of non linear -systems using a uryson model[END_REF]. Also, the uncertainty measurement is presented through the Bootstrap procedure in Sec. [START_REF] Farina | Simultaneous measurement of impulse response and distorsion with a swept-sine technique[END_REF]. The obtained results on a vibrating plate are finally given in Sec. 6.

THE SINE SWEEP METHOD

Mathematical foundations

Cascade of Hammerstein models

Volterra series are a convenient tool to express analytically the relationship between the input e(t) and the output s(t) of a weakly non-linear system [START_REF] Boyd | Fading memory and the problem of approximating nonlinear operators with volterra series[END_REF] which is fully characterized by the knowledge of its Volterra Kernels in the frequency domain

{V k (f 1 , • • • , f k )} k∈N * .
Cascade of Hammerstein models constitute an interesting subclass of Volterra systems whose Kernels possess the following property:

∀k, ∃ Ṽk : ∀(f 1 , • • • , f k ), V k (f 1 , • • • , f k ) = Ṽk (f 1 + • • • + f k ) (1)
Volterra Kernels of cascade of Hammerstein models can thus be expressed as functions of only one frequency variable and are in practice easy to estimate experimentally [START_REF] Rébillat | Identification of cascade of hammerstein models for the description of nonlinearities in vibrating devices[END_REF][START_REF] Novak | Modeling of non linear audio systems using swept-sine signals: application to audio effects[END_REF]. This simple method is based on a phase property of exponential sine sweeps and the Kernels of such a model can be estimated from only one measured response of the system.

In a cascade of Hammerstein models [START_REF] Gallman | Iterative method for identification of non linear -systems using a uryson model[END_REF], each branch is composed of one non linear static polynomial element followed by a linear one h n (t), as shown fig. 1. Mathematically, the relation between the input e(t) and the output s(t) of such a system is given by Eq. 2 ,where * denotes the convolution:

s(t) = N n=1 (h n * e) n (t) (2)
In this model, each impulse response h n (t) is convolved with the input signal elevated to its n th power and the output s(t) is the sum of these convolutions. The first impulse response h 1 (t) represents the linear response of the system. The other impulse responses {h n (t)} n∈{2•••N } model the non linearities. The family {h n (t)} n∈{1•••N } will be referred to as the Kernels of the model. These Kernels are assumed to be integrable. Any cascade of Hammerstein models is fully represented by its Kernels.

As can be seen in Eq. 2, cascade of Hammerstein models are linear in the parameters to be estimated, i.e. the output of the system is a linear combination of the Kernels

{h n (t)} n∈{1•••N } .
The method studied here gives direct mathematical access to all the Kernels {h n (t)} n∈{1•••N } from the contributions of the different orders of non linearity obtained as in Ref. [START_REF] Farina | Simultaneous measurement of impulse response and distorsion with a swept-sine technique[END_REF]. The main advantage of the proposed approach is that it provides a direct evaluation of the N impulse responses h n (t) of the system. The foundations and the key implementation of this simple method are explained in detail in [START_REF] Rébillat | Identification of cascade of hammerstein models for the description of nonlinearities in vibrating devices[END_REF], where the whole procedure is validated on a simulated system and on two real systems. Estimating each kernel h n (t) of a cascade of Hammerstein models is not a straightforward task. A simple estimation method that has been proposed previously by the authors [START_REF] Rébillat | Identification of cascade of hammerstein models for the description of nonlinearities in vibrating devices[END_REF] for this purpose is briefly recalled here. To experimentally cover the frequency range over which the system under study has to be identified, cosines with time-varying frequencies are commonly used. When the instantaneous frequency of e(t) = cos[f (t)] is increasing exponentially from f 1 to f 2 in a time interval T , this signal is called an "Exponential Sine Sweep". It can be shown [START_REF] Rébillat | Identification of cascade of hammerstein models for the description of nonlinearities in vibrating devices[END_REF][START_REF] Novak | Modeling of non linear audio systems using swept-sine signals: application to audio effects[END_REF] that by choosing

T m = (2mπ -π/2)ln( f 2 f 1 )/2πf 1 with m ∈ N * one obtains: ∀k ∈ N * , cos(kφ(t + ∆t k )) with ∆t k = T m lnk ln( f 2 f 1 ) (3)
Eq. 3 states that for any exponential sine sweep of duration T , multiplying the phase by a factor k yields to the same signal, advanced in time by ∆t k .

Kernel recovery in the time domain

If an exponential sine sweep is presented at the input of a cascade of Hammerstein models, by combining Eq. 2 and Eq. 3 we obtain the following relation:

s(t) = N n=1 (γ n * e) n (t + ∆t n ) with γ n (t) = N n=1 c k,n h k (t) (4)
where g n (t) is the contribution of the different kernels to the n th harmonic. In order to identify each kernel h n (t) separately, a signal y(t) operating as the inverse of the input signal e(t) in the convolution sense can be built as shown in [START_REF] Rébillat | Identification of cascade of hammerstein models for the description of nonlinearities in vibrating devices[END_REF]. After convolving the output of the cascade of Hammerstein models s(t) given in Eq. 4 with y(t), one obtains Eq. ( 4):

(y * s)(t) = N n=1 γ n (t + ∆t n ) (5)
Because ∆t n ∝ ln(n) and f 2 > f 1 , the higher the order of non-linearity n, the more advanced is the corresponding γ n (t). Thus, if T m is chosen long enough, the different g n (t) do not overlap in time and can be separated by simply windowing them in the time domain. Using Eq. 5, the family {h n (t)} n∈{1•••N } of the kernels of the cascade of Hammerstein models under study can then be fully extracted. Details of the computation of the matrix C are provided in [START_REF] Rébillat | Identification of cascade of hammerstein models for the description of nonlinearities in vibrating devices[END_REF].

[h 1 (t) • • • h N (t)] T = C[γ 1 (t) • • • γ N (t)] T (6)

An original method for detection of non linearities

The main advantage of this sine sweep based method is to be quick and simple; one can have direct access to the kernels of a vibrating device. However, this method needed to be improved to some extent. The main problem was the impossibility to make the difference between non linearities and the experimental noise, which resulted in an overestimation of the non linear behavior. Also, in order to fit better to the mathematical model, the number of kernels taken into account to describe the behavior of the vibrating device had to be more precise. Finally, in order to study the vibrating plate through the "Sweep" method, a simple tool was implemented to the previous method to provides the user with an indicator of the variability on the experimentations.

NOISE ESTIMATION

The multi-sine sweep described here is basically a signal containing the repetition of the exact same sweep. The sweep is defined by a starting frequency, an ending frequency and a time laps.

In the sweep-based method developed in [START_REF] Rébillat | Identification of cascade of hammerstein models for the description of nonlinearities in vibrating devices[END_REF], the input signal is a single sweep. The main hypothesis relies on the fact that the experimental noise follows a stochastic pattern. Therefore the aim to take the average out of the multi sine sweep from the result signal.

Mathematical foundations

The signal can be defined as the sum of response of the system s and the experimental noise n :

x(t) = s(t) + n(t) (7)
For a number K of exact same sweep of period T , the synchronous average corresponds to:

s(t) = 1 K K-1 k=0 [x(t -kT )] (8)
The noise can then be estimated by:

n(t) = x(t) -s(t) (9)
Finally, the signal to noise ratio (SNR) is:

SNR = 20 • log 10 ( RM S(s) RM S(n) ) (10)

Experimental validation

We can see on figure 2 the SNR levels of the result of different experiments, in comparison with the theoretical curve:

y = SNR(2) + 3 • ln 2 (2/N ) (11)
Experimental curves are likely to follow the same pattern as the theoretical one: a decrease of 3 dB as one double the number of repetitions.

Noise canceling and SNR Indicator

Once the noise is extracted from the result signal, it becomes then possible to assess the signal to noise ratio (SNR) of the experiment. This ratio gives one a precise indication on the validity of the experiments. The signal applied to the vibrating device must not be too low, in order to avoid the amplitude of the experimental noise reaching the amplitude of the global signal. However, a great SNR ratio also means important displacement, which involve chaotic behavior from the plate, instead of being weakly non linear; the model would become false.

SETTING OF THE RIGHT NUMBER OF KERNELS

The previous method could estimate a certain number of kernels. This number had to be defined by the user and could be wrong. Actually, a too small number of kernels taken into account resulted in a severe misunderstanding of the non linearities of the vibrating device. This is due to the fact that the lower kernels are depending on the higher ones. On the other side, a too big number of kernels implicated to consider experimental noise as a kernel. Therefore, the procedure has to be autonomous, in its way of assessing the right number of kernels called Nopt. Hence the idea of having the computer estimating a great number of kernels, and to compare the RMS level of each of the estimated kernels. Once the procedure detects that the RMS level of a certain kernel do not reach the RMS level of experimental noise, this kernel must not be taken into account, and the method only study kernels from lower orders.

UNCERTAINTY OF MEASUREMENT

Statistical analysis of the experience

For a multi sine sweep of K periods, one does B random draws with replacement of K periods among the original K periods. For each value b (b = 1...B) the procedure calculates the synchronous average over the K periods randomly chosen. Then the B signals are implemented in the procedure in order to obtain different kernels. One gets then B number of each of the kernels estimated.

The idea then is to evaluate the accuracy of the experiments by calculating the synchronous average for each kernel, and then the confidence interval.

A more robust method

By synthetically increasing the number of experiments, one can assess the device behavior in a more reasonable time, as the Bootstrap implementation is not a really time consuming. Thus, this procedure provides the user with a more accurate view on the non linear behavior of the device under study. The device under study is a 1mm steel plate (see fig. 3), in free boundaries conditions. A shaker glued to the plate excites the plate with a multi sine sweep. In order to damp the system, a foam sheet was added to the plate.

APPLICATION TO A VIBRATING PLATE

Noise estimation and canceling

Fig. 4 shows how the noise perturbs the signal. This comparison between few random periods chosen among the result signal and the synchronous averaged signal gives a precise view on the spectral density of the stochastic noise; it is mainly present in the mid-frequencies.

Also, canceling the noise from the signal provides one with the certainty that resonance peaks really are a vibrating issue, and not a low quality power supply problems. The cancellation of the noise in the experimental signal is particularly important for the assessment of high orders kernels, which are very sensitive to noise due to there low level.

SNR indicator and setting of the right number of kernels

As depicted on fig. 5 non linearities level increases so as the global gain of the input signal. From a very low gain where the plate barely vibrates to great displacements (higher than the plate's thickness) where non linearities level reaches linear level, one can observe several things concerning the improvements of the method. First, SNR is increasing as long as the input gain increases. In this case, SNR should be between 35 dB and 45 dB, so that the non linear response is comparable to the linear response, considering a weakly non linear behavior. Secondly, the method is now able to adapt to any situation, taking the right number of kernels into account. This number is getting bigger so as the input gain. Figure 5: Hammerstein kernels for different amplitudes. Nopt is the optimal number of kernels calculated by the method.

Uncertainty of measurement

One can see on fig. 6 the application of the Bootstrap analysis on an experimental signal. Thanks to the Bootstrap procedure, the user is now able to detect artefacts in the experimental signal, particularly in mid-frequencies. This point is very important in the case of experiments in noisy environment.

CONCLUSION AND PERSPECTIVES

In this article an original vibro-acoustical method is presented to assess the proportion of nonlinearities among the resulting signal. This method, based on a sine sweep, is able to assess whether the vibrating device behaves linearly or not. Moreover, the repetition of the same sweep signal considerably improved the accuracy of the method, as it is now possible to extract noise from the resulting signal, and then to estimate the signal to noise ratio (SNR). The Bootstrap analysis really enhanced the view on the frequency response of the vibrating plate. The improvements on the method developed by Rébillat et al. [START_REF] Rébillat | Identification of cascade of hammerstein models for the description of nonlinearities in vibrating devices[END_REF] now provides the user with a more precise tool to study vibration and acoustic response of a vibrating device. The experiments being now K times longer (with K the number of periods of the multi sine sweep) than before, a next step would be to find a compromise between accuracy and time consuming.

Figure 1 :

 1 Figure 1: Block diagram, representation of a cascade of Hammerstein models

Figure 2 :

 2 Figure 2: Multi sine sweep (--) compared to the theoretical curve (-) the SNR ratio correspond here to the ratio between the RMS amplitude of the signal on the extracted RMS amplitude of the noise.

Figure 3 :

 3 Figure 3: Steel plate excited with a shaker.

Figure 4 :

 4 Figure 4: Frequency response of the system (linear part) without noise (-) (synchronous averaging) and with noise (-):1 period randomly chosen from the signal.

Figure 6 :

 6 Figure 6: Frequency response of the system (linear part) with 10 Bootstrap samples.

ACKNOWLEDGMENTS

Since 2011, the Laboratoire Vibrations Acoustique is part of the LabEx CeLyA of Université de Lyon, operated by the French National Research Agency (ANR-10-LABX-0060/ ANR-11-IDEX-0007).


Sylvain Lobry

Loïc Denis

Sparse + smooth decomposition models for multi-temporal SAR images

SAR images have distinctive characteristics compared to optical images: speckle phenomenon produces strong fluctuations, and strong scatterers have radar signatures several orders of magnitude larger than others. We propose to use an image decomposition approach to account for these peculiarities.

Several methods have been proposed in the field of image processing to decompose an image into components of different nature, such as a geometrical part and a textural part. They are generally stated as an energy minimization problem where specific penalty terms are applied to each component of the sought decomposition.

We decompose temporal series of SAR images into three components: speckle, strong scatterers and background. Our decomposition method is based on a discrete optimization technique by graph-cut. We apply it to change detection tasks.

I. INTRODUCTION

Denoising and restoration is one of the oldest challenges of image processing. Like other coherent imaging techniques, SAR imagery suffers from strong fluctuations due to speckle phenomenon. Speckle can be described as a multiplicative noise, while most algorithms developed in the field of image processing are designed to deal with an additive noise. This calls for an adaptation of algorithms developed for optical images or the development of new techniques.

The simplest way to reduce the noise is multi-looking which amounts to averaging pixel values within a small window. With a window big enough, noise can be strongly reduced, at the cost of a resolution loss. This kind of method is not well suited to areas with bright scatterers. Point-like strong scatterers have a high contrast compared to the background and are typically caused by man-made structures, such as buildings. Local averaging methods need to deal with these points separately in order to prevent from spreading them out on the whole window. Strong scatterers can be identified using for instance a likelihood ratio test as in [START_REF] Lopes | Adaptive processing of multilook complex SAR images[END_REF] and [START_REF] Lopes | Structure detection and statistical adaptive speckle filtering in SAR images[END_REF]. The corresponding detector compares the radiometry of the point under test to the average value computed over a surrounding window (e.g., cross-shaped). While being simple, this test fails when other scatterers are present in the window.

To cope with these problems related to the specific nature of the SAR images, we propose to model separately these two contributions (background and bright scatterers) and to jointly estimate them. This paper extends the idea of SAR image decomposition introduced in [START_REF] Aujol | Image decomposition application to SAR images[END_REF] and in [START_REF] Denis | Exact discrete minimization for TV+L0 image decomposition models[END_REF] to temporal series and applies it to change detection. This decomposition is then applied to detect changes between SAR images taken at different times. Detecting and analyzing changes is one of the primary application for SAR data, mainly thanks to the reliability of acquisition regarding the weather. The problem has first been tackled by [START_REF] Singh | Digital change detection techniques using remotely-sensed data[END_REF] using a difference operator. State-of-the art methods now use likelihood ratio test such as in [START_REF] Su | Change detection and classification of multi-temporal SAR series based on generalized likelihood ratio comparing-and-recognizing[END_REF] and in [START_REF] Quin | Mimosa: An automatic change detection method for sar time series[END_REF].

In section II, we describe a model to detect strong scatterers and perform change detection. This model is stated as a minimization problem. In section III, we provide a way of obtaining the exact solution for this problem. Finally, we describe and evaluate our models in an application of change detection in section IV.

II. MULTI-TEMPORAL SPARSE + SMOOTH DECOMPOSITION MODEL

Total variation (TV) minimization [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] is a method widely used for noise reduction. It strongly penalizes signal fluctuations while preserving edges by minimizing the total variation of the signal (i.e., the sum of the spatial gradient magnitude).

Although this method has been applied by several authors to speckle noise reduction in SAR imagery (see [START_REF] Aubert | A variational approach to removing multiplicative noise[END_REF]- [START_REF] Palsson | SAR image denoising using total variation based regularization with SUREbased optimization of the regularization parameter[END_REF]), it is not well adapted to urban areas where punctual bright scatterers are frequent. Indeed, bright scatterers are points with a radiometric value several orders of magnitude higher than the surrounding background. Using total variation regularization for such signals has the effect of removing the point-like scatterers and of biasing the estimation of the neighboring area (spreading effect). On the other hand, these points are sparsely distributed in the image and are thus present in a limited number. In this section, we propose a decomposition model combining regularity of the background and sparsity of the bright scatterers in multi-temporal series.

We consider a temporal series of observed images V = {v 1 , . . . , v n } where v t is the observed image at the t-th date. The underlying (noiseless) scene is U = {u 1 , . . . , u n }. The scene u t at the t-th date is modeled as a sum of 2 components: u BVt (component with bounded variations, representing the background) and u St (sparse component representing the bright scatterers). The observed image v t is related to the underlying scene u t by:

v t = (u BVt + u St ) × n t , (1)
where n t is a speckle noise realization corresponding to a SAR acquisition system with multiplicative noise model. In the case of SAR images, noise follows a Rayleigh distribution (when considering amplitude images) and is weakly correlated (if images are not over-sampled) so that the likelihood distribution can be considered separable:

p(V|U) = n t=1 i∈Ω p(v t (i)|u t (i)) = n t=1 i∈Ω 2 v t (i) u 2 t (i) exp -v 2 t (i) u 2 t (i) , ( 2
)
where Ω is the set of the pixels in an image and u t (i) is the value at pixel i in image t.

We define next the prior p(U) in order to enforce the decomposition of the radar scene into its two independent components U BV and U S :

-log p(U) = -β BV log p(U BV ) -β S log p(U S ) , (3)
where the distributions p(U BV ) and p(U S ) model each component. The background is consider to be formed by homogeneous regions with sharp boundaries (piece wise constant model), with possibly some abrupt changes between two dates. Images following such models have a low 3D total variation, we therefore choose to define:

-log p(U BV ) = TV 3D (U BV ) = n t=1   i,j∈C |u BVt (i) -u BVt (j)|   + n-1 t=1 i∈Ω |u BVt (i) -u BVt+1 (i)| ,
where C is the set of all cliques in Ω depending of the chosen neighborhood. The choice of an anisotropic version of total variation (the gradient magnitude is defined in the L1 norm sense) will find its justification in the numerical optimization method described in the next section. The component containing strong scatterers is sparse, which is classically enforced by an energy that increases with the number of non-zero elements, such as the L0 pseudo-norm:

-log p(U S ) = U S 0 . (4)
Up to a constant term, the (opposite) log posterior distribution corresponds to the following energy:

E(U) = -log(p(V|U)) -log(p(U)), (5) with:
-

log(p(V|U)) = -log n t=1 i∈Ω 2 v t (i) u 2 t (i) exp -v 2 t (i) u 2 t (i) = n t=1 i∈Ω -log(2v t (i)) + 2 log(u BVt (i) + u St (i)) + v 2 t (i) (u BVt + u St ) 2 (i) = DT(v t , u BVt , u St ), (6)
with the first term (-log(2v t (i))) independent of U, which will thus be discarded in the following. In summary, the energy is given by:

E(U) = n t=1 β S L0(u St ) + β BV TV 3D (U BV ) + DT(v t , u BVt , u St ) . (7)

III. EXACT DISCRETE OPTIMIZATION BY GRAPH-CUTS

The cost function introduced in equation 7 is highly nonconvex. In this section, we show how to find the global optimum of the optimization problem (up to a given precision defined by the quantization step size) without resorting to an approximation by L1 relaxation of the L0-norm, as typically done when solving such optimization problem. First, we consider the sub-problem of solving arg min

U S E(U)
for a fixed U BV :

U S (U BV ) = arg min U S n t=1 (DT(v t , u BVt , u St ) + β S L0(u St )) + β BV TV 3D (U BV ) .
Since all terms are separable, we can solve the minimization independently for each pixel:

u St (i) (u BVt (i)) = u St (i) if DT(v t (i), u BVt (i), u St (i) ) +β S < DT(v t (i), u BVt (i), 0) 0 otherwise With u St (i) = arg min u S t (i) DT(v t (i), u BVt (i), u St (i)).
We can now rewrite the minimization problem under a form involving only U BV :

arg min u BV E(U) = arg min u BV n t=1 DT(v t , u BVt , u St (u BVt )) + β S u S (u BV ) 0 + β BV TV 3D (U BV ) (8)
The first two terms of equation 8 are separable (i.e., they correspond to a sum of terms involving a single pixel at a time) and the last one is convex and involves only pairwise terms (i.e., pairs of pixel values). This minimization problem can thus be exactly minimized by graph-cuts using the method described in [START_REF] Ishikawa | Exact optimization for Markov random fields with convex priors[END_REF]. This optimization method consists if seen as a discrete optimization problem, i.e., for a given quantization of the component u BV in searching the min-cut in a graph where every possible value for each pixel is represented by a node. Neighboring nodes (spatially or in time) are interconnected by vertices with weights corresponding to β BV . Nodes representing the same pixel i at two consecutive possible values α and α + 1 are interconnected by vertices with weights p(v t (i)|α). We show an example of a decomposition obtained with our model using this optimization technique in figure 1.

Note that the graph construction used here presents a high memory usage: it requires V = n × |Ω| × |A| vertices and E = 7n × |Ω| × |A| edges where n is the number of images in the series, |Ω| is the number of pixels in each image and |A| the number of possible values for the background image. In the implementation of the min-cut algorithm introduced in [START_REF] Boykov | An experimental comparison of mincut/max-flow algorithms for energy minimization in vision[END_REF], a vertex takes 48 bytes and an edge takes 32 bytes. Therefore, the graph construction requires 272 × n × |Ω| × |E| bytes, limiting the size of the images that can be processed. On a larger image, we could use algorithms such as those described in [START_REF] Denis | SAR image regularization with fast approximate discrete minimization[END_REF] or in [START_REF] Shabou | A graph-cut based algorithm for approximate MRF optimization[END_REF] that sub-sample the set of amplitudes to lower the memory occupation of the optimization technique and obtain an approximate solution to the problem.

Regarding computational cost, it has a worst-case complexity of O(EV 2 |C|) where |C| is the minimum cut and is also the complexity of our method. On a computer with an Intel Xeon(R) CPU E5-1620 with 16Gb of RAM, the algorithm takes 52,04s to compute on 2 images of 300 × 400 pixels with 50 levels of quantification. Note that we do not fully benefit from the power of the processor as the implementation of the algorithm is single-core.

IV. SCATTERERS CHANGE DETECTION

A. Method

For a series of observed images V, the model described in section II provides a serie of decompositions {(u BV1 , u S1 ), . . . , (u BVn , u Sn )}. A first change detection approach consists of exploiting the scatterer detection to find changes between them. This is not a straightforward task since the value and the exact position of the detected scatterers can greatly vary when a small change happens in the scene. Therefore, we start by binarizing the scatterers images {u St } The main difficulty to overcome when performing change detection based on point-like scatterers is that these points are not necessarily stable. To handle this situation, we look for pixels that are detected as scatterers in an image without any scatterer detected in a neighborhood in the other image. This can be formalized by:

T (i) = δ u S bin 1 (i + δ) - δ u S bin 2 (i + δ) (9)
where T is a temporary image, and T (i) represents the chances of having a change at position i. δ is used to iterate over the neighborhood of the pixel i. This image is then thresholded to find the changes between the two images.

B. Results

We ran this algorithm on a time serie of SAR images of Saint-Gervais (France). We compare the results of this algorithm with other change detection methods in figure 2 and we show an example of results in image 3. While the results are not as good as dedicated state-of-the-art change detection algorithms such as the ones presented in [START_REF] Su | NORCAMA: Change Analysis in SAR Time Series by Likelihood Ratio Change Matrix Clustering[END_REF] our method performs better than the ones presented in [START_REF] Lombardo | Maximum likelihood approach to the detection of changes between multitemporal SAR images[END_REF] and in [START_REF] Krylov | Change detection with synthetic aperture radar images by Wilcoxon statistic likelihood ratio test[END_REF].

V. CONCLUSION

This paper introduces a decomposition model suitable to multi-temporal series of SAR images. The proposed model combines a TV regularization with an L0 pseudo-norm to achieve a good estimation of the radiometry of the scene even when point-like bright scatterers are present. It provides a decomposition between the background of the scene and the scatterers. We proposed a simple application of this decomposition framework to change detection using the images of scatterers computed by the algorithm on a SAR series. The method yields results close to that of state-of-the-art algorithms.

Further work includes more complex change detection approaches using all components (background and targets). A contrario framework could be used to detect changes as in [START_REF] Dellinger | Change detection for high resolution satellite images, based on SIFT descriptors and an a contrario approach[END_REF] from the scatterers images. Furthermore, strategies should be investigated to tackle the practical problem caused by the heavy consumption of memory of our optimization method.

Fig. 1 :

 1 Fig. 1: Decomposition of an image of Saint-Gervais acquired by TerraSAR-X. Thanks to the German Aerospace Agency (DLR) for the images (project MTH0232 and LAN1746).

  and Oliver, 2001] Wilcoxon Test [Krylov et al., 2012] R GLRT [Su et al. 2014] The proposed TV+L0 change detection

Fig. 2 :

 2 Fig. 2: False positive alarm versus true positive curves of various change detection algorithms

  (a) Input image 1. (b) Input image 26. (c) Changes on image 1. (d) Changes on image 26.

Fig. 3 :

 3 Fig.3: Change detection results using the proposed method on images of Saint-Gervais series. Regions with changes that have been detected are indicated in red.

This work has been funded by the CNES and the "Future et Ruptures" program.