N Forcadel
email: forcadel@cermics.enpc.fr

An error estimate for a new scheme for mean curvature motion

Keywords: AMS Classification: 65M06, 65M12, 65M15, 49L25, 53C44 Mean curvature motion, error estimate, numerical scheme, dislocations dynamics

In this work, we propose a new numerical scheme for the anisotropic mean curvature equation. The solution of the scheme is not unique, but for all numerical solutions, we provide an error estimate between the continuous solution and the numerical approximation. This error estimate is not optimal, but as far as we know, this is the first one for mean curvature type equation. Our scheme is also applicable to compute the solution to dislocations dynamics equation.

Introduction

Mean curvature motion has been largely studied in terms of both theory and computation, in particular due to the large number of applications like front propagation, image processing, fluid dynamics...(see for instance Sethian [START_REF] Sethian | Level set methods and fast marching methods[END_REF] and Osher, Paragios [START_REF]Geometric level set methods in imaging, vision, and graphics[END_REF]).

The level set framework has been used in both theoretical and numerical problems. Theoretically, the mean curvature equation has been well understood using the framework of viscosity solutions by Chen, Giga, Goto [START_REF] Chen | Uniqueness and existence of viscosity solutions of generalized mean curvature flow equations[END_REF] and Evans, Spruck [START_REF] Evans | Motion of level sets by mean curvature. I[END_REF]. However, this equation has serious problems for the question of numerical approximations. Nevertheless, there are several works on this question. First, let us mention the numerical method of Osher, Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF]. This method is very used in practice but, as far as we know, there is no convergence result. Another algorithm is the Merriman, Bence, Osher scheme [START_REF] Merriman | Diffusion generated motion by mean curvature[END_REF] in which motion by mean curvature is viewed as singular limit of a diffusion equation with threshold. The convergence of this scheme has been proved by Barles,Georgelin [4] and Evans [START_REF] Evans | Convergence of an algorithm for mean curvature motion[END_REF] (see also Ishii [16], Ishii, Pires, Souganidis [START_REF] Ishii | Threshold dynamics type approximation schemes for propagating fronts[END_REF], and Chambolle, Novaga [START_REF] Chambolle | Approximation of the anosotropic mean curvature flow[END_REF]). A class of convergent schemes for nonlinear parabolic equations including mean curvature motion have been proposed by Crandall, Lions [START_REF] Crandall | Convergent difference schemes for nonlinear parabolic equations and mean curvature motion[END_REF]. Let us mention also the recent work of Oberman [START_REF] Oberman | A convergent monotone difference scheme for motion of level sets by mean curvature[END_REF]. In this last two works, two different scales are used. The first one is the space step ∆x and the second one is the size of the stencil ε. As it was point out in [START_REF] Crandall | Convergent difference schemes for nonlinear parabolic equations and mean curvature motion[END_REF], these two scales are very important to approximate degenerate equations like mean curvature equation.

The goal of our work is to propose a new scheme for mean curvature motion and to prove an error estimate between the continuous solution and its numerical approximation. This error estimate is not optimal, but as far as we know, this is the first one for complete discretized scheme for mean curvature type motion.

The idea is to use a recent work of Da Lio, Monneau and the author [START_REF] Da Lio | Convergence of a non-local eikonal equation to anisotropic mean curvature motion. Application to dislocations dynamics[END_REF] concerning the convergence of dislocations dynamics to mean curvature motion. Dislocations are linear defects which move in crystals. Their dynamics can be represented by a non-local first order Hamilton-Jacobi equation (see Alvarez, Hoch, Le Bouar and Monneau [START_REF] Alvarez | Dislocation dynamics: short time existence and uniqueness of the solution[END_REF]). The first goal of this work is to prove an error estimate between the solution of dislocations dynamics and the solution of mean curvature motion. To do this, we will use in a more quantitative way the definition of viscosity solutions for mean curvature motion proposed by Barles, Georgelin [START_REF] Barles | A simple proof of convergence for an approximation scheme for computing motions by mean curvature[END_REF] by considering a regularization with quartics.

The second goal of this work is to propose a numerical scheme for dislocations dynamics. The main properties of this scheme is that it is implicit (so there is no CFL condition) and uses two different scales. Moreover, this scheme is not monotone and does not admit a unique solution. It is only "almost monotone" (see Lemma 5.2). Nevertheless, the fact that it is implicit and the monotonicity of the velocity, allows us to "freeze" the velocity and to prove a Crandall Lions type [START_REF] Crandall | Two approximations of solutions of Hamilton-Jacobi equations[END_REF] error estimate for any solutions (we refer to Alvarez, Carlini, Monneau and Rouy [START_REF] Alvarez | Convergence of a first order scheme for a non local eikonal equation[END_REF][START_REF]A convergent scheme for a nonlocal hamilton-jacobi equation, modeling dislocation dynamics[END_REF] for error estimate for dislocations dynamics in the non monotone case). It is also possible to explicit the scheme. In this case, we are able to prove the convergence of the scheme under a CFL condition, but error estimates in this case are still open. This comes from the fact that there is no consistency error for the scheme (see Proposition 5.6).

Finally, let us mention some works on error estimate for numerical scheme for Hamilton-Jacobi equations. For the first order case, we refer to Souganidis [START_REF] Souganidis | Approximation schemes for viscosity solutions of Hamilton-Jacobi equations[END_REF]. For the second order case, when the Hamiltonian is convex, we refer to Krylov [START_REF] Krylov | On the rate of convergence of finite-difference approximations for Bellman's equations[END_REF][START_REF]On the rate of convergence of finite-difference approximations for Bellman's equations with variable coefficients[END_REF] and Barles, Jakobsen [START_REF] Barles | On the convergence rate of approximation schemes for Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Barles | Error bounds for monotone approximation schemes for Hamilton-Jacobi-Bellman equations[END_REF]. When the Hamiltonian is uniformly elliptic, we refer to Caffarelli, Souganidis [START_REF] Caffarelli | A rate of convergence for monotone finite difference approximations to fully nonlinear uniformly elliptic pde[END_REF].

Let us now explain how this paper is organized. In Section 2, we state the main results of this work. In Section 3, we prove the main result concerning mean curvature type motion. Section 4 is devoted to prove the error estimate between the solution of dislocations dynamics and the one of mean curvature motion. In Section 5, we study the numerical scheme for dislocations dynamics and we prove an error estimate between the continuous solution and its numerical approximation. Some numerical simulations are provided in Section 6.

Main Results

Error estimate for Mean Curvature Motion

In order to propose a numerical scheme for anisotropic mean curvature motion, we will use the work of Da Lio et al. [START_REF] Da Lio | Convergence of a non-local eikonal equation to anisotropic mean curvature motion. Application to dislocations dynamics[END_REF] that we briefly recall here. Given a function g defined on the unit sphere

S N -1 of R N by (2.1) g ∈ Lip(S N -1 ), g(-θ) = g(θ) ≥ 0, ∀θ ∈ S N -1 we consider kernels c 0 ∈ L ∞ (R N ) ∩ W 1,1 (R N ) satisfying (2.2)        c 0 (x) = 1 |x| N +1 g x |x| if |x| ≥ 1, c 0 (-x) = c 0 (x) ≥ 0, ∀x ∈ R N
and we use the following scaling for 0 < ε < e -1

(2.3)

c ε 0 (x) = 1 ε N +1 | ln ε| c 0 x ε .
The term | ln ε| comes from the bad decay at infinity of the kernel c 0 (see [START_REF] Da Lio | Convergence of a non-local eikonal equation to anisotropic mean curvature motion. Application to dislocations dynamics[END_REF]Section 4.1]).

We then consider the following auxiliary problem (2.4)

       u ε t (x, t) = (c ε 0 ⋆ 1 {u ε (•,t)≥u ε (x,t)} )(x) - 1 2 R N c ε 0 |Du ε (x, t)| in R N × (0, T ), u ε (•, 0) = u ε 0 in R N
where u ε t denotes the derivative with respect to the time variable, Du ε indicates the gradient of u ε with respect to the space variables, the convolution is done in space only and 1 {u ε (•,t)≥u ε (x,t)} is the characteristic function of the set {u ε (•, t) ≥ u ε (x, t)} (which is equal to 1 on the set and 0 outside). This equation arises in the theory of dislocations dynamics (see Alvarez, Hoch, Le Bouar and Monneau [START_REF] Alvarez | Dislocation dynamics: short time existence and uniqueness of the solution[END_REF]) and uses the Slepčev formulation [START_REF] Slepčev | Approximation schemes for propagation of fronts with nonlocal velocities and Neumann boundary conditions[END_REF] to consider the simultaneous evolutions of all the level sets of the function u ε . We refer to Da Lio et al. [START_REF] Da Lio | Convergence of a non-local eikonal equation to anisotropic mean curvature motion. Application to dislocations dynamics[END_REF] for the study of this model and in particular to [START_REF] Da Lio | Convergence of a non-local eikonal equation to anisotropic mean curvature motion. Application to dislocations dynamics[END_REF]Definition 2.1] for the definition of viscosity solutions for Problem (2.4). In particular, we recall that we take the indicatrice of {u ε (•, t) ≥ u ε (x, t)} for sub-solution and the indicatrice of {u ε (•, t) > u ε (x, t)} for super-solution.

The main result of [START_REF] Da Lio | Convergence of a non-local eikonal equation to anisotropic mean curvature motion. Application to dislocations dynamics[END_REF] is that, if u ε 0 = u 0 , then the unique solution u ε of (2.4) converges locally uniformly on compact sets to the solution u 0 of the following limit equation

(2.5)        u 0 t -trace D 2 u 0 • A Du 0 |Du 0 | = 0 in R N × (0, T ) u 0 (•, 0) = u 0 in R N with (2.6) A p |p| = θ∈ S N -2 = S N -1 ∩{ x, p |p| =0} 1 2 g(θ) θ ⊗ θ dθ
and where M • A and •, • denote respectively the product between the two matrices and the usual scalar product. Our first main result is an error estimate between u ε and u 0 :

Theorem 2.1 (Error estimate for Mean Curvature Motion) Let N ≥ 1 and T ≤ 1. Assume that u 0 , u ε 0 ∈ Lip(R N ), g ∈ Lip(S N -1 ) and that c 0 given in (2.2) satisfies c 0 ∈ W 1,1 (R N ). Then, there exists a constant K 1 depending on N, sup R N c 0 , |Dg| L ∞ (S N -1 ) and |Du 0 | L ∞ (R N )
such that the difference between the solution u ε of (2.4) and the solution u 0 of (2.5) is given by

(2.7) sup R N ×(0,T ) |u ε -u 0 | ≤ K 1 T | ln ε| 1 6 + sup R N |u ε 0 -u 0 |.

Discrete-continuous error estimate for dislocations dynamics

To approximate the solution u 0 of (2.5), we then have to approximate the solution u ε of (2.4). Up to a change of variable (see Corollary 3.1), it suffices to approximate the solution u of (2.8)

       u t (x, t) = (c 0 ⋆ 1 {u(•,t)≥u(x,t)} )(x) - 1 2 c 0 |Du(x, t)| in R N × (0, T ), u(•, 0) = ū0 in R N .
Given a mesh size ∆x 1 , ..., ∆x N , ∆t and a lattice

Q ∆ T = Q ∆ × {0, ..., (∆t)N T } where Q ∆ = {(i 1 ∆x 1 , ..., i N ∆x N ), I = (i 1 , ..., i N ) ∈ Z N }
and N T is the integer part of T /∆t, we will denote by (x 1 , ..., x N , t n ) = (x I , t n ) the node (i 1 ∆x 1 , ..., i N ∆x N , n∆t) and by v n I the value of a numerical approximation of the exact solution u(x I , t n ). We set ∆x = ∆x 2 1 + ... + ∆x 2 N the space mesh size. We shall assume throughout that ∆x + ∆t ≤ 1.

The discrete solution v is computed iteratively by solving the implicit scheme (2.9)

v 0 I = ũ0 (x I ), v n+1 I -v n I ∆t = c ∆ [v] n+1 I G(v n+1 ) I
where ũ0 is an approximation of ū0 and G(v n+1 ) I is a suitable approximation of the gradient of v n+1 taken at point x I . The non-local velocity is the discrete convolution

(2.10) c ∆ [v] n+1 I = c ∆ [v](x I , t n+1 ) = J∈Z N c0 I-J 1 {v n+1 J ≥v n+1 I } ∆x 1 ...∆x N - 1 2 J∈Z N c0 J ∆x 1 ...∆x N with (2.11) c0 I = 1 |Q I | QI c 0 (x) dx
where Q I is the square cell centred in x I

(2.12)

Q I = [x i1 -∆x 1 /2, x i1 + ∆x 1 /2] × ... × [x iN -∆x N /2, x iN + ∆x N /2].
Finally, let us define

(2.13) v # (y, t n ) = I v(x I , t n )χ QI (y)
where χ QI is the indicator function of Q I . The approximation of the gradient is obtained using the Osher, Sethian scheme [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF] (we can also use the one proposed by Rouy, Tourin [START_REF] Rouy | A viscosity solutions approach to shape-from-shading[END_REF]). It is monotone, consistent and depends on the sign of the non-local velocity. Its precised definition is recalled in Section 5.

Since the velocity c ∆ is non-local and not continuous, we have to give a sense to the equality in the scheme (2.9). In fact, we will use the analogue of the Slepčev formulation [START_REF] Slepčev | Approximation schemes for propagation of fronts with nonlocal velocities and Neumann boundary conditions[END_REF] for discrete sub and super-solution (see Definition 5.1) and we will use a discrete version of the Perron's method to construct a discrete solution. The solution of the scheme is not unique but, for any solutions, we have the following Crandall-Lions type [START_REF] Crandall | Two approximations of solutions of Hamilton-Jacobi equations[END_REF] error estimate: Theorem 2.2 (Discrete-continuous error estimate for (2.8)) Let N ≥ 1 and T ≤ 1. Assume that ū0 , ũ0 ∈ W 1,∞ (R N ) and that c 0 given in (2.2) satisfies c 0 ∈ W 1,1 (R N ). Assume that ∆x + ∆t ≤ 1. Then there exists a constant K 2 > 0 depending only on N ,

|c 0 | W 1,1 (R N ) , |Dū 0 | L ∞ (R N ) and |Dũ 0 | L ∞ (R N )
such that the error estimate between the continuous solution u of (2.8) and any discrete solution v of the finite difference scheme (2.9) is given by sup

R N ×{0,...,tN T } |u -v # | ≤ K 2 √ T (∆x + ∆t) 1/2 + sup R N |ū 0 -(ũ 0 ) # | provided ∆x + ∆t ≤ 1 K 2 2 .
Remark 2.3 It is possible to explicit the computation of the gradient, i.e., to replace the term G(v n+1 ) I by G(v n ) I in the scheme (2.9) and to consider the solution v of (2.14)

v 0 I = ũ0 (x I ), v n+1 I -v n I ∆t = c ∆ [v] n+1 I G(v n ) I .
In this case, as usual, we have to satisfy a CFL condition like for instance

∆t ≤ ∆x 2|c 0 | L 1 (R N )
for the Osher Sethian discretisation of the gradient. Under this additional assumption, Theorem 2.2 remains true with v solution of the scheme (2.14).

Discrete continuous error estimate for Mean Curvature Motion

Using the above results, we will prove in Section 3 the following theorem: Theorem 2.4 (Discrete-continuous error estimate for the mean curvature motion) Let N ≥ 1 and T ≤ 1. Let us denote by v ε a solution of (2.9)-(2.10)-(2.11) with initial condition ũε 0 (which is an approximation of u 0 ) and with c ε 0 in the place of c 0 . Assume that u 0 , ũε 0 ∈ W 1,∞ (R N ), g ∈ Lip(S N -1 ) and that c 0 given in (2.2) satisfies c 0 ∈ W 1,1 (R N ). Assume also that ∆x + ∆t ≤ 1. Then there exists a constant

K 3 > 0 depending only on N, sup R N c 0 , |Dg| L ∞ (S N -1 ) , |c 0 | W 1,1 (R N ) , |Du 0 | L ∞ (R N ) and |Dũ ε 0 | L ∞ (R N )
such that the error estimate between the continuous solution u 0 of (2.5) and its numerical approximation v ε is given by sup Proposition 5.4 for the exact setting), we can iterate the process and get a linear estimate in

R N ×{0,..,tN T } |u 0 -v ε # | ≤ K 3 T | ln ε| 1/6 + sup R N |u 0 -(ũ ε 0 ) # | where ε ≥ K 3 (∆x + √ ∆t). Remark 2.5 If T ≥ 1, since |Du 0 | L ∞ (R N ×(0,T )) ≤ |Du 0 | L ∞ (R N ) and "|Dv ε | L ∞ (R N ) "≤ |Dũ ε 0 | L ∞ (R N ) (see
T , i.e sup R N ×{0,..,tN T } |u 0 -v ε # | ≤ K 3 | ln ε| 1/6 T + sup R N |u 0 -(ũ ε 0 ) # |.
Remark 2.6 We can truncate the kernel c 0 at infinity and consider

cR 0 (x) = c 0 (x) if |x| ≤ R, 0 else
In this case, we make an error of order

R N \BR(0) c 0 ≤ K R
and we can make the computation on a finite stencil, even if ∆x goes to zero. This is possible, if we choose ε of the same order of ∆x. The condition ε ≥ K(∆x + √ ∆t) in Theorem 2.4 then implies that we need to impose a CFL condition ∆t ≤ K∆x 2 (which is classical for second order equation).

At the opposite, if we do not impose any CFL condition, we can choose ∆t larger than ∆x 2 , but we have to choose ε of the same order of √ ∆t and so to make the convolution on larger and larger stencils as ∆x goes to zero.

As we mention in the introduction, it is also possible to explicit the computation of the velocity in the scheme (2.9) and to consider (2.15)

v 0 I = ũ0 (x I ), v n+1 I -v n I ∆t = c ∆ [v] n I G(v n ) I ,
with c ε 0 in the place of c 0 in the definition of c ∆ . We can prove the convergence of the scheme (see Theorem 2.7) under the CFL condition

(2.16) ∆t ≤ ε| ln ε| 2|c 0 | L 1 (R N ) ∆x.
But, in this case, we are not able to prove an error estimate. The reason is that when we implicit the velocity, it "freezes" it, and then we can use the consistency error on the scheme with given velocity. At the opposite, for the explicit scheme, we have to control the consistency error of the velocity which is not possible as we point out in Proposition 5.6. We define

v ε,δ (x, t) = v ε (x I , t n ) if x ∈ Q I , t ∈ [t n , t n+1 )
where δ = (∆x, ∆t) and v ε is the solution of the scheme (2.15)-(2.10)-(2.11) with the kernel c ε 0 in the place of c 0 . Then we have the following convergence result: Theorem 2.7 (Convergence of the explicit scheme) Assume (2.16). Under the assumptions of Theorem 2.4, the function v ε,δ converges locally uniformly on compact sets as δ → 0 to u ε solution of (2.4).

Remark 2.8 In Theorem 2.7, if we take the limit δ → 0 and ε → 0 with δ ≪ ε, we will approach the solution of (2.5). The condition δ ≪ ε implies that, in practice, we have to make the convolution in larger and larger stencil as δ → 0. Remark 2.9 In the scheme (2.15), it is possible to implicit the computation of the gradient and thus to withdraw the CFL condition (2.16) to get the same result as Theorem 2.7.

Notation In what follows, we will denote by K a generic constant, which will then satisfy K + K = K, K • K = K, and so on.

3 Numerical scheme for mean curvature motion

Proof of Theorem 2.4

Using a rescaling argument, we will prove the following corollary of Theorem 2.2: Corollary 3.1 (Discrete-continuous error estimate for (2.4)) Let us denote by v ε a solution of (2.9)-(2.10)-(2.11) with initial condition ũε 0 (which is an approximation of u ε 0 ) and with c ε 0 in the place of c 0 . Under the assumptions of Theorem 2.2, there exists a constant K > 0

depending only on |c 0 | W 1,1 (R N ) , |Du ε 0 | L ∞ (R N ) and |Dũ ε 0 | L ∞ (R N )
such that the error estimate between the continuous solution u ε of (2.4) and its numerical approximation v ε is given by sup

R N ×{0,...,tN T } |u ε -v ε # | ≤ K 1 ε T | ln ε| ε∆x + ∆t | ln ε| + sup R N |u ε 0 -(ũ ε 0 ) # | provided ε∆x + ∆t | ln ε| ≤ ε 2 K 2 .
Proof of Corollary 3.1 First we remark that a simple change of variable gives (for u and v solutions respectively of (2.8) and (2.9))

u ε (y, τ ) = εu y ε , τ ε 2 | ln ε| , u ε 0 (y) = εū 0 y ε v ε (y, τ ) = εv y ε , τ ε 2 | ln ε| , ũε 0 (y) = εũ 0 y ε . Let us denote by x = y ε , t = τ ε 2 | ln ε| , ∆x = ∆y ε , ∆t = ∆τ ε 2 | ln ε| , T = Γ ε 2 | ln ε|
, s n = n∆τ and by N Γ the integer part of Γ/∆τ . Then, the following inequality holds: sup

R N ×{0,...,sN Γ } |u ε -v ε # | =ε sup R N ×{0,...,tN T } |u -v # | ≤K √ Γ | ln ε| ∆y ε + ∆τ ε 2 | ln ε| + ε sup R N |ū 0 -(ũ 0 ) # | ≤K √ Γ ε | ln ε| ε∆y + ∆τ | ln ε| + sup R N |u ε 0 -(ũ ε 0 ) # | provided that ∆y ε + ∆τ ε 2 | ln ε| ≤ 1 K 2 , i.e. ε∆y + ∆τ | ln ε| ≤ ε 2 K 2 .
This ends the proof of the corollary.

Proof of Theorem 2.4

The proof of Theorem 2.4 is now very easy. Indeed, using Theorem 2.1 and Corollary 3.1, we deduce that

|u 0 -v ε | ≤|u 0 -u ε | + |u ε -v ε | ≤K T | ln ε| 1 6 + K ε |lnε| √ T √ ε∆x + ∆t + sup R N |u ε 0 -(ũ ε 0 ) # | + sup R N |u ε 0 -u 0 |. Taking ε ≥ K(∆x + √ ∆t)
implies the prevalence of the first term with respect to the second one and so this implies the desired estimate for the choice of initial condition u ε 0 = u 0 .

Proof of Theorem 2.7

The idea of the proof is borrowed from Barles, Souganidis [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF]. Let us set

v ε = lim sup δ→0 * v ε,δ , v ε = lim inf δ→0 * v ε,δ .
We now prove that v ε and v ε are respectively sub and super-solution of (2.4

). Let ϕ ∈ C ∞ (R N × [0, T ))
such that v ε -ϕ reaches a strict maximum at the point (x 0 , t 0 ) with t 0 ∈ (0, T ). Then, there exists δ → 0, (x δ , t δ ) → (x 0 , t 0 ) such that v ε,δ -ϕ reaches a maximum at (x δ , t δ ). We denote by (

x δ I , t δ n ) ∈ Q ∆ T the node such that v ε,δ (x δ , t δ ) = v ε,δ (x δ I , t δ n ). This implies that, for all (x, t): v ε,δ (x, t) ≤ v ε,δ (x δ I , t δ n ) -ϕ(x δ , t δ ) + ϕ(x, t).
Using the monotony of the scheme for given velocity and the fact that v is solution of the scheme, we get

ϕ(x δ , t δ ) -ϕ(x δ , t δ -∆t) ∆t ≤c ε [v ε,δ ](x δ I , t δ -∆t) G + (D + ϕ(x δ , t δ -∆t), D -ϕ(x δ , t δ -∆t)) if c[v ε,δ ](x δ I , t δ -∆t) ≥ 0 G -(D + ϕ(x δ , t δ -∆t), D -ϕ(x δ , t δ -∆t)) if c[v ε,δ ](x δ I , t δ -∆t) ≤ 0 with c ε [u](x, t) = c ε 0 ⋆ 1 {u(•,t)≥u(x,t)} (x) - 1 2 R N c ε 0
and where we have used Lemma 5.5 to obtain the velocity. Sending δ → 0, using Slepčev Lemma [26, equation ( 5)] and the consistency of the scheme for given velocity, we get

ϕ t (x 0 , t 0 ) ≤ c ε 0 ⋆ 1 {v ε (•,t0)≥v ε (x0,t0)} (x 0 ) - 1 2 R N c ε 0 |Dϕ(x 0 , t 0 )|.
This proves that v ε is a sub-solution of (2.4). The proof for v ε is the same and we skip it. Moreover, using a barrier argument (using the equivalent of Proposition 5.4) we get that (2.4). This ends the proof of the theorem. Remark 3.2 In the scheme (2.15), we have made the choice to take the velocity c ∆ [v] n I . An another possibility is to take

v ε (•, 0) = v ε (•, 0) = u ε 0 . Then we have that v ε = v ε = v ε is solution of
c∆ [v] n+1 I = J∈Z N c0 I-J 1 {v n+1 J >v n+1 I } ∆x 1 ...∆x N - 1 2 J∈Z N c0 J ∆x 1 ...∆x N .
This construct two discrete solution v and ṽ. So, we can take any discrete function comprised between ṽ and v. In fact, this is equivalent to define a discrete solution as a sub and a super-solution as in Definition 5.1. With this definition, v will be the greater sub-solution and ṽ will be the lower super-solution and Theorem 2.7 will be true for every solutions.

What happens if we change the kernel?

A natural question is what happens if we change the kernel and if we take a kernel which decrease more quickly at infinity. For this kind of kernel K 0 , the natural scaling is the following one (see [13, Section 4.1])

K ε 0 (x) = 1 ε n+1 K 0 x ε .
In fact, using the same arguments as the one we use for the proof of Theorem 2.1, we can prove the same kind of error estimate. For example for K 0 (x) = 1 |x| n+p for |x| ≥ 1, with p ≥ 3, we get

|u ε -u 0 | ≤ K √ εT 1/4 .
This is the best estimate for ε small that we can obtain for general kernel, as we can see in Step 3, Case 1 of the proof of Theorem 2.1.

The main difference is that the estimate of Corollary 3.1 is replaced by sup

R N ×{0,...,tN T } |u ε -v ε # | ≤ K √ T ε √ ε∆x + ∆t + sup R N |u ε 0 -(ũ ε 0 ) # | for ε∆x + ∆t ≤ ε 2 K 2 .
Finally, we obtain for the choice

u ε 0 = u 0 sup R N ×{0,...,tN T } |u 0 -(v ε ) # | ≤ K √ εT 1 4 + sup R N |u 0 -ũε 0 | if ε ≥ K √ ∆x + (∆t) 1 3
. This implies that ∆x ε → 0 as ∆x → 0. So, when ∆x goes to zero, we have to make the convolution on all the space which can be very expensive in practice. This approach can also be used for Gaussian kernel and should give the same estimates. In particular, this could give an error estimate for the equivalent version of the classical Bence, Merriman, Osher [START_REF] Merriman | Diffusion generated motion by mean curvature[END_REF] algorithm.

Proof of Theorem 2.1

Before to prove Theorem 2.1, we need some notation. Let us define

F (M, p) = trace M • A p |p| and G(M, p) = -1 |p| F (M, p)
where we recall that A(p) • p = 0 ∀p ∈ S N -1 . Then we have the following fundamental estimate for balls:

Lemma 4.1 (Error estimate for a ball) Let ϕ ∈ C 2 (R N ) with Dϕ(x 0 ) = 0, be such that the set {ϕ(x) ≥ ϕ(x 0 )} is a ball of radius R. For

c ε 0 (•) = 1 ε n+1 | ln ε| c 0 • ε , let us define c ε = (c ε 0 ⋆ 1 {ϕ(•)>ϕ(x0)} )(x 0 ) - 1 2 R n c ε 0 .
Then, there exists a constant

K = K(N, sup R n c 0 , |Dg| L ∞ (S N -1 ) ) > 0 such that for 0 < ε < δ with 0 < δ ≤ R/2, we have |c ε + G(D 2 ϕ(x 0 ), Dϕ(x 0 ))| ≤ K • e (ε, δ, R) with e (ε, δ, R) = 1 | ln ε| 1 δ + 1 R | ln δ| + δ R 2 .
This is a straightforward consequence of Da Lio et al. [START_REF] Da Lio | Convergence of a non-local eikonal equation to anisotropic mean curvature motion. Application to dislocations dynamics[END_REF]Proposition 4.1] where their Lemma 4.3 is replaced by a direct estimate for an explicit function h

(x ′ ) = R -R 2 -|x ′ | 2 with x ′ ∈ R N -1 and |x ′ | ≤ R.
Proof of Theorem 2.1 Let u 0 be the solution of the mean curvature motion (2.5). The idea of the proof is to regularise the function u 0 by a kind of sup-convolution but using quartic penalization and then to plug the regularised function u 0,α into (2.4). This regularisation allows us to control quantitatively the first and the second derivatives in space of u 0,α . This is a quantitative version of the definition of Barles and Georgelin [START_REF] Barles | A simple proof of convergence for an approximation scheme for computing motions by mean curvature[END_REF] for viscosity solutions of mean curvature motion where they used this kind of regularisation to prove that one can take test functions such that Dϕ = 0 or Dϕ = 0 and D 2 ϕ = 0 in the definition of viscosity solutions. This kind of arguments is also used to obtained the comparison principle for mean curvature type equations.

The proof is decomposed into five steps:

Step 1. Regularisation of u 0 By classical estimates for mean curvature type equations, we have (4.17)

Du 0 L ∞ (R N ×(0,T )) ≤ Du 0 L ∞ (R N ) .
We regularise u 0 by considering a spatial sup-convolution by quartics

u 0,α (x, t) = sup y∈R N u 0 (x -y, t) - 1 4α |y| 4 (4.18) = sup z∈R N u 0 (z, t) - 1 4α |z -x| 4 . (4.19)
Since u 0 is Lipschitz continuous, the supremum is reached.

Step 2. Estimate between u 0 and u 0,α For x ∈ R N and t ∈ (0, T ), we denote by x = x(x, t) a point where the maximum is reached in (4.19). Then (4.20)

u 0,α (x, t) = u 0 (x, t) - 1 4α |x -x| 4 ≥ u 0 (x,
|u 0,α (x, t) -u 0 (x, t)| ≤ Kα 1/3
Before to continue the proof, we need some notation. We now want to compare u ε and u 0,α . To do that, we will prove that u 0,α ± (f (ε, α)t + Kα

1 3 + sup R N |u ε 0 -u 0 |)
are respectively super and sub-solutions of Problem (2.4) for

f (ε, α) = K 1 | ln ε| + ln | ln ε| α 1/3 | ln ε| + 1 α 2/3 1 | ln ε| + ε . Let ψ ∈ C ∞ (R N × (0, T )
) such that u 0,α -ψ reaches a global strict maximum at (x 0 , t 0 ) ∈ R N × (0, T ). Since u 0,α is semi-convex in the space variable, the functions u 0,α is derivable with respect to x at the point (x 0 , t 0 ).

Let us denote by

ϕ z (x) = u 0 (z, t 0 ) - 1 4α |z -x| 4 . Then u 0,α (x, t 0 ) = sup z∈R N ϕ z (x).
We denote by x0 a point where the maximum is reached for x = x 0 and by ϕ = ϕ x0 . We then deduce that (4.25) ψ(x, t 0 ) ≥ u 0,α (x, t 0 ) ≥ ϕ(x)

with equality for x = x 0 . This implies that

e(ϕ) = |c ε [ϕ](x 0 )|Dϕ(x 0 )| -F (D 2 ϕ(x 0 ), Dϕ(x 0 ))| where c ε [ϕ](x 0 ) = c ε 0 ⋆ 1 {ϕ(•)≥ϕ(x0)} (x 0 ) - 1 2 R N c ε 0 .
Step 3. e(ϕ) ≤ f (ε, α)

We distinguish three cases:

Case 1. h ≤ 2ε. In this case, we directly have the estimate (using (4.29))

|c ε [ϕ](x 0 )|Dϕ(x 0 )|| ≤ |c ε 0 | L 1 (R N ) h 3 α ≤ h 2 |c 0 | L 1 (R N ) ε| ln ε| h α ≤ K α 2/3 ε | ln ε|
where we have also used (4.28). Using moreover the fact that

|F (D 2 ϕ(x 0 ), Dϕ(x 0 ))| ≤ K h |Dϕ(x 0 )| ≤ K h 2 α ,
we then deduce that e(ϕ) ≤ K α 2/3 ε ≤ f (ε, α).

Case 2. 2ε < h ≤ 2

| ln ε| . Using (4.29) and Lemma 4.1, with δ = h/2, we then deduce that

e(ϕ) ≤ Kh 3 α 1 | ln ε| 1 h + 1 h | ln δ| + 1 h ≤ Kh α 2/3 ≤ K α 2/3 1 | ln ε| ≤ f (ε, α). Case 3. h > 2 | ln ε| . We set δ = 1 | ln ε|
. We then have ε < δ < h/2. So by Lemma 4.1, we deduce that

e(ϕ) ≤ Kh 3 α 1 | ln ε| + 1 h ln | ln ε| | ln ε| + 1 h 2 1 | ln ε| ≤K 1 | ln ε| + 1 α 1/3 ln | ln ε| | ln ε| + 1 α 2/3 1 | ln ε| ≤f (ε, α).
Step 4. Estimate between u 0,α and u ε

By construction, we have, since u 0,α is a sub-solution of (2.5) (but does not satisfy the initial condition) ψ t (x 0 , t 0 ) ≤F (D 2 ψ(x 0 , t 0 ), Dψ(x 0 , t 0 ))

≤c ε [ϕ](x 0 )|Dψ(x 0 , t 0 )| + F (D 2 ϕ(x 0 ), Dϕ(x 0 )) -c ε [ϕ](x 0 )|Dϕ(x 0 )| ≤c ε [u 0,α (•, t 0 )](x 0 )|Dψ(x 0 , t 0 )| + f (ε, α)
where we have used (4.25) for the last line. On the other hand, we have

u 0,α (•, 0) -Kα 1/3 -sup R N |u ε 0 -u 0 | ≤ u ε 0
where we have used (4.24) at the limit t = 0. We then deduce that ũ0,α (x, t) = u 0,α (x, t) -f (ε, α)t - (2.4). This implies that ũ0,α ≤ u ε and so

Kα 1/3 -sup R N |u ε 0 -u 0 | is sub-solution of
(4.30) u 0,α -u ε ≤ f (ε, α)t + Kα 1/3 + sup R N |u ε 0 -u 0 |.
Step 5. Estimate between u 0 and u ε Using (4.24) and (4.30), we then get

u 0 -u ε =u 0 -u 0,α + u 0,α -u ε ≤ f (ε, α)t + Kα 1/3 + sup R N |u ε 0 -u 0 |. But f (ε, α)t + Kα 1/3 =K 1 | ln ε| + ln | ln ε| α 1/3 | ln ε| + 1 α 2/3 1 | ln ε| + ε t + Kα 1/3 ≤K T | ln ε| 1 3
for α = Kt | ln ε| and t ≤ T ≤ 1. We then get

u 0 -u ε ≤ K T | ln ε| 1 3 + sup R N |u ε 0 -u 0 |.
A lower bound is proved similarly (using the same result as Lemma 4.1 where the set {ϕ(x) ≥ ϕ(x 0 )} is replaced by {ϕ(x) > ϕ(x 0 )}). This implies (2.7) for T ≤ 1. This ends the proof of the theorem.

5 Numerical scheme for dislocations dynamics

Definitions and preliminary results

We recall here the notation used in the scheme. The discrete solution v is computed iteratively by solving the implicit scheme (5.31)

v 0 I = ũ0 (x I ), v n+1 I -v n I ∆t = c ∆ [v] n+1 I G(v n+1 ) I
where the non-local velocity is defined in (2.10). The approximation of the gradient of v n+1 at the point x I is given by

G(v n+1 ) I = G + (D + v n+1 I , D -v n+1 I ) if c ∆ [v] n+1 I ≥ 0, G -(D + v n+1 I , D -v n+1 I ) if c ∆ [v] n+1 I < 0
where G ± is a suitable approximation of the Euclidean norm and

D ± v n (x I ) = (D ± x1 v n (x I ), ..., D ± xN v n (x I
)) are the discrete gradients. The terms D ± xi v n (x I ) are the standard forward and backward first order differences, i.e. for a general function f (x I ):

D + xi f (x I ) = f (x I i,+ ) -f (x I ) ∆x i , D - xi f (x I ) = f (x I ) -f (x I i,-) ∆x i ,
where (5.32)

I k,± = (i 1 , ..., i k-1 , i k ± 1, i k+1 , ..., i N ).
The approximations of the Euclidean norm G ± are those proposed by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF] (we can also use the ones proposed by Rouy, Tourin [START_REF] Rouy | A viscosity solutions approach to shape-from-shading[END_REF]):

(5.33)

G + (P, Q) =   i=1,..,N max(P i , 0) 2 + i=1,..,N min(Q i , 0) 2   1 2 , G -(P, Q) =   i=1,..,N min(P i , 0) 2 + i=1,..,N max(Q i , 0) 2   1 2
.

We recall that the functions G ± are Lipschitz continuous with respect to the discrete gradients, i.e.

(5.34)

|G ± (P 1 , P 2 ) -G ± (P ′ 1 , P ′ 2 )| ≤ (|P 1 -P ′ 1 | + |P 2 -P ′ 2 |).
They are consistent with the Euclidean norm

(5.35) G ± (P, P ) = |P |,
and G ± = G ± (P + 1 , ..., P + N , P - 1 , ..., P - N ) enjoy suitable monotonicity with respect to each variable (5.36)

∂G + ∂p + i ≥ 0, ∂G + ∂p - i ≥ 0, ∂G - ∂p + i ≥ 0, ∂G - ∂p - i ≥ 0.
To define discrete solution, we need the following notation

(5.37) c∆ [v] n+1 I = J∈Z N c0 I-J 1 {v n+1 J >v n+1 I } ∆x 1 ...∆x N - 1 2 J∈Z N c0 J ∆x 1 ...∆x N , G(v n+1 ) I = G + (D + v n+1 I , D -v n+1 I ) if c∆ [v] n+1 I ≥ 0, G -(D + v n+1 I , D -v n+1 I ) if c∆ [v] n+1 I < 0
where c0 is defined in (2.11). Finally, for simplicity of presentation, let us denote by

s[v] n+1 I = Sign(c ∆ [v] n+1 I ), s[v] n+1 I = Sign(c ∆ [v] n+1 I ).
Definition 5.1 (Numerical sub, super and solution of the scheme) We say that v is a discrete sub-solution (resp super-solution) of the scheme (2.9) if for all I ∈ Z N , n ∈ N, we have

v n+1 I ≤ v n I + ∆t c ∆ [v] n+1 I G(v n+1 ) I resp. v n+1 I ≥ v n I + ∆t c∆ [v] n+1 I G(v n+1 ) I .
We say that v is a discrete solution if and only if it is a sub and a super-solution.

Lemma 5.2 (Almost monotonicity of the scheme)

The scheme (2.9) is almost monotone in the following sense. Let v, w be two discrete functions and assume that there is I such that w I = v I and w J ≥ v J for J = I. Then

c ∆ [v] I G(v) I ≤ c ∆ [w] I G(w) I and c∆ [v] I G(v) I ≤ c∆ [w] I G(w) I .
Proof of Lemma 5.2 First, we remark that c ∆ [v] I ≤ c ∆ [w] I . Then, there is three cases:

1. c ∆ [v] I ≤ 0 ≤ c ∆ [w] I .
In this case, the result is trivial.

2. 0 ≤ c ∆ [v] I ≤ c ∆ [w] I .
Using the monotonicity of G, we get that G(v) I ≤ G(w) I . This implies the result.

3. c ∆ [v] I ≤ c ∆ [w] I ≤ 0.
Using the monotonicity of G, we get that G(v) I ≥ G(w) I . This implies the result.

The proof for c∆ is the same and we skip it. This ends the proof of the lemma.

The existence of a solution for the scheme is not trivial (since it is implicit and non-local). This is the subject of the following proposition: Proposition 5.3 (Existence of solution for the scheme (2.9)) There exists, at least, one discrete solution v of the scheme (2.9) in the sense of Definition 5.1.

Proof of Proposition 5.3

We assume that there exists a solution v n at step n and we will construct a solution v n+1 at step n + 1. First, we remark that if (v n+1,i ) i is a family indexed by i of discrete sub-solution at step n + 1, then v n+1 = max v n+1,i is still a sub-solution (it suffices to use Lemma 5.2). Moreover,

w + = sup Q ∆ |v n | and w -= -sup Q ∆ |v n |
are respectively discrete super and sub-solution of the scheme (2.9). Then, let us define

v n+1 = max{w subsolution at step n + 1 s.t w ≤ w + }.
Then v n+1 is a discrete sub-solution at step n + 1. Now let us prove that v n+1 is a super-solution. By contradiction, assume that there is I such that

v n+1 I -v n I ∆t < c[v] n+1 I G s[v] n I (D + v n+1 I , D -v n+1 I ).
This implies in particular that v n+1 I < w + . Now, let us consider, the solution w I (in a sense similar to Definition 5.1) of

w I -v n I ∆t = c ∆ v n+1 [w] I G + (D + v n+1 w I , D - v n+1 w I ) if c ∆ v n+1 [w] I ≥ 0 G -(D + v n+1 w I , D - v n+1 w I ) if c ∆ v n+1 [w] I < 0 , where c ∆ v n+1 [w] I = J∈Z N c0 I-J 1 {v n+1 J ≥wI } ∆x 1 ...∆x N - 1 2 J∈Z N c0 J ∆x 1 ...∆x N and D ± v n+1 w I = (D ± x1,v n+1 w I , ..., D ± xN ,v n+1 w I ) with D + xi,v n+1 w I = v n+1 I i,+ -w I ∆x i , D - xi,v n+1 w I = w I -v n+1 I i,- ∆x i ,
where I k,± is defined in (5.32). The existence of such a solution comes from the fact that the left hand side is non-decreasing in w I and the right hand side is non-increasing. Then, it is easy to prove (using Lemma 5.2) that w I > v n+1 I and that w defined by

w J = w I if J = I v n+1
J otherwise is a discrete sub-solution of (2.9) at step n + 1. This contradicts the definition of v n+1 and ends the proof of the proposition.

Proposition 5.4 (Properties of the discrete solutions)

Assume that ũ0 ∈ W 1,∞ (R N ). Let us denote by L = |Dũ 0 | L ∞ (R N ) and (5.38) v = sup{v subsolution of (2.9) with initial condition v 0 = ũ0 , v ≤ |ũ 0 | L ∞ (R N ) }.
For k ∈ Z N such that k • ∆X > 0 (with ∆X = (∆x 1 , ..., ∆x N )), the following estimates hold

(5.39) v(x I , t n ) -v(x I + k • ∆X, t n ) k • ∆X ≤ L, (5.40) v(x I , t n ) -v(x I + k • ∆X, t n ) k • ∆X ≥ -L, (5.41) v n+1 I -v n I ∆t ≤ √ N 2 |c 0 | L 1 (R N ) L.
The same estimates hold for

(5.42) v = inf{v supersolution of (2.9) with initial condition v 0 = ũ0 , v ≥ -|u 0 | L ∞ (R N ) }. Proof of Proposition 5.4 For k ∈ Z N , let us denote by ũ0,k (x I ) = ũ0 (x I + k • ∆X) + Lk • ∆X ≥ ũ0 (x I )
and by v k the greater sub-solution, with initial condition ũ0,k , i.e.

v k = sup{v k subsolution of (2.9) with initial condition v 0 k = ũ0,k , v k ≤ |u 0 | L ∞ (R N ) }.
First, since the equation does not see the constants, we remark that v -Lk • ∆X is a sub-solution of (2.9) with initial condition v 0 -Lk • ∆X ≤ ũ0,k . So, by definition of v k , we have

v -Lk • ∆X ≤ v k .
Moreover, using the fact that the scheme is invariant by translation in space, we deduce that v(x I +k•∆X, t n ) is the greater sub-solution with initial condition v(x I + k • ∆X, 0) = ũ0,k (x I ). We then deduce that

v k (x I , t n ) = v(x I + k • ∆X, t n ). This implies that for k • ∆X > 0 v(x I , t n ) -v(x I + k • ∆X, t n ) k • ∆X ≤ L.
The proof of (5.40) is similar and we skip it.

To prove (5.41), we use the fact that v is a solution and the two previous estimates. This implies

- √ N 2 |c 0 | L 1 (R N ) L ≤ v n+1 I -v n I ∆t ≤ √ N 2 |c 0 | L 1 (R N ) L
which is the desired result.

Before to prove Theorem 2.2, we need the following lemma:

Lemma 5.5 (Equivalent formulation for the discrete velocity)

The discrete velocity c ∆ [v] can be rewritten as

c ∆ [v] n I = c 0 ⋆ 1 {v # (.,tn)≥v # (xI ,tn)} (x I ) - 1 2 |c 0 | L 1 (R N )
where v # is defined in (2.13).

Proof of Lemma 5.5

The idea of the proof is borrowed from Alvarez et al. [START_REF]A convergent scheme for a nonlocal hamilton-jacobi equation, modeling dislocation dynamics[END_REF]. Using the definition of the discrete velocity and c0 , we get

c ∆ [v] n I = J∈Z N c0 I-J 1 {v n J ≥v n I } ∆x 1 ...∆x N - 1 2 J∈Z N c0 J ∆x 1 ...∆x N = J∈Z N QJ-I c 0 (y) dy 1 {v n J ≥v n I } - 1 2 J∈Z N QJ c 0 (y) dy = J∈Z N QJ c 0 (x I -y) 1 {v n J ≥v n I } dy - 1 2 |c 0 | L 1 (R N ) = J∈Z N QJ c 0 (x I -y) 1 {v # (y,tn)≥v # (xI ,tn)} dy - 1 2 |c 0 | L 1 (R N ) = c 0 ⋆ 1 {v # (•,tn)≥v # (xI ,tn)} (x I ) - 1 2 |c 0 | L 1 (R N ) .
This ends the proof of the lemma.

We now prove Theorem 2.2: Proof of Theorem 2.2

The proof is an adaptation of the one of Crandall Lions [START_REF] Crandall | Two approximations of solutions of Hamilton-Jacobi equations[END_REF], revisited by Alvarez et al. [START_REF] Alvarez | Convergence of a first order scheme for a non local eikonal equation[END_REF]. Nevertheless, for the reader's convenience, we give the main steps in order to show the new difficulties due to the non-local term. The main idea of the proof is the same as the one of comparison principles, i.e. to consider the maximum of u -v # , to duplicate the variable and to use the viscosity inequalities to get the result.

The proof splits into four steps. We first assume that (5.43) ū0 (x I ) ≥ (ũ 0 (x I )) # , for all I ∈ Z N . and we set (5.44)

µ 0 = sup R N (ū 0 -(ũ 0 ) # ) ≥ 0.
We denote throughout by K various constant depending only on N ,

|c 0 | W 1,1 (R N ) , |Dū 0 | L ∞ (R N ) , |Dũ 0 | L ∞ (R N )
and µ 0 .

Step 1: Estimate on v We have the following estimate for the discrete solution

(5.45) -Kt n ≤ ū0 (x) -v # (x, t n ) ≤ Kt n + µ 0 .
To show this, it suffices to use the estimate (5.41) of Proposition 5.4. This implies that, for x ∈ Q I (with the notation v and v defined in Proposition 5.4)

(ũ 0 ) # (x) -Kt n = v 0 I -Kt n ≤ v n I ≤ v n I = v # (x, t n ) ≤ v n I ≤ Kt n + v 0 I = Kt n + (ũ 0 ) # (x).
This implies the desired estimate.

Before continuing the proof, we need a few notation. We put

µ = sup R N (u -v # ).
We want to bound from above µ by µ 0 plus a constant. For every 0 < α ≤ 1, 0 < γ ≤ 1 and 0 < η ≤ 1, we set M α,γ η = sup

R N ×R N ×(0,T )×{0,...,tN T } Ψ α,γ η (x, y, t, t n ) with Ψ α,γ η (x, y, t, t n ) = u(x, t) -v # (y, t n ) - |x -y| 2 2γ - |t -t n | 2 2γ -ηt -α(|x| 2 + |y| 2 ).
We shall drop the super and subscripts on Ψ when no ambiguity arises as concerning the value of the parameter.

The main difference with the classical Crandall Lions proof, is to consider the function v # in the place of v in the separation of variables. This allows us to treat the non-local velocity.

Since u is Lipschitz continuous and T ≤ 1, we have (5.46) |u(x, t)| ≤ K(1 + |x|).

Moreover by

Step 1, we have

|v(x i , t n )| ≤|v(x i , t n ) -ū0 (x i )| + |ū 0 (x i )| ≤Kt n + K(1 + |x i |) ≤K(1 + |x i |).
We then deduce that Ψ achieves its maximum at some point that we denote by (x * , y * , t * , t * n ).

Step 2: Estimates for the maximum point of Ψ

The maximum point of Ψ enjoys the following estimates (5.47)

α|x * | + α|y * | ≤ K and (5.48) |x * -y * | ≤ Kγ, |t * -t * n | ≤ (K + 2η)γ.
Indeed, by inequality Ψ(x * , y * , t * , t * n ) ≥ Ψ(0, 0, 0, 0) ≥ 0, we obtain

α|x * | 2 + α|y * | 2 ≤u(x * , t * ) -v # (y * , t * n ) ≤ K(1 + |x * | + |y * |) ≤K + K 2 α + α 2 |x * | 2 + α 2 |y * | 2 .
This implies (5.47), since α ≤ 1.

The first bound of (5.48) follows from the Lipschitz regularity in space of u, from the inequality Ψ(x * , y * , t * , t * n ) ≥ Ψ(y * , y * , t * , t * n ) and from (5.47). The second bound of (5.48) is obtained in the same way, using the inequality Ψ(x * , y * , t * , t * n ) ≥ Ψ(x * , y * , t * n , t * n ).

Step 3: A better estimate for the maximum point of Ψ Inequality (5.47) can be strengthened to (5.49)

α|x * | 2 + α|y * | 2 ≤ K.
Indeed, using the Lipschitz regularity of u, the inequality Ψ(x * , y * , t * , t * n ) ≥ Ψ(0, 0, 0, 0), Step 1 and equation (5.48), yields

α|x * | 2 + α|x * i | 2 ≤u(x * , t * ) -v # (y * , t * n ) + u 0 (y * ) -u 0 (y * ) ≤K(|x * -y * | + t * ) + Kt * n + µ 0 ≤ K.
Step 4 : Upper bound of µ We have the bound µ ≤ K √ T (∆x + ∆t)

1 2 + µ 0 if ∆x + ∆t ≤ 1 K 2 .
First, we claim that for η large enough, we have either t * = 0 or t * n = 0. We argue by contradiction. Then the function (x, t) → Ψ(x, y * , t, t * n ) achieves its maximum at a point of R N × (0, T ). Then, using the fact that u is a sub-solution of (2.8), we deduce that (5.50)

η + p * t ≤ c[u](x * , t * )|p * x + 2αx * | where c[u](x, t) = c 0 ⋆ 1 {u(•,t)≥u(x,t)} (x) - 1 2 R N c 0 and p * t = t * -t * n γ , p * x = x * -y * γ .
Since t * n > 0, we also have Ψ(x * , y * , t * , t * n ) ≥ Ψ(x * , y, t * , t n ) for t n ≥ t * n -∆t. This implies

(5.51) v # (y * , t * n ) -v # (y, t n ) ≤ ϕ(y * , t * n ) -ϕ(y, t n ) for t n ≥ t * n -∆t for ϕ(y, t n ) = -|x * -y| 2 2γ + |t * -tn| 2 2γ
+ α|y| 2 . We denote by x * I the node such that y * ∈ Q I * (the unit cube centred in x * I ). We then deduce ϕ(y

* , t * n ) -ϕ(y * , t * n -∆t) ∆t ≥ v # (y * , t * n ) -v # (y * , t * n -∆t) ∆t = v(x * I , t * n ) -v(x * I , t * n -∆t) ∆t ≥c ∆ [v](x * I , t * n )G s[v](x * I ,t * n ) (D + v(x * I , t * n ), D -v(x * I , t * n )) ≥c ∆ [v](x * I , t * n )G s[v](x * I ,t * n ) (D + ϕ(y * , t * n ), D -ϕ(y * , t * n ))
where we have use the monotonicity of G ± and (5.51) with t n = t * n for the last line. We denote by s

* = s[v](x * I , t * n ). Straightforward computations of the discrete derivative of ϕ yield (5.52) p * t + ∆t 2γ ≥ c∆ [v](x * I , t * n )G s * p * x - ∆x 2γ -α(2y * + ∆x), p * x + ∆x 2γ -α(2y * -∆x) .
Subtracting the above inequality to (5.50) yields

η ≤ ∆t 2γ + c[u](x * , t * )|p * x + 2αx * | -c∆ [v](x * I , t * n )G s * p * x - ∆x 2γ -α(2y * + ∆x), p * x + ∆x 2γ -α(2y * -∆x) ≤ ∆t 2γ + (c[u](x * , t * ) -c∆ [v](x * I , t * n ))|p * x | + αK|x * | + c∆ [v](x * I , t * n ) G s * p * x - ∆x 2γ -α(2y * + ∆x), p * x + ∆x 2γ -α(2y * -∆x) -G s * (p * x , p * x ) ≤ ∆t 2γ + (c[u](x * , t * ) -c∆ [v](x * I , t * n ))|p * x | + αK|x * | + K ∆x γ + 2α|y * | + 2α∆x
We now have to bound the term

I =(c[u](x * , t * ) -c∆ [v](x * I , t * n ))|p * x | =(c 0 ⋆ 1 {u(•,t * )≥u(x * ,t * )} (x * ) -c 0 ⋆ 1 {v # (•,t * n )>v # (y * ,t * n )} (x * I ))|p * x |.
We have to distinguish two cases: Case 1. Assume that for γ > 0 fixed |x * -y * | 2 γ ≥ C 2 γ for some constant C γ > 0 and for all α small. By inequality Ψ(x * , y * , t * , t * n ) ≥ Ψ(x, x, t * , t * n ), we deduce that

u(x, t * ) -u(x * , t * ) ≤ v # (x, t * n ) -v # (y * , t * n ) - |x * -y * | 2 2γ -α(|x * | 2 + |y * | 2 ) + 2α|x| 2 .
We then get the inclusion

{u(•, t * ) ≥ u(x * , t * )} ∩ {v # (•, t * n ) ≤ v # (y * , t * n )} ⊂ {|x| 2 ≥ R 2 α,γ } where R 2 α,γ = 1 2α |x * -y * | 2 2γ + α(|x * | 2 + |y * | 2 ) . We also have {|x -x * | ≥ R α,γ } ⊂ {|x| ≥ Rα,γ } where Rα,γ = R α,γ -|x * | → ∞ as α → 0 (see Da Lio et al. [13, Lemma 2.5]
). We then obtain

I ≤(c 0 ⋆ 1 {v # (•,t * n )>v # (y * ,t * n )} (x * ) -c 0 ⋆ 1 {v # (•,t * n )>v # (y * ,t * n )} (x * I ) + c 0 ⋆ 1 B c (0,Rα,γ ) (x * ))|p * x | ≤K |Dc 0 | L 1 (R N ) |x * -x * I | + B c (0, Rα,γ ) c 0 (x)dx <K(γ + ∆x) + o α (1)
where we have used (5.48) to bound |p * x | in the second line. Case 2. Assume that there exists a subsequence α n > 0 which we still denote by α such that |x * -y * | 2 γ → 0 as α → 0.

We then get that p * x → 0 as α → 0. Since the velocities are bounded, we get that I = o α (1). To sum up, we have shown that I < K(γ + ∆x) + o α (1).

We then get that

η < K ∆x + ∆t γ + γ + o α (1). Putting η * = K ∆x + ∆t γ + γ + o α (1)
, we then conclude that either t * = 0 or t * n = 0 for η ≥ η * .

Whenever t * = 0, we get for η

* ≤ η ≤ 1 M =Ψ(x * , y * , 0, t * n ) ≤ u 0 (x * ) -v # (y * , t * n ) ≤K|x * -y * | + Kt * n + µ 0 ≤Kγ + µ 0
The same result holds whenever t * n = 0. We then deduce that for η = η * and (x,

t n ) ∈ R N × {0, ..., t NT } u(x, t n ) -v # (x, t n ) -η * T -2α|x I | 2 ≤ M ≤ Kγ + µ 0 .
Sending α → 0, taking the supremum over (x, t n ) and choosing γ = T 1/2 (∆x + ∆t) 1/2 , we get

µ ≤ K(∆x + ∆t) 1/2 √ T + µ 0 provided ∆x + ∆t ≤ 1 K 2 .
Using the same arguments of Alvarez et al.[1, Theorem 2], we easily deduce the result in the general case. This ends the proof of Theorem 2.2.

We now point out that there is no consistency error for the scheme (2.15). If there is a consistency error, then there will be a consistency error for the velocities, i.e. for all ϕ ∈

C 1 (R N × (0, T )) with ϕ C 1 (R N ×(0,T )) ≤ C then (5.53) c ∆ [ϕ |Q ∆ T ] n I -c[ϕ](x I , t n ) ≤ f (∆x, ∆t)
with f (∆x, ∆t) → 0 as (∆x, ∆t) → 0.

Proposition 5.6 (No consistency error)

There is no consistency error for the scheme (2.15), i.e. equation (5.53) does not hold.

Proof of Proposition 5.6

We have to prove that there exists a constant C 0 > 0 such that for all ∆x, ∆t > 0 there exists a function ϕ such that c ∆ [ϕ |Q ∆ T ] n I -c[ϕ](x I , t n ) ≥ C 0 and |ϕ| C 1 (R N ×(0,T )) ≤ C. To see this, it suffices to take a function ϕ which oscillates as shown in Figure 1. We now point out that the two scheme (2.9) and (2.15) are not monotone:

Proposition 5.8 (Non monotonicity of the scheme)

The two scheme (2.9) and (2.15) are not monotone for general c0 .

Proof of Proposition 5.8

1. The implicit scheme (2.9)

We give a simple counter-example in the one dimensional case. We take ũn i = v n i = i∆x. The goal is to construct a sub-solution u n+1 and a super-solution v n+1 such that

u n+1 i = i∆x + c 1 v n+1 i = i∆x + c 2 .
A simple computation gives We then deduce that c 1 -c 2 = ∆tc 0 0 > 0 (if c0 0 > 0). This implies that we can construct two different solutions (the supremum of sub-solutions and the infimum of super-solutions). This implies that the scheme (2.9) is not monotone.

The explicit scheme (2.15)

We also give a counter-example in the one dimensional case. We take

u n i =    -δ if i = 0 1 -δ if i = -1 0 otherwise and v n i =    δ if i = 0 1 -δ if i = -1 0 otherwise
We then have u n i ≤ v n i for all i ∈ Z. We then prove that u n+1 0 > v n+1 0 for δ small enough. First we remark that c ∆ [v n ] 0 ≤ 0 (if i =0,-1 c0 j ≥ c0 0 + c0 -1 ), so v n+1 0 ≤ δ. Moreover,

u n+1 0 = -δ + ∆t ∆x   j c0 j - 1 2 j c0 j   1 + δ 2 ≥ -δ + ∆t ∆x 1 2 |c 0 | L 1 (R) >δ if δ < ∆t 2∆x 1 2 |c 0 | L 1 (R) .
This contradicts the monotonicity of the scheme.

This ends the proof of the Proposition.

Numerical Simulations

In this section, we provide some numerical simulations. In a first subsection, we explain a method to solve numerically the implicit scheme (2.9)-(2.10)-(2.11). Then we provide a simple simulation concerning a collapsing circle, with the implicit scheme and an another one, to highlight the fattening phenomena, with the explicit scheme.

How to solve the implicit scheme?

To solve (2.9)-(2.10)-(2.11), we will use an iterative process. Assume that we have a solution v n at step n.

We now want to compute a solution v n+1 at step n + 1. To do this, for w ∈ R Q ∆ , we denote by w = Φ( w) the solution of the following auxiliary scheme: where I k,± is defined in (5.32). Since the left hand side of (6.54) is non-decreasing in w I and the right hand side is non-increasing, there exists a unique solution (in a sense similar of Definition 5.1) which can be computed using a dichotomy process. The important point is that if w is a sub-solution (resp. super-solution) then w = (w I ) I is still a sub-solution (resp. super-solution) and satisfies w I ≥ wI (resp. w I ≤ wI ) (see the proof of Proposition 5.3).

The idea is then to define w 0 = v n -C where the constant C is such that w 0 is a sub-solution of (2.9) and then to construct iteratively w k+1 = Φ(w k ). Setting v n+1 = lim k→∞ w k , we have that v n+1 is a solution of (2.9).

(4. 26 )h 3 α

 263 Dϕ(x 0 ) = Du 0,α (x 0 , t 0 ) = Dψ(x 0 , t 0 ) and (4.27)D 2 ϕ(x 0 ) ≤ D 2 ψ(x 0 , t 0 ).Moreover, with h = |x 0 -x 0 |, we have (by(4.22))and the set {ϕ(x) ≥ ϕ(x 0 )} is a ball of radius R = h. Let us define the error

Figure 1 :Remark 5 . 7

 157 Figure 1: Graph of the function ϕ defined in R.

( 6 . 2

 62 54)w I -v n I ∆t = c ∆ w [w] I G + (D + w w I , D - w w I ) if c ∆ w [w] I ≥ 0 G -(D + w w I , D - w w I ) if c ∆ w [w] I < 0 where c ∆ w [w] I = J∈Z N c0 I-J 1 { wJ ≥wI } ∆x 1 ...∆x N -1 J∈Z N c0 J ∆x 1 ...∆x N and D ± w w I = (D ± x1, ww I , ..., D ± xN , w w I ) with D + xi, ww I = wI i,+ -w I ∆x i , D - xi, w w I = w I -wI i,- ∆x i ,

Figure 3 :

 3 Figure 3: Evolution of the 8 at time 0, 0.05, 0.1 and 0.15.

  t). x| 4 ≤ u 0 (x, t) -u 0 (x, t) ≤ K|x -x|

	Therefore
	(4.21) |x -and so 1 4α
	(4.22)	|x -x| ≤ Kα 1/3 .
	Moreover, by (4.20), we get
	(4.23)	0 ≤ u 0,α (x, t) -u 0 (x, t) ≤ u 0 (x, t) -u 0 (x, t) ≤ K|x -x|
	Using (4.22) and (4.23), we then deduce for x ∈ R N and t ∈ (0, T )
	(4.24)

Acknowledgements

The author would like to thanks R. Monneau for proposing him the problem and for fruitful discussions in the preparation of this paper and E. Carlini for enlighting discussions about numerical simulations. The author was supported by the contract JC 1025 called "ACI jeunes chercheuses et jeunes chercheurs" of the French Ministry of Research (2003Research ( -2006)).

A collapsing circle

In this subsection, we provide a simple test with the implicit scheme concerning the evolution of a circle. The goal of this simple simulation is just to check that the circle will disappear with the good time. We take a circle of radius 1. The parameters are ∆x = 0, 05, ∆t = 0, 01 and ε = 0, 3. Moreover, we take the kernel

The initial condition is the distance to the circle. The result is shown in Figure 2.

Fattening phenomena

The second test is concerning with the evolution of the 8 to point out the fattening phenomena. This test have been made with the explicit scheme. We take two circles of radius 0, 58 such that they are tangent in one point and we look at the evolution of the level set 0 and -0, 06. The parameters are ∆x 1 = ∆x 2 = 0, 01, ∆t = 0, 0001 and ε = 0, 1. Moreover, we take the kernel

The results are provided in Figure 3.


Daniel F Bernardes

Matthieu Latapy

Fabien Tarissan

Relevance of SIR Model for Real-world Spreading Phenomena: Experiments on a Large-scale P2P System

Understanding the spread of information on complex networks is a key issue from a theoretical and applied perspective. Despite the effort in developing theoretical models for this phenomenon, gauging them with large-scale real-world data remains an important challenge due to the scarcity of open, extensive and detailed data. In this paper, we explain how traces of peer-to-peer file sharing may be used to this goal. We also perform simulations to assess the relevance of the standard SIR model to mimic key properties of spreading cascade. We examine the impact of the network topology on observed properties and finally turn to the evaluation of two heterogeneous versions of the SIR model. We conclude that all the models tested failed to reproduce key properties of such cascades: typically real spreading cascades are relatively "elongated" compared to simulated ones. We have also observed some interesting similarities common to all SIR models tested.

I. INTRODUCTION

Diffusion phenomena in complex networks -such as the spread of virus on contact networks, gossip on social networks and files in peer-to-peer (P2P) networks -have spawned an increasing interest in recent years. The boost of computer networks and online social network platforms offers data and new insights on these phenomena in large scale networks; the possibility to validate and refine current models might lead to breakthroughs in the field.

Although large scale diffusion phenomena have always attracted considerable interest, it has been historically

challenging to obtain open, extensive and detailed real-world data at this level. Despite this obstacle, diffusion models emerged, notably in epidemiology. The early models, both discrete and continuous (see [START_REF] Andersson | Stochastic Epidemic Models and Their Statistical Analysis[END_REF], [START_REF] Anderson | Infectious Diseases of Humans: Dynamics and Control[END_REF] for a survey), focused primarily on macroscopic aspects of diffusion -such as the evolution of the number of infected individuals in a population -overlooking the microscopic dynamic of the epidemic -i.e., how (by whom) individuals become infected. The advent of network analysis in various contexts has pushed for a more detailed description of the diffusion process. Indeed, models based on the detailed interactions of agents on a network have blossomed in sociology [START_REF] Granovetter | Threshold Models of Collective Behavior[END_REF], computer science [START_REF] Easley | Networks, Crowds, and Markets -Reasoning About a Highly Connected World[END_REF] and economics [START_REF] Jackson | Social and Economic Networks[END_REF], among others. New epidemic models inspired by the classical approaches featuring a detailed dynamic description in the context of networks also appeared (see [START_REF] Barrat | Dynamical Processes on Complex Networks[END_REF], [START_REF] Draief | Epidemics and rumours in complex networks[END_REF] for a survey). In particular the network version of the SIR family of models has established itself as reference model in the study of information diffusion [START_REF] Cointet | How realistic should knowledge diffusion models be?[END_REF], [START_REF] Newman | The structure and function of complex networks[END_REF], [START_REF] Colizza | The role of the airline transportation network in the prediction and predictability of global epidemics[END_REF], [START_REF] Leskovec | Cascading Behavior in Large Blog Graphs[END_REF], [START_REF] Onnela | Structure and tie strengths in mobile communication networks[END_REF].

In this context, assessing the pertinence of such models to describe real-world data is critical. In order to validate this model a comprehensive empirical spreading trace, consisting of (1) detailed chronological data of who transmitted the information to whom and (2) data describing the underlying network on which the diffusion process takes place. Indeed, the network version of the SIR model (henceforth called simply SIR model) is based on local rules of transmission which take into account the network topology. In large epidemic bursts the available data often provides the evolution of large aggregate quantity, such as the number of touched individuals, but rarely uncover the local trail of the epidemic. Conversely, other empirical studies feature transmission events, but lack complete information of the underlying network structure on which the diffusion takes place [START_REF] Adar | Implicit structure and the dynamics of blogspace[END_REF], [START_REF] Iribarren | Impact of Human Activity Patterns on the Dynamics of Information Diffusion[END_REF]. This work analyses the relevance of the SIR model for real-world diffusions, using data obtained measuring the activity on a peer-to-peer file sharing network. This rich dataset allows one to reconstruct both the underlying network and the detailed diffusion trail at a remarkable scale.

We begin with a description of our dataset and framework in section II. In section III we define the spreading cascade. Next, in section IV, we simulate the spreading of files as a standard SIR process and confront it with the observed spreading; we also investigate the interplay between this process and structural properties of the underlying network where the spreading takes place. In section V we examine the spreading pattern when we modify the SIR model to account for heterogeneity in the behavior of the peers and in the popularity of files. We conclude the paper with of future work perspectives.

II. DATASET AND FRAMEWORK

The data used in this study comes from file sharing in an eDonkey server, obtained from a measurement of six hours of activity (akin to [START_REF] Aidouni | Ten weeks in the life of an edonkey server[END_REF]). In this setting, peers query the eDonkey server indexing files and for each file they get a list of available peers in the network possessing the requested file. Next, peers contact potential providers directly and transmission between them ensues. This dataset is a collection of answers to these queries, encoded as 4-tuples of integers in the following format: (t, P, C, F ), where capital letters represent unique ids (e.g. in Fig. 2). Each tuple accounts for a query made at time t of the file F by the peer C, satisfied by the peer P -that is, P has provided F to the peer C at time t. Let D be the set of all recorded tuples, P the set of all peers in these tuples and F the set of all files exchanged. In our dataset we have |P| = 1 908 500 peers, |F| = 801 280 files and |D| = 22 944 800 file transfers.

A. Underlying network

The trace D naturally induces a relationship between files and peers (who request or provide them), which we encode in a bipartite graph B = (P, F, A) on the disjoint sets of peers P and F files respectively. Let (t, P, X, F ) ∈ D be a recorded transmission of the file F by the peer P to some peer X at some time t, which we denote simply by (•, P, •, F ). Likewise, let (•, •, P, F ) ∈ D be a recorded transmission of the file F to the peer P , provided by some peer at some time instant. Hence:

A = {(P, F ) ∈ P × F : (•, P, •, F ) ∈ D ∨ (•, •, P, F ) ∈ D}
To study the diffusion, it is necessary to define the underlying network on which spreading takes place. Focusing on information content diffusion among peers, it is natural to consider the interest graph in which each node represents a peer and each edge joining two peers stand for common interest. Interests connecting peers may include broad subjects such as open source software, folk rock or French literature or narrow ones such as movies by Quentin Tarantino, a particular computer game or pictures of Beijing. It is reasonable to suppose that peers store and share content related to their interests and, likewise, peers will search for content matching their interests. Hence the diffusion of files among peers takes place on the interest graph and occurs from neighbor to neighbor. Indeed, if a peer P provides a file F (corresponding to a music album for example) to another peer P then there is link between them in the interest graph, since both are interested in the same content, namely F .

It is beyond doubt extremely difficult in a large scale interaction network to know precisely whether any two individuals have a common interest. Nonetheless, it is possible to approximate this graph using the data in D: the inferred interest graph is given by the projection G = (P, E) of B on P, connecting the peers who belong to the neighborhood of a common file in the bipartite graph, for each file: E = {(P, P ) ∈ P ×P : ∃F ∈ F, (P, F ) ∈ A ∧ (P , F ) ∈ A} See example in Fig. 1. For the sake of readability the inferred interest graph will be henceforth called simply interest graph.

B. Observed network structure

We begin examining properties of the bipartite graph B constructed from the P2P diffusion trace. In order to estimate the typical number of interested peers per file we have calculated the median degree of the files in the bipartite graph, 5, and the average degree, 14.73, with standard deviation of 34.74. Likewise, we have calculated the same statistics for the peers, to estimate the number of files commonly shared by peers: its median degree in the bipartite graph is 3 and the average degree is 6.19, with corresponding standard deviation of 12.66. The degree distribution of both peers and files is however heterogeneous (Fig. 3b) and mostly concentrated on small values; all degree values for peers and files remain below 10 4 .

The interest graph obtained from the observed bipartite graph (as explained above and in Fig. 1) has a single giant component containing almost all nodes (99.99%) and density 2.62 × 10 -4 . In Fig. 3a we have plotted the degree distribution for the peers: considering the set of all peers, the median degree is 118 and the mean value is 500.11, with corresponding standard deviation of 1271.42. We proceed to a finer analysis of the degree distribution, grouping peers in categories (Fig. 3a). Let us consider first the set of clients C ∈ P such that (•, •, C, •) ∈ D: i.e., peers having requested files during our measurements. Their degree distribution superposes the degree distribution of all nodes. This is due to the fact that 99.63% of peers in our observations have Fraction of cascades with a num. of initial providers > x Initial providers (c) Complementary cumulative distribution of the number of initial providers in the spreading cascades Fig. 3. Properties of the underlying network and observed spreading cascades requested at least one file, so the clients degree distribution is essentially the global degree distribution. A much more restrictive category is the set of providers P such that (•, P, •, •) ∈ D, i.e., peers having supplied files during our measurements. They account for 4.33% of the peers in P. Their degree distribution has a similar shape, but it is concentrated on larger values, indicated by a median of 1821 and an average degree of 2906.54 -with corresponding standard deviation of 3471.80. The last curve, superposing the curve corresponding to the providers, represents the degree distribution of a particular subset of providers called initial providers, which will be detailed in the next section.

We close this section with a brief summary of our dataset. Using the framework introduced, we were able to reconstruct the interest graph of the peers, where the spreading of files takes place. This graph connects essentially all peers, which can be grouped in two categories: providers and clients. Most peers in our observations are clients, but only a small fraction supply files: this revels a high proportion of free-riders (peers obtain files and do not share back) in the P2P-network observed. Furthermore, there is a sharp distinction between clients and providers in terms of their degree distribution.

III. SPREADING IN OUR DATA

In this work we analyze the spreading cascade representing the diffusion of each file in the P2P network. For a file F , the spreading cascade is a directed graph featuring the set P F of peers who have participated in the spread of F (as clients and/or providers) and links P → C, connecting each client C with the first peer(s) who provided F to it. More formally, let τ F (C) = inf{t : (t, •, C, F ) ∈ D} be the first instant C obtained F and let the directed graph K F = (P F , L F ) be the spreading cascade of F , with

P F = {P ∈ P : (P, F ) ∈ A} L F = ∪ C∈P F {(P, C) ∈ P F × P F : (τ F (C), P, C, F ) ∈ D}
A client requesting a file may receive a response from potentially several providers simultaneously, which implies that nodes in the cascade graph not only have multiple outgoing links, but also multiple incoming links in general. The causality induced by the fact that we only consider the links corresponding to the first time a node received F prevents the appearance of cycles. Hence the cascade is in fact a directed acyclic graph (DAG).

The first key property encoded in the spreading cascade of a given file F is the number of nodes who possess it at the end of the observed period, which is given by the size of the cascade |P F |. We also explore two other key topological properties of the cascade, namely its depth and number of links. The former is defined as the length of the longest path on the cascade and captures the maximum number of hops from peer to peer that the file has undergo before it was relayed from a provider to a client. The number of links, given by |L F |, combined with the size of the cascade gives information on the sharing pattern of the network. An example of observed trace and constructed spreading cascade is given in Fig. 2: the spreading cascade has size 7, depth 3 and 6 links.

Another relevant spreading data concerns the initial providers for each file F , namely the set of peers that possessed it prior to any transfer activity on the observed trace. These nodes are the origin of the spreading cascade, triggering the diffusion of the file F . This information can also be inferred from the request log and be determined in the following way. Let C F (t) = {C ∈ P : (t , •, C, F ) ∈ D, t < t} be the set of peers who requested F prior to t. We define the set of initial providers of F as the set of peers P who have provided F at some time t, without having obtained it before t from another peer in the network:

I F = {P ∈ P : (t, P, •, F ) ∈ D, P / ∈ C F (t)}
Plotting the complementary cumulative distribution of the number of initial providers for the spreading cascades (Fig. 3c) we obtain an interesting curve, revealing a scalefree distribution. This means that although most spreading cascades in our observation have few initial providers, there is a non negligible fraction of cascades with a large number of initial providers.

IV. SIMPLE SIR MODEL

As mentioned in the introduction, we have decided to investigate the file spreading in the light of the simple SIR model. In our setting, each file spreading corresponds to an independent epidemic in the interest graph, in which each node is in one of the following states: susceptible, infected or non-interacting (sometimes denoted removed, hence the acronym SIR). Susceptible nodes do not possess the file and may receive it from an infected node, thus becoming infected. Infected nodes, in turn, spread the file to each of its neighbors, independently, with probability p and become promptly non-interacting thereafter. Although non-interacting nodes remain in this state, infected nodes may unsuccessfully try to infect them sending the file.

Supposing the observed diffusion trace was the result of such a simple SIR epidemic we may estimate the spreading parameter p. Each neighbor-to-neighbor transmission trial can be seen as a Bernoulli random variable, whose value is 1 in case of success and 0 otherwise and whose expected value is p. Assuming each trial independent and the parameter p is homogeneous for each P and F , we may estimate it by the empirical proportion of successes over all trials. Since each tuple in D accounts for a successful neighbor-to-neighbor transmission, |D| is the number of successful trials for all diffusion cascades. The total number of trials, in turn, is given by the sum of the degrees of all nodes involved in the spreading of each file. Hence, we obtain the following estimate, with a 95% confidence interval p ± 10 -6 :

p = |D| / F ∈F P ∈P F d(P ) = 1.063 × 10 -3
Since the simple SIR model depends upon a single parameter, namely the spreading probability p, we have fully characterized it with the preceding estimation.

A. The underlying network influence

The goal of simulating the standard SIR model and comparing the simulated cascades with the observed ones is primarily to assess how realistic this model would perform on the interest graph, in terms of size, depth and number of links of the spreading cascades. Secondly, we wish to compare the results with simulations on random networks to understand the role of the network topological structure on the shape of the spreading cascades generated with the SIR model. With this aim, we have considered the spreading of files in a sequence of random networks derived from the interest graph, with increasing topological complexity. More precisely we begin considering an Erdös-Rényi (ER) random graph with the same density of our interest graph, the simplest random graph in our sequence. Then we have chosen a random graph with the same density and degree distribution using the Configuration Model (CM) approach [START_REF] Newman | The structure and function of complex networks[END_REF]. Next we have generated a Random Bipartite (RB) graph, with the same density and degree distribution as our original bipartite graph B of peers and files [START_REF] Guillaume | Bipartite structure of all complex networks[END_REF]. Compared to the interest graph, the projection of this random bipartite graph has similar density, degree distribution and clustering coefficient. In sum, for each new element of this sequence of (uniformly chosen) random graphs we introduce a new constraint to make it more realistic -in the sense that its topological properties will be closer to the interest graph.

B. File spreading simulation

Combining the network topology, the initial condition information (the list of initial providers I F calculated for each file F ) and the calibrated spreading parameter p we can proceed to the simulations for each underlying network: for each F , we begin with the initial providers in an infected state and the other nodes in a susceptible state. At each step, infected nodes will infect each of its neighbors with probability p, becoming non-interacting afterwards. The epidemic continues as long as there are active infected nodes.

The first observation concerning the model simulation is that the observed time (measured in seconds) has no direct relation with the simulation time (number of steps). Furthermore, our dataset corresponds to an observation in a bounded window of time of six hours, so that we have no reason to suppose that the file spreading cascades we observe correspond to the whole spreading cascade of a file. In other words, if we had measured a longer time window we would likely observe bigger cascades (in terms of size and depth) for the same files -due to, among other reasons, new users who could eventually request the same files. This is also true for our SIR model: we observe increasingly bigger cascades as time increases. In fact performing unconstrained simulations we have obtained a distribution of significantly bigger cascades than the ones we have observed in the real trace. Thus, in order to perform a suitable comparison with the observed cascades, we have decided to hold one property fixed and compare the other properties. More precisely for each file we generate a simulated cascade with the same size (resp. depth) as the corresponding observed cascade and compare the depth (resp. size) and number of links. In practice, for each file we simulate the SIR epidemic as described earlier and halt it when it reaches the size (resp. depth) of the corresponding observed cascade.

We have generated populations of simulated cascades for each underlying network and constraint (on depth and size). We have performed 801 280 file spreading simulations (one for each file in F) for each network and have selected every simulated file spreading cascade which attained the depth (resp. size) of the real spreading cascade for the same file - and have rejected the others for purpose of comparison. With this procedure, each underlying network yields a different population of file spreading cascades, since the rejected cascades may be different in each case. However 93.80% of the files have generated simulated cascades with the same depth as the corresponding real cascades, for all networks. Similarly, 85.64% of the files have generated simulated cascades with the same size as the corresponding real cascades, for all networks -except the ER network. Indeed, only 21.76% of the files have generated the contemplated size in the ER graph. Furthermore the properties of these simulated cascades on the ER graph deviated significantly from the properties of the cascades on the other graphs. Hence, in the following analysis we do not include the simulations for the ER graph. Rather, we focus on the properties of the of files with comparable spreading cascade depth (resp. size) on all networks but ER.

In Fig. 4a we plotted the complementary cumulative distribution of the size of cascades with comparable depth.

We observe a divergence of the cascade size from the observed cascades: simulated cascades are typically much bigger in size for a given depth compared to real cascades. The range of values in both categories is also striking: the biggest real cascade is at least two orders of magnitude smaller than the biggest simulated ones. Among the simulated cascades, there is a remarkable matching in size values for the simulation on the CM and the interest graph (curves are superposed). In Fig. 4c we plot the complementary cumulative distribution of the depth of cascades with fixed size. Real cascades feature a much higher depth compared to simulations, holding cascade size constant. In particular there is a cutoff on the cascade depth for the simulations: we do not observe any cascade depth bigger than 11 in the simulations. As for the number of links, we have two interesting situations. If we fix the depth (Fig. 4b) the number of links distribution resembles closely the size distribution (Fig. 4a). This is not completely surprising, since the two quantities are related. In this case we observe a larger number of links for all simulations compared to the number of links in the real cascades since the simulated cascades itselves are bigger. If, in constrast, we fix the cascade size to fit the observed cascades size (Fig. 4d), we observe a typically smaller number of links. Combining these observations on both plots we conclude that real spreading cascades are denser than simulated ones, a clear qualitative feature not captured by the simple SIR model. Finally we note that most cascades are simple, featuring depth equal to one and correspondingly small size.

To sum up, we have compared simple topological properties of real spreading cascades and simulated cascades from a calibrated SIR model, with comparable depth and size. We have observed that simulated cascades are relatively "wider" whereas real cascades are relatively "elongated", that is, real cascades have a smaller size per depth ratio. Moreover, real cascades are typically denser than simulated ones. In terms of interplay between underlying network structure and the simple SIR spreading cascades, we have observed that respecting the interest graph degree distribution was the only property that caused a striking change in simulations behavior on the considered random networks. Indeed we have observed sharp qualitative dissimilarities between the simulations on the ER graph (different degree distribution) and no sensible dissimilarities between the simulations on the CM, RB and the interest graphs.

V. HETEROGENEOUS SIR MODELS

In the previous section we have examined the adequacy of the simple SIR model to generate verisimilar file spreading cascades. We have also inspected the interplay between the underlying network and the model simulating file spreading in different networks. In this section we perform a complementary analysis, focusing on a single underlying network and examining different extensions of the SIR model considered previously. In particular we consider two heterogeneous versions of the SIR model, characterized by a distribution of spreading probabilities, instead of a single homogeneous parameter. The natural choice in this case for the underlying network is the interest graph, which is the most complete and realistic graph among the ones tested in the previous section.

A. File popularity

A first refinement of the simple SIR model consists in introducing different spreading probabilities according to the file being spread. The rationale in this case is to account for different levels of popularity depending on the file. Exogenous reasons -such as a movie release or the death of an an artist -can change the supply and demand of a given file and consequently alter its spreading probability. The knowledge of the actual reasons that explain the heterogeneity in file popularity are irrelevant to the characterization of this model, if we know the spreading probabilities for each file, i.e., {p(F ) : F ∈ F }. An estimate of these probabilities, in turn, can be obtained from the trace D if we suppose it was generated by a process following this extended SIR model. Indeed, since each file spreading is independent of the others, it is possible to estimate p(F ) for each F separately, with the same method used to derive the homogeneous parameter. Restricting the calculations to the spreading cascade of F , p(F ) will be given by the empirical proportion of successful transmissions of F over all possible transmissions of F :

p(F ) = |{(•, •, •, F ) ∈ D}| / P ∈P F d(P )
In Fig. 5a we plot the distribution of the heterogeneous spreading parameters depending on the files. The values of p are concentrated on the range 10 -5 to 10 -2 , indicating that there is a considerable fraction of cascades with a significantly different spreading regime (bigger than one order of magnitude). This distribution characterizes the extended SIR model we use in the following simulations.

B. Peer behavior

A second possible refinement is motivated by the fact that peers might have intrinsically distinct levels of "generosity" regarding file sharing. Under this hypothesis we extend the standard SIR model assigning an heterogeneous spreading probability to each peer, regardless of which file it is sharing. Thus, we do not need any other information but the spreading probability distribution to characterize the model. In this context altruistic peers, who typically spread files to a large proportion of their neighbors, would feature a bigger spreading probability compared to the homogeneous spreading probability corresponding to the diffusion aggregates of all peers. By the same token, the extreme case of free-riders would have their spreading probability assigned to zero. Again we can study transmissions as outcomes of Bernoulli trials to estimate the spreading probabilities. Let F P = {F ∈ F : (P, F ) ∈ A} be the files carried by the peer P ; for each such file the number of transmission trials P could perform corresponds to its degree in the interest graph, namely d(P ). Hence, to obtain p(P ) for each peer P we divide the number of successful transmissions of P to other peers (of any file carried by P ) over the total number of potential trials: We have plotted the distribution of the positive spreading probabilities estimates in this case (Fig. 5b). They account for small fraction of all the peers, since the only peers who have a positive spreading probability are those who provided a file at least once -4.33% cf. observations made in section II. Conversely, a large fraction of the peers do not share the file in this model. We observe a marked range of values, which is significantly greater than the one calculated for the homogeneous SIR.

C. File spreading simulation

Our aim is to generate simulated cascades following both extensions of the SIR model presented -with heterogeneous spreading probability depending on the files and on the peers -and compare their properties with the simulated cascade of the simple SIR model and the real observed cascades. In this sense, we apply the same methodology of the previous simulations: we fix the depth (resp. size) for the simulated cascades and examine the other two properties -the idea is to compare similar spreading cascades in terms of the chosen property. As discussed previously, the great majority of the cascades is simple, with depth equal to one and a small size. Hence the simulated cascades corresponding to the simple observed cascades will likely correspond in terms of depth, size and number of links. For this reason, we have decided in this section to focus on the spreading cascades with depth greater than one.

The simulation results are plotted in Fig. 6: we have plotted the complementary cumulative distributions of the spreading cascade depth, size and number of links. Imposing a constrain on the depth for the simulated cascades and comparing their size (Fig. 6a) we observe the contrast between the simulated and the real observed cascades with the same depth: the former have a typically bigger size compared to latter. What is remarkable, however, is the agreement among all the simulated cascade distributionscurves superposed in Fig. 6a. Next, if we fix the size for the simulated cascades and examine their depth, we are faced with the same qualitative similarity among simulated curves.

Indeed, the curves corresponding to the heterogeneous SIR models also feature a cutoff in depth, failing to reproduce the scale-free curve representing the depth of the observed real cascades. Finally, the cascade links distribution plotted in Fig. 6b and Fig. 6d reveals the pattern observed previously, namely that the observed spreading cascades are typically denser than corresponding simulated cascades.

Inspite of the improvements in the SIR model, introducing an heterogeneous spreading parameter to account for different profile of files (respectively peers), the simulations indicate that this refinement does not change qualitatively the basic properties of the simulated spreading cascades. Indeed we observe a surprising agreement between the three SIR models compared, notwithstanding the particularities of each model.

VI. CONCLUSION AND PERSPECTIVES

We have presented a large-scale dataset from a real-world peer-to-peer network, featuring diffusion of files among peers. We have proposed a framework to study this dataset which allows us to obtain, simultaneously, the interest graph of peers -where the diffusion of content takes place -and the spreading cascade. Guided by simulations we have examined spreading cascades generated by the simple SIR model and have analyzed the interplay between this model and the network topology. We concluded that simulated file spreadings do not capture key qualitative properties of the observed spreading cascades. Furthermore, in terms of the studied properties, the simple SIR model generates similar cascades on random networks having the same degree distribution as the interest graph. Next we have focused on the spreading of files on the interest graph and studied extended versions of the SIR model featuring an heterogeneous spreading parameter. Surprisingly enough, simulated cascades using both extensions of the SIR model show similar properties as the simple homogeneous SIR model -and thus, fail to reproduce qualitative features of the observed cascades.

The SIR model is an attractive choice to model the information spreading in complex networks: it is based on classical epidemiological models, it is based upon few assumptions and can be characterized with one parameter. However, the results on this paper suggest that this model might not be suited to describe file spreading in our data. Furthermore, extensions of this epidemic model to make it more realistic, featuring heterogeneous spreading probabilities do not offer a better alternative in terms of the properties we observed. At this point, we consider two main exploration tracks. The first possibility consists in constructing a weighted interest graph, which takes into account the number of interactions (file exchanges) between peers. In this case the same analysis may be performed and a comparison with the results of this paper would be pertinent. The second possibility is to contrast epidemiological models to adoption/threshold models [START_REF] Cointet | How realistic should knowledge diffusion models be?[END_REF], [START_REF] Granovetter | Threshold Models of Collective Behavior[END_REF].

Fig. 1 .

 1 Fig. 1. Interest graph as a projection of the bipartite graph of peers and files contructed from the trace D.

Fig. 2 .

 2 Fig. 2. Trace log example with corresponding spreading cascade in black and underlying network in light gray.

  Degree distributions on the interest graph. Superposed curves: all peers and clients, providers and initial providers Fraction of nodes with degree > x Bipartite graph of peers and files Peers Files (b) Complementary cumulative degree distributions of peers and files on the bipartite graph B

  Size of cascades with fixed depth. Curves corresponding to the interest graph and CM superposed.

  Number of links of cascades with fixed depth. Curves corresponding to the interest graph and CM superposed.

  Fraction of cascades with number of links > x Number of links of cascades with fixed size Config. Model Random bip Interest graph Observed trace (d) Number of links of cascades with fixed size. Curves corresponding to the interest graph, RB and CM superposed.

Fig. 4 .

 4 Fig. 4. Simulation of file spreading on different underlying networks: complementary cumulative distribution of cascade properties

Fig. 5 .

 5 Fig. 5. Heterogeneous spreading parameter distributions

Fig. 6 .

 6 Fig. 6. Simulation of file spreading on the interest graph with different SIR processes: complementary cumulative distribution of cascade properties

ACKNOWLEDGMENT

This work is partly funded by the European Commission through the FP7-FIRE project EULER (Grant No.258307) and by the City of Paris Émergence program through the DiRe project.