Pascal Azerad
email: azerad@math.univ-montp2.fr

Afaf Bouharguane

Finite difference approximations for a fractional diffusion/anti-diffusion equation

Keywords: Anti-diffusive fractional operator, Finite difference approximations, Von Neumann stability analysis, Error analysis. Mathematics Subject Classification: 35L65, 45K05, 65M12, 65M06

A class of finite difference schemes for solving a fractional anti-diffusive equation, recently proposed by Andrew C. Fowler to describe the dynamics of dunes, is considered. Their linear stability is analyzed using the standard Von Neumann analysis: stability criteria are found and checked numerically. Moreover, we investigate the consistency and convergence of these schemes.

Introduction

Partial Differential Equations with nonlocal or fractional operators are widely used to model scientific problems in mechanics, physics, signal processing, see for example [START_REF] Azerad | Simultaneous denoising and enhancement of signals by a fractal conservation law[END_REF] and references therein. We consider in this chapter a nonlocal conservation law which appears in the formation and dynamics of sand structures such as dunes and ripples [START_REF] Fowler | Dunes and drumlins[END_REF][START_REF] Lagree | Stability of an erodible bed in various shear flows[END_REF]. Since it is generally impossible to obtain analytical solutions of these nonlocal models, one must rely on numerical solutions. In the last few decades, significant advances in numerical analysis and computational implementation of numerical methods for nonlocal/fractional PDEs have been made. For instance, [START_REF] Droniou | A numerical method for fractal conservation laws[END_REF] propose a finite volume method to approximate the solutions of a fractal scalar conservation law, that is to say a conservation law regularized by a diffusive fractional power of the Laplacian operator and [START_REF] Meerschaert | Finite difference approximations for fractional advectiondispersion flow equations[END_REF][START_REF] Yuste | Weighted average finite difference methods for fractional diffusion equations[END_REF] use finite difference methods to approximate fractional diffusive equations. In this chapter, we develop the basic numerical analysis of the following evolution equation proposed by A.C. Fowler (see [START_REF] Fowler | Dunes and drumlins[END_REF], [START_REF] Fowler | Evolution equations for dunes and drumlins[END_REF] and [START_REF] Notes H T T P ; W W W2.M A T H S.O X.A C.U K/˜f O W L E R/C O U R S E S/M A T H E N V O.H T M L | Mathematics and the environment[END_REF] for more details) to study the nonlinear dune formation:

∂ t u(t, x) + ∂ x u 2 2 (t, x) + η I[u(t, •)](x) -ǫ ∂ 2 xx u(t, x) = 0, (1)
where u = u(t, x) represents the dune height and I is a nonlocal operator defined as follows: for any Schwartz function ϕ ∈ S(R) and any x ∈ R,

I[ϕ](x) := +∞ 0 |ξ| -1 3 ϕ ′′ (x -ξ) dξ. (2)
The second and fourth terms of equation ( 1) correspond to the nonlinear and dissipative terms respectively, while the third term is the nonlocal term, which is anti-dissipative as we will show later on. The positive parameters ǫ (resp. η) quantify the amount of local diffusion (resp. nonlocal anti-diffusion).

Remark 1. For causal functions (i.e. ϕ(x) = 0 for x < 0), this operator is, up to a multiplicative constant, the Riemann-Liouville integral which is defined as follows:

1 Γ( 2 3 ) +∞ 0 ϕ ′′ (x -ξ) |ξ| 1/3 dξ = d -2/3 dx -2/3 ϕ ′′ (x) = d 4/3 dx 4/3 ϕ(x), (3)
with Γ the Euler function.

Many numerical methods for the evaluation of fractional order integrals and the solution of fractional order equations are proposed in the literature. Usually, time and spatial fractional derivatives are considered: we refer for instance to [START_REF] Diethelm | Algorithms for the fractional calculus: a selection of numerical methods[END_REF][START_REF]A fast time stepping method for evaluating fractional integrals[END_REF][START_REF] Yuste | Weighted average finite difference methods for fractional diffusion equations[END_REF]. In our case, the integral operator I can be seen as a fractional power of order 2/3 of the Laplacian with the bad sign. Indeed, it has been proved that I has the following Fourier transform [START_REF] Alibaud | A non-monotone nonlocal conservation law for dune morphodynamics[END_REF]:

F(I[ϕ])(ξ) = ψ I (ξ)Fϕ(ξ), (4)
where ψ I (ξ) = -a I |ξ| 4 3 + ib I ξ|ξ| 1 3 with a I = 2 π 2 Γ( 2 3 ), b I = 2 π 2 √ 3 Γ( 2 3 ) and F denotes the Fourier transform defined for f ∈ L 1 (R) by: for all ξ ∈ R Ff (ξ) = R e -2iπxξ f (x) dx. Formula (4) stems from the following integral formula [START_REF] Alibaud | A non-monotone nonlocal conservation law for dune morphodynamics[END_REF]:

I[ϕ](x) = 4 9 0 -∞ ϕ(x + z) -ϕ(x) -ϕ ′ (x)z |z| 7/3 dz. (5)
Finally, equation [START_REF] Alibaud | A non-monotone nonlocal conservation law for dune morphodynamics[END_REF] involves two antagonistic terms: the anti-diffusive operator I which creates instabilities and the diffusion operator -∂ 2 xx which controls these perturbations. Recently, some theoretical results regarding the Fowler model (1) have been obtained, namely, existence of travelling-waves, the global well-posedness, the failure of the maximum principle and the instability of constant solutions [START_REF] Alibaud | A non-monotone nonlocal conservation law for dune morphodynamics[END_REF]2,[START_REF] Bouharguane | On the instability of a nonlocal conservation law[END_REF]. The last two results are a consequence of the non-positivity of the kernel K of I -∂ 2 xx defined for t > 0 and x ∈ R by

K(t, •)(x) := F -1 (e -t(4π 2 |•| 2 +ψ I (•)) )(x). (6)
These two "bad properties" show that the discrete problem must be handled with care. Indeed, for monotone models, a classical way to get numerical stability criteria for explicit scheme is to ensure that the approximated problem satisfies the discrete maximum principle, which cannot be true for Equation [START_REF] Alibaud | A non-monotone nonlocal conservation law for dune morphodynamics[END_REF].

In [START_REF] Alibaud | A non-monotone nonlocal conservation law for dune morphodynamics[END_REF], some numerical results regarding this equation have been obtained using an explicit finite difference scheme but the detailed numerical study was not performed. Hence, in this chapter, we would like to go one step further investigating the numerical stability, consistency and convergence of a class of explicit finite difference schemes approximating the Fowler equation. The numerical stability is specially interesting here because the growth of the solution depends on frequencies and time. Hence, the notion of A-stability, also called strong stability, is not suitable nor desirable. In the literature, some authors use another definition of stability, less restrictive than the Astability: the C-stability. This is an abbreviation for convergence stability and is linked with stability in the Lax-Richtmyer sense. In this definition, a numerical scheme is stable for the norm || • || if for all T > 0, there exists a constant K(T ) > 0 independent of the time and space steps δx, δt such that for all initial data u 0

||u n || ≤ K(T )||u 0 ||, ∀ 0 ≤ n ≤ T δt ,
where u n represents the approximated solution at the time t n = nδt. This definition allows the solution to grow with time, which is the case for example for the equation u t -u xx = cu. For the L 2 -stability, a simple way to prove the numerical stability and specially to get stability criteria is Fourier analysis, see Section 3. Hence, considering the C-stability, the Von Neumann condition is written as

∃C > 0, ∃δt * > 0, such that ∀δt ∈]0, δt * ]; ∀k ∈ Z |g(k)| ≤ 1 + Cδt, (7)
where g is the discrete amplification factor, k the wave number and C is a positive constant independent of δx and δt. If C = 0, the Von Neumann condition coincides with the A-stability. As we will see later in Section 2, the amplification of solutions of the Fowler equation also depends on frequencies: low frequencies are slowly amplified whereas the high frequencies are dampened. Hence, the notion of C-stability is not adapted for this model because it considers only the amplification due to time. To take into account this phenomenon, the "constant" C introduced in the Von Neumann condition [START_REF] Fowler | Dunes and drumlins[END_REF] should also depend on the space step in order to be able to control the amplification w.r.t. different frequencies and this is not possible for a constant, by definition. Since high frequencies are usually responsible of numerical instabilities, we are going to focus our attention on them. Thereafter, the idea is to exhibit numerical stability conditions to ensure the validity of simulations. We then seek numerical stability criteria such that the amplification factor satisfies:

∀|k| ≥ k 0 , |g(k)| ≤ 1, (8)
where k 0 is some threshold frequency. To ensure this inequality, we will exhibit two sufficient conditions. The first one is rather unusual: it imposes to the space step δx to be smaller than a given positive constant which depends on the ratio ǫ/η of local diffusion to non-local anti-diffusion. We will in fact check numerically that this condition is not necessary. The second one looks more familiar. It is a classical CFL-type condition modified by a η δt/δx 4/3 term, which stems from the nonlocal operator. We will see in the numerical simulations that this condition is both necessary and sufficient to ensure numerical stability.

For a comprehensive study, we also carry out an error analysis: we compute the truncation and phase errors of several finite difference schemes. We finally investigate the convergence of these schemes.

The remaining of this chapter is organized as follows: in the next section, we present finite difference schemes with some discrete version for the fractional derivative and we study the continuous amplifica-tion factor of the linearized Fowler model. Sections 3 and 4 are, respectively, devoted to the stability and error analysis. The paper ends with some remarks in section 5.

Preliminaries

Finite difference approximations

The spatial discretization is given by a set of points x j ; j = 0, ..., N and the discretization in time is represented by a sequence of times t 0 = 0 < ... < t n < ... < T . For the sake of simplicity we will assume constant step sizes δx and δt in space and time respectively. The discrete solution at a point will be represented by u n j ≈ u(t n , x j ). The schemes consist in computing approximate values u n j of solution to (1) on [nδt, (n + 1)δt[×[jδx, (j + 1)δx[ for n ∈ N and j ∈ N thanks to the following relation:

u n+1 j -u n j δt + F (u n j-1 , u n j , u n j+1 ) -ǫ u n j+1 -2u n j + u n j-1 δx 2 + η I δx [u n ] j = 0, (9)
where I δx and F are, respectively, the discretizations of the nonlocal and nonlinear terms. Note that the Laplacian term is discretized using centred finite difference approximation. We begin by considering two discretizations I 1 δx , I 2 δx for the operator I corresponding to formulae (2) and ( 5), respectively. In both cases, we use a basic quadrature rule on the mesh ([jδx, (j + 1)δx[) j∈N to approximate each integral and we use a finite difference approximation of the derivative:

I 1 δx [ϕ] j = δx -4/3 +∞ l=1 l -1/3 (ϕ j-l+1 -2ϕ j-l + ϕ j-l-1 ) , (10)
I 2 δx [ϕ] j = 4 9 δx -4/3 +∞ l=1 l -7/3 ϕ j-l -ϕ j + l ϕ j+1 -ϕ j-1 2 . ( 11
)
Let us remark that we begin the sums at l = 1 in order to avoid the singularity of 1/|z| 1/3 and 1/|z| 7/3 at z = 0. We will comment later on the truncation of the series, see Section 4. Let us simply note that if ϕ j = 0 for all j < 0 then the series [START_REF] Kouakou | Evolution of a model dune in a shear flow[END_REF] is in fact a finite sum. Since the spatial mesh is given by ([jδx, (j + 1)δx[) j∈N , we will indeed assume that ϕ j = 0 for all j < 0.

Remark 2. Using fractional calculus, we could also consider, for any causal function ϕ, the standard Grünwald-Letnikov formula for the fractional derivative I. Indeed, using the expression (3), I can be approximated by the following two formulae

I 3 δx [ϕ] j = Γ(2/3) δx 4/3 l≥0 (-1) l 4/3 l ϕ j-l = Γ(2/3) δx 4/3 l≥0 Γ(l -4/3) Γ(l + 1)Γ(-4/3) ϕ j-l , (12)
and

I 3 δx [ϕ] j = Γ(2/3) δx 4/3 l≥0 -2/3 l (ϕ j-l+1 -2ϕ j-l + ϕ j-l-1 ), (13)
where, for all α > 0 and k ∈ N we denote by α k the binomial coefficient defined by

α k := α(α -1) . . . (α -k + 1) k! = (-1) k Γ(k -α) Γ(-α)Γ(k + 1)
and p k denotes the negative binomial given by

p k = p(p + 1) • • • (p + k -1) k! = (-1) k -p k .
For more details about Grünwald-Letnikov derivatives, we refer the reader to the book [START_REF] Podlubny | An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications[END_REF].

To analyze the stability of the discrete problem (9) using Fourier analysis, we investigate the following linearized explicit scheme

u n+1 j -u n j δt + v u n j -u n j-1 δx -ǫ u n j+1 -2u n j + u n j-1 δx 2 + η I δx [u n ] j = 0, ( 14
)
where v is a positive constant.

Remark 3.

In the case where we consider that v is a non-positive constant, ∂ x u is discretized using a downstream finite difference approximation and so F is given by

F (u n j-1 , u n j , u n j+1 ) = v u n j+1 -u n j δx .
Therefore, taking into account the discretization (10), the numerical scheme is written as follows:

u n+1 j = ǫ δt δx 2 u n j+1 + 1 - vδt δx -2 ǫ δt δx 2 u n j + v δt δx + ǫ δt δx 2 u n j-1 - η δt δx 4/3 +∞ l=1 l -1/3 u n j-l+1 -2u n j-l + u n j-l-1 , (15)
and since

+∞ l=1 l -1/3 u n j-l+1 -2u n j-l + u n j-l-1 = +∞ l=2 (l + 1) -1/3 -2l -1/3 + (l -1) -1/3 u n j-l -u n j -(2 -2 -1/3 )u n j-1 ,
the numerical scheme (15) reads

u n+1 j = ǫ δt δx 2 u n j+1 + 1 - v δt δx -2 ǫ δt δx 2 - η δt δx 4/3 u n j + v δt δx + ǫ δt δx 2 + (2 -2 -1/3 ) η δt δx 4/3 u n j-1 - η δt δx 4/3 +∞ l=2 (l + 1) -1/3 -2l -1/3 + (l -1) -1/3 u n j-l . (16)
Considering now the discretization (11), the numerical scheme ( 14) can be written as follows:

u n+1 j = ǫ δt δx 2 u n j+1 + (1 - v δt δx -2 ǫ δt δx 2 )u n j + ( v δt δx + ǫ δt δx 2 )u n j-1 - 4 9 η δt δx 4/3 +∞ l=1 l -7/3 u n j-l -u n j + l u n j+1 -u n j-1 2 .
Recall that the Riemann zeta function, for Re(s) > 1

ζ(s) = ∞ n=1 n -s . Since +∞ l=1 l -7/3 u n j-l -u n j + l u n j+1 -u n j-1 2 = 1 2 ζ( 4 3
)

u n j+1 -ζ( 7 3
)

u n j - 1 2 ζ( 4 3
)

-1 u n j-1 + +∞ l=2 l -7/3 u n j-l ,
with ζ( 43 ) ≈ 3.601, ζ( 73 ) ≈ 1.415, the numerical scheme reads

u n+1 j = ǫ δt δx 2 - 4 9 η δt δx 4/3 1 2 ζ( 4 3
)

u n j+1 + 1 - v δt δx -2 ǫ δt δx 2 + 4 9 η δt δx 4/3 ζ( 7 3
)

u n j + v δt δx + ǫ δt δx 2 + 4 9 η δt δx 4/3 ( 1 2 ζ( 4 3 ) -1) u n j-1 - 4 9 η δt δx 4/3 +∞ l=2 l -7/3 u n j-l . (17)
Remark 4. If the Fowler equation (1) satisfied the maximum principle, a classical way to get sufficient conditions for the L ∞ -stability of the scheme would be to ensure that u n+1 j is a convex combination of (u n j ) j∈N . Though one can easily check that all coefficients sum up to 1, we remark that (l + 1) -1/3 -2l -1/3 + (l -1) -1/3 > 0 because the function x → x -1/3 is convex and -4 9 η δt δx 4/3 l -7/3 < 0 for all l > 1. Thus, u n+1 j is not a convex combination of (u n j ) j∈N . To get conditions of numerical stability we have to rely on the Von Neumann method.

The continuous amplification factor

In this section, we are going to study the amplification factor of the following equation

∂ t u(t, x) + v ∂ x u(t, x) -ǫ ∂ 2 xx u(t, x) + η I[u(t, •)](x) = 0. (18)
Then, u(t, x) = e ikx+σt is a solution to (18) if and only if the following dispersion relation is satisfied

σ + ivk + ǫk 2 -η|k| 4/3 1 2 Γ( 2 3 ) 1 -i √ 3 sign(k) = 0,
where k ∈ R and σ ∈ C. Indeed, we have

u t (t, x) = σu(t, x), u x (t, x) = iku(t, x), u xx (t, x) = -k 2 u(t, x)
and

I[u(t, •)](x) = +∞ 0 |ξ| -1/3 (-k 2 )e ik(x-ξ)+σt dξ, = -k 2 u(t, x) +∞ 0 |ξ| -1/3 e -ikξ dξ, = -k 2 u(t, x) +∞ 0 |ξ| -1/3 cos(kξ) dξ + i +∞ 0 |ξ| -1/3 sin(kξ) dξ , = -|k| 4/3 1 2 Γ( 2 3 ) + k|k| 1/3 √ 3 2 Γ( 2 3 ) u(t, x),
where we have used Fresnel integrals. Hence the multiplicative factor which enables to get the solution at the time t n+1 from the solution at the time

t n is G cont (k) = e -δt φ(k) , (19)
where

φ(k) = ǫ k 2 -η 1 2 Γ( 2 3 ) |k| 4/3 + i η √ 3 2 Γ( 2 3 ) k|k| 1/3 + vk . Therefore |G cont (k)| = e -δt(ǫ k 2 -η 1 2 Γ( 2
3 ) |k| 4/3 ) .

k * = 1 3 Γ 2 3 η ǫ 3/2 .
Thereby, the exact continuous amplification is maximum for frequency k * , and its modulus is bigger than 1 only for frequencies in the range (0, k 0 ]. The magnitude of this amplification during one time step will also be proportional to δt. Obviously, this phenomenon affects only the low frequencies in the range (0, k 0 ] and strongly depends on the choice of parameters η and ǫ. This is why the standard definitions of stability are not adapted for this model because they do not take into account the possibility of amplification of certain frequencies. And since high frequencies are usually responsible of numerical instabilities, we are going to focus our attention on the high frequencies which are quickly dampened in Fowler's continuous model in order to exhibit numerical stability conditions.

Stability analysis

The purpose of this section is to study the numerical stability of schemes introduced in the previous section and to exhibit stability criteria. We recall that the numerical stability enables to ensure that the difference between the approximated solution and the exact solution remains bounded for all T > 0 with δx, δt given. To get numerical stability criteria, we consider the Von Neumann or Fourier method. In this approach, we assume that the discrete solution is written in as a single Fourier mode

u n j = ûn k e ikx j , (20)
where k ∈ Z is the wave number. Injecting (20) in the numerical scheme [START_REF] Podlubny | An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications[END_REF], we get

ûn+1 k = g(δx, δt, k)û n k , ( 21
)
where g is the discrete amplification factor. In what follows, for simplicity, we denote indifferently g(δx, δt, k) = g(δx, δt, θ), where θ = kδx.

Remark 5. Note that due to the aliasing phenomenon it is enough to study the discrete amplification factor for θ ∈ [0, π].

Following the previous discussion concerning the notion of numerical stability (see Section 2), we introduce the following definition: Definition 1. We say that a numerical scheme which approximates the linearized Fowler equation problem is stable if the high frequencies are strongly stable that is to say:

∃ 0 < θ 0 < π such that ∀θ ∈ (θ 0 , π], |g(δx, δt, θ)| < 1,
where g is the discrete amplification factor. Proof. We can easily check this property, see Figure 2. Proposition 1. The finite difference scheme [START_REF] Podlubny | An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications[END_REF] is stable in the sense of Definition 1 if δx and δt satisfy the following conditions:

• For I 1 δx : v δt δx + 2 ǫ δt δx 2 + (2 -2 -1/3 ) η δt δx 4/3 ≤ 1, (22)
• For I 2 δx : v δt δx + 2 ǫ δt δx 2 + 4 9 ζ( 4 3 ) -1 η δt δx 4/3 ≤ 1, (23)
and if moreover, the space-step δx is small enough in order that • For I 1 δx :

(1 -2 -1/3 ) η δt δx 4/3 ≤ 2 ǫ δt δx 2 sin 2 ( θ 0 2 ), (24)
• For I 2 δx :

4 9 ζ( 7 3
)

-1 + ζ( 4 3 ) η δt δx 4/3 ≤ 2 ǫ δt δx 2 sin 2 ( θ 0 2 ), ( 25
)
where θ 0 designates the stability threshold frequency.

Proof. For I 1 δx . For the numerical scheme (16), the amplification factor is given by:

g 1 (δx, δt, θ) = 1 - v δt δx -2 ǫ δt δx 2 (1 -cos θ) - η δt δx 4/3 + v δt δx + (2 -2 -1/3 ) η δt δx 4/3 e -iθ - η δt δx 4/3 ∞ l=2 (l + 1) -1/3 -2l -1/3 + (l -1) -1/3 e -ilθ , (26)
where θ = kδx. Since, for all

N ∈ N N l=2 (l + 1) -1/3 -2l -1/3 + (l -1) -1/3 = (N + 1) -1/3 -N -1/3 -2 -1/3 + 1, then +∞ l=2 (l + 1) -1/3 -2l -1/3 + (l -1) -1/3 = 1 -2 -1/3 > 0. (27)
Thus, from (26), to have |g 1 (δx, δt, θ)| ≤ 1 it is sufficient to have

1 - v δt δx -4 sin 2 ( θ 2 ) ǫ δt δx 2 - η δt δx 4/3 + v δt δx + (2 -2 -1/3 ) η δt δx 4/3 e -iθ ≤ 1 - η δt δx 4/3 (1 -2 -1/3 ), (28)
where we assume that η δt

δx 4/3 (1 -2 -1/3 ) < 1. (29)
Next from Lemma 1, (28) is satisfied if and only if we have

       1 -4 ǫ δt δx 2 sin 2 ( θ 2 ) + η δt δx 4/3 (1 -2 -1/3 ) ≤ 1 -(1 -2 -1/3 ) η δt δx 4/3 , 1 -2 v δt δx -4 ǫ δt δx 2 sin 2 ( θ 2 ) -(3 -2 -1/3 ) η δt δx 4/3 ≥ -1 -(1 -2 -1/3 ) η δt δx 4/3 .
A sufficient condition is then

(1 -2 -1/3 ) η δt δx 4/3 ≤ 2 ǫ δt δx 2 sin 2 ( θ 2 ), (30)
v δt δx + 2 ǫ δt δx 2 + (2 -2 -1/3 ) η δt δx 4/3 ≤ 1.
Let 0 < θ 0 < π. Then, for all θ ∈ (θ 0 , π], condition (30) can be rewritten as

(1 -2 -1/3 ) η δt δx 4/3 ≤ 2 ǫ δt δx 2 sin 2 ( θ 0 2 ).
Therefore, the numerical scheme [START_REF] Podlubny | An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications[END_REF] with the discretization I 1 δx is stable in the sense of Definition 1 if the space and time steps δt, δx satisfy the following conditions

δx 2/3 ≤ 2 (1 -2 -1/3 ) sin 2 ( θ 0 2 ) ǫ η , (31)
v δt δx + 2 ǫ δt δx 2 + (2 -2 -1/3 ) η δt δx 4/3 ≤ 1. (32)
Note that from condition (32), we can see that hypothesis (29) is satisfied.

For I 2 δx . Injecting (20) in (17), the amplification factor g 2 associated to this scheme is given by:

g 2 (δx, δt, θ) := e iθ ǫ δt δx 2 - 1 2 ζ( 4 3 ) 4 9
η δt δx 4/3 + 1 -

v δt δx -2 ǫ δt δx 2 + 4 9 ζ( 7 3
)

η δt δx 4/3 + e -iθ v δt δx + ǫ δt δx 2 + 4 9 ( 1 2 ζ( 4 3
) -1) η δt δx 4/3 -4 9

η δt δx 4/3 l≥2 l -7/3 e -iθl , = 1 - v δt δx -4 ǫ δt δx 2 sin 2 ( θ 2 ) + 4 9 ζ( 7 3
)

η δt δx 4/3 -i 4 9 ζ( 4 3
) η δt δx 4/3 sin θ

+ v δt δx - 4 9
η δt δx 4/3 e -iθ -4 9

η δt δx 4/3 +∞ l=2 l -7/3 e -iθl , = 1 - v δt δx -4 ǫ δt δx 2 sin 2 ( θ 2 ) + 4 9 ζ( 7 3
)

-ζ( 4 3
) cos θ η δt δx 4/3

+ v δt δx + 4 9 ζ( 4 3
) -1 η δt δx 4/3 e -iθ -4 9

η δt δx 4/3 +∞ l=2 l -7/3 e -iθl . (33)
Since

+∞ l=2 l -7/3 = ζ( 7 3
) -1 ≈ 0.415,

from (33), |g 2 (δx, δt, θ)| ≤ 1 if 1 - v δt δx -4 ǫ δt δx 2 sin 2 ( θ 2 ) + 4 9 ζ( 7 3
)

-ζ( 4 3
) cos θ η δt

δx 4/3 + v δt δx + 4 9 ζ( 4 3
)

-1 η δt δx 4/3 e -iθ ≤ 1 - 4 9 ζ( 7 3
) -1 η δt δx 4/3 , (34

)
where we assume that

4 9 ζ( 7 3
)

-1 η δt δx 4/3 < 1. (35)
From Lemma 1, we have that (34) is satisfied if and only if

       1 -4 ǫ δt δx 2 sin 2 ( θ 2 ) + 4 9 ζ( 7 3
)

-1 + ζ( 4 3 )(1 -cos θ) η δt δx 4/3 ≤ 1 - 4 9 ζ( 7 3
)

-1 η δt δx 4/3 , 1 -2 v δt δx -4 ǫ δt δx 2 sin 2 ( θ 2 ) + 4 9 ζ( 7 3 ) + 1 -ζ( 4
3 )(1 + cos θ) η δt δx 4/3 ≥ -1 + 4 9 ζ( 7 3
) -1 η δt δx 4/3 .

A sufficient condition is then

4 9 ζ( 7 3
)

-1 + ζ( 4 3 ) η δt δx 4/3 ≤ 2 ǫ δt δx 2 sin 2 ( θ 2 ), (36)
v δt δx + 2 ǫ δt δx 2 + 4 9 ζ( 4 3 ) -1 η δt δx 4/3 ≤ 1. (37)
Let 0 < θ 0 < π. Then, for all θ ∈ (θ 0 , π], condition (36) is rewritten as

4 9 ζ( 7 3
)

-1 + ζ( 4 3 ) η δt δx 4/3 ≤ 2 ǫ δt δx 2 sin 2 ( θ 0 2 ), (38)
where ζ( 73 ) -1 + ζ( 43 ) ≈ 4.02. Note again that from condition (37), we can see that hypothesis (35) is satisfied.

Notations. We will denote by CF L 1 mod and CF L 2 mod the following modified Courant-Friedrichs-Lewy conditions

CF L 1 mod = v δt δx + 2 ǫ δt δx 2 + (2 -2 -1/3 ) η δt δx 4/3 ≤ 1, CF L 2 mod = v δt δx + 2 ǫ δt δx 2 + 4 9 ζ( 4 3
) -1 η δt δx 4/3 ≤ 1.

Some remarks.

1. Condition (22) (resp. ( 23)) can be seen as an extension of the classical CFL condition with in addition the anti-diffusive term η δt δx 4/3 . This criterion is not more restrictive than the usual condition of stability without the nonlocal operator which corresponds to the linearized Burgers equation with viscous term. This condition is very restrictive on the space and time steps in particular because of the term ǫ δt δx 2 which stems from the explicit discretization of the Laplacian. In order to have less restrictive conditions, we can implicit some terms. For instance, if we decide to implicit the nonlocal and the Laplacian terms, condition (22) (resp. ( 23) ) is reduced to v δt δx < 1.

We find again the well-known CFL condition.

Figure 3 shows the behaviour of amplification factors for I 1 δx and I 2 δx . We can see, for I 1 δx , that the maximal value of δt which ensures the numerical stability is δt max ≈ 0.042 and that for this value we have CF L 1 mod ≈ 0.99. discretization I 1 δx as a function of θ. We can notice that the high frequencies are strongly amplified. This phenomenon illustrates the numerical instability because high frequencies should be quickly dampened. Figure 5 shows that the low frequencies are slowly amplified. This phenomenon is not due to the instability of numerical schemes but stems from the model. In Tables 1, 2 and 3 (see Section 4), we have studied the quotient |g i | |Gcont| , i = 1, 2. We can see that globally the discrete schemes dampen more than the continuous problem when the stability conditions ( 22) and ( 23) are satisfied. represents the amount of nonlocal anti-diffusion while the term proportional to ǫ δt δx 2 corresponds to the amount of classical diffusion. Both conditions simply mean that, for frequencies above the threshold θ 0 , diffusion should control nonlocal anti-diffusion. We can see that conditions (24) and (25) cannot be satisfied for low frequencies. Indeed, for θ 0 close to 0, these criteria impose to the space step to vanish, which is not possible. Let us note that this is coherent because the low frequencies are not "strongly stable", they are slowly amplified by the continuous problem. We can see in Figure 6 that condition (24) is not necessary. Indeed, if we choose the threshold θ 0 = π/2 "large enough", condition (24) reads

δx ≤ 0.25, (39)
and we have plotted |g 1 | in function of θ for δx = 0.5 which does not satisfy the condition (39) but we can still notice that the numerical scheme is stable. All numerical simulations that we performed confirm this statement. This leads us to think that condition (24) (resp. ( 25)) is too pessimistic. In fact to estimate the magnitude of sums

∞ l=2 (l + 1) -1/3 -2l -1/3 + (l -1) -1/3 e -ilθ ,
(resp. +∞ l=2 l -7/3 e -ilθ ), we just controlled the sum of the modulus

∞ l=2 (l + 1) -1/3 -2l -1/3 + (l -1) -1/3 , (resp.
+∞ l=2 l -7/3 ). In this manner, we probably miss some cancellation effect of the e -ilθ . But we could not find any other way to estimate these polylogarithm series.

3.

Finally, in practice, the single condition (22) (resp. ( 23)) can be used to ensure the numerical stability of the scheme ( 14) with I 1 δx (resp. I 2 δx ). We saw in Figures 5 and6 that the scheme with the discretization I 1 δx is stable if condition ( 22) is satisfied. Figure 7 shows that the high frequencies are amplified, when condition (23) is violated. This phenomenon is only due to numerical instability because the continuous problem quickly dampens the high frequencies.

Error analysis 4.1 Truncation error

In this section, we analyze the truncation error. Finite difference scheme ( 14) is consistent with the linearized partial differential equation if for any smooth function φ(t, x) the local error E δt,δx satisfies

E δt,δx = P φ -P i δt,δx φ → 0, (40)
as δt, δx → 0 with

     P φ = φ t + v φ x -ǫ φ xx + η I[φ], P i δt,δx φ = φ n+1 j -φ n j δt + v φ n j -φ n j-1 δx -ǫ φ n j+1 -2φ n j + φ n j-1 δx 2 + η I i δx [φ],
for i = 1, 2.

Remark 6. The practical implementation of the schemes requires to make some truncations. First, we consider a bounded domain [0, T ] × [0, D] and to simplify, we also assume that δt = T /N δt and δx = D/N δx for some integers N δx and N δt . Another truncation concerns the integral operator for the nonlocal term I. We replace +∞ 0 with A 0 and in the finite difference approximations (10) and (11) series ∞ l=1 are replaced with partial sums A δx l=1 , where A = A δx δx. However, the truncation parameter A has to be chosen judiciously. A "short memory" principle has been investigated to choose this parameter. This principle is based on the fact that terms l -7/3 and l -1/3 in discretizations [START_REF] Kouakou | Evolution of a model dune in a shear flow[END_REF] and [START_REF] Lagree | Stability of an erodible bed in various shear flows[END_REF] decrease with l therefore, we have to take into account the behaviour of ϕ(x) only in the recent past, i.e. in the interval [x -L, x], where L > 0 is called the "memory length". Finally, the use of the short-memory principle leads to the simple replacement of +∞ l=1 by A δx l=1 , where

A δx = [ L δx ] [14]
. Note that the truncation parameter A also strongly depends on the discretization of the nonlocal term I because l -7/3 decreases more quickly than l -1/3 . For the sake of simplicity, we will denote by A the truncation parameter for the discretizations [START_REF] Kouakou | Evolution of a model dune in a shear flow[END_REF] and [START_REF] Lagree | Stability of an erodible bed in various shear flows[END_REF].

Proposition 2 (Local error). The local error of the numerical scheme [START_REF] Podlubny | An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications[END_REF] satisfies:

• For I 1 δx : |E 1 δt,δx,A | ≤ O(δt) + O(δx 2/3 ) + O(A -1/3 ) + O A -1/3 δx + O A 2/3 δx 2 . ( 41
)
• For I 2 δx :

|E 2 δt,δx,A | ≤ O(δt) + O(δx 2/3 ) + O(A -1/3 ) + O A -4/3 δx + O A 2 δx 3 + O A δx 2 . ( 42
)
Proof. From Taylor series, we have

φ t (t n , x i ) - φ n+1 i -φ n i δt = O(δt), (43)
φ x (t n , x i ) - φ n i -φ n i-1 δx = O(δx), (44)
φ xx (t n , x i ) - φ n i+1 -2φ n i + φ n i-1 δx 2 = O(δx 2 ). ( 45
)
Let us now study the truncation error for the nonlocal term I.

For I 1 δx : We rewrite (2) as follows

I[φ(t n , •)](x i ) = A δx j=1 (j+1/2)δx (j-1/2)δx ξ -1/3 φ xx (t n , x i -ξ) dξ + δx 2 0 ξ -1/3 φ xx (t n , x i -ξ) dξ + A A+ δx 2 ξ -1/3 φ xx (t n , x i -ξ) dξ + +∞ A ξ -1/3 φ xx (t n , x i -ξ) dξ, (46)
and the discretization (10) becomes

I 1 δx [φ(t n , •)] i := A δx j=1 δx ξ -1/3 δx Φ n δx = A δx j=1 (j+1/2)δx (j-1/2)δx ξ -1/3 δx Φ n δx dξ, with ξ δx := jδx and Φ n δx = φ n i-j+1 -2φ n i-j +φ n i-j-1 δx 2
. Using (46), we then get the following relation:

I 1 δx [φ(t n , •)] i -I[φ(t n , •)](x i ) = A δx j=1 (j+1/2)δx (j-1/2)δx ξ -1/3 δx Φ n δx -ξ -1/3 φ xx (t n , x i -ξ) dξ - δx 2 0 ξ -1/3 φ xx (t n , x i -ξ) dξ + A+ δx 2 A ξ -1/3 φ xx (t n , x i -ξ) dξ - +∞ A ξ -1/3 φ xx (t n , x i -ξ) dξ, = T 1 -T 2 + T 3 -T 4 .
Let us study the term T 1 . Since

ξ -1/3 δx Φ n δx -ξ -1/3 φ xx (t n , x i -ξ) = ξ -1/3 δx -ξ -1/3 φ xx (t n , x i -ξ) + ξ -1/3 δx (Φ n δx -φ xx (t n , x i -ξ)) , then T 1 = A δx j=1 (j+1/2)δx (j-1/2)δx ξ -1/3 δx -ξ -1/3 φ xx (t n , x i -ξ) dξ + A δx j=1 (j+1/2)δx (j-1/2)δx ξ -1/3 δx (Φ n δx -φ xx (t n , x i -ξ)) dξ, = T 1,1 + T 1,2 .
By the mean value theorem applied to z → |z| -1/3 , we have for all ξ ∈ [(j

-1 2 )δx; (j + 1 2 )δx] |ξ -1/3 δx -ξ -1/3 | ≤ sup z∈[(j-1 2 )δx;(j+ 1 2 )δx] | 1 3 z -4/3 | |ξ δx -ξ|, ≤ 1 3 |(j - 1 2 )δx| -4/3 |ξ δx -ξ|, ≤ 1 6 |(j - 1 2
)δx| -4/3 δx.

Thus, integrating over [(j -1 2 )δx; (j + 1 2 )δx] we get |T 1,1 | ≤ C δx 2/3 A δx j=1 1 (j -1/2) 4/3 ≤ C δx 2/3 , (47) because j≥1 1
(j-1/2) 4/3 < +∞ and C is a positive constant which depends on ||φ xx || L ∞ ((0,T )×R) . Moreover, by classical midpoint quadrature rule

(j+1/2)δx (j-1/2)δx φ xx (t n , x i -ξ) dξ = δx φ xx (t n , x i -ξ δx ) + δx 3 24 φ 4x (t n , x i -η j ), (48)
with

η j ∈ [(j -1/2)δx; (j + 1/2)δx] then T 1,2 = A δx j=1 (jδx) -1/3 δxΦ n δx - (j+1/2)δx (j-1/2)δx φ xx (t n , x i -ξ) dξ , = A δx j=1 (jδx) -1/3 δxΦ n δx -δx φ xx (t n , x i -ξ δx ) - δx 3 24 φ 4x (t n , x i -η j ) .
From Taylor series, we have

φ xx (t n , x i -ξ δx ) = Φ n δx - δx 2 12 φ 4x (t n , x i -ξ δx ) + O(δx 4 ), (49)
thus we obtain

T 1,2 = A δx j=1 (jδx) -1/3 δx 3 12 φ 4x (t n , x i -ξ δx ) - δx 3 24 φ 4x (t n , x i -η j ) + O(δx 5 ) .
We finally get

|T 1,2 | ≤ O(δx 8/3 ) A δx j=1 1 j 1/3 = O(δx 8/3 )   1 + A δx j=2 1 j 1/3   , = O(δx 8/3 ) 1 + A δx 0 y -1/3 dy , which implies that |T 1,2 | ≤ O(δx 8/3 ) + O(A 2/3 δx 2 ), (50)
because A δx = A δx . We next control the term T 2 by

|T 2 | ≤ ||φ xx || L ∞ ((0,T )×R) δx 2 0 |ξ| -1/3 dξ = Cδx 2/3 . (51)
We estimate T 3 as follows:

|T 3 | ≤ ||φ xx || L ∞ ((0,T )×R) A+δx A |ξ| -1/3 dξ, ≤ ||φ xx || L ∞ ((0,T )×R) A -1/3 δx. (52)
Finally, using an integration by parts, the term T 4 is written as

T 4 = +∞ A ξ -1/3 φ xx (t n , x i -ξ) dξ, = -A -1/3 φ x (t n , x i -A) + 1 3 +∞ A ξ -4/3 φ x (t n , x i -ξ) dξ, hence, we obtain |T 4 | ≤ CA -1/3 , (53)
where C is a positive constant which depends on ||φ x || L ∞ ((0,T )×R) . Hence, using relations (43), ( 44), ( 45), ( 47), ( 50), (51), ( 52) and (53), we obtain

|E 1 δx,δt,A | = |P φ(t n , x i ) -P 1 δt,δx φ| ≤ O(δx 2/3 ) + O(δt) + O(A -1/3 ) + O(A 2/3 δx 2 ) + O A -1/3 δx ,
which completes the proof for I 1 δx . For I 2 δx : As previously, we rewrite (5) as

I[φ(t n , •)](x i ) = A δx j=1 (j+1/2)δx (j-1/2)δx Φ(t n , ξ) |ξ| -7/3 dξ + δx 2 0 Φ(t n , ξ) |ξ| -7/3 dξ + A A+ δx 2 Φ(t n , ξ) |ξ| -7/3 dξ + +∞ A Φ(t n , ξ) |ξ| -7/3 dξ,
with Φ(t n , ξ) = 4 9 (φ(t n , x i -ξ) -φ(t n , x i ) + φ x (t n , x i ) ξ) and the approximated integral (11) becomes

I 2 δx [φ(t n , •)] i := A δx j=1 (j+1/2)δx (j-1/2)δx ξ -7/3 δx Φ n δx dξ, with Φ n δx = 4 9 φ n i-j -φ n i + φ n i+1 -φ n i-1
2 j and ξ δx = jδx. Let us now estimate the error on the nonlocal term.

I 2 δx [φ(t n , •)] i -I[φ(t n , •)](x i ) = A δx j=1 (j+1/2)δx (j-1/2)δx Φ n δx ξ -7/3 δx -Φ(t n , ξ)|ξ| -7/3 dξ - δx 2 0 Φ(t n , ξ) |ξ| -7/3 dξ + A+ δx 2 A Φ(t n , ξ) |ξ| -7/3 dξ - +∞ A Φ(t n , ξ) |ξ| -7/3 dξ, = T 1 -T 2 + T 3 -T 4 .
Let us study the term T 1 . As previously for I 1 δx , we rewrite T 1 as

T 1 = A δx j=1 (j+1/2)δx (j-1/2)δx (ξ -7/3 δx -|ξ| -7/3 )Φ(t n , ξ) + ξ -7/3 δx (Φ n δx -Φ(t n , ξ)) dξ, = T 1,1 + T 1,2 .
By the mean value theorem applied to z → |z| Next, by Taylor-Lagrange formula, we have

|T 1,1 | = 4 9 A δx j=1 (j+1/2)δx (j-1/2)δx |ξ -7/3 δx -|ξ| -7/3 | |φ(t n , x i -ξ) -φ(t n , x i ) + φ x (t n , x i ) ξ| , ≤ C A δx j=1 (j+1/2)δx (j-1/2)δx δx -7/3 1 (j -1/2) 10/3 ||φ xx || L ∞ ((0,T )×R) ξ 2 2 dξ, ≤ Cδx -7/3 A δx j=1 1 (j -1/2) 10/3 (j+1/2)δx (j-1/2)δx ξ 2 dξ, ≤ Cδx -7/3 A δx j=1 1 (j -1/2) 10/3 δx 3 j 2 + δx 3 , ≤ Cδx 2/3 , because ∞ j=1 j 2
(j-1/2) 10/3 < ∞. C denotes a positive constant which depends on ||φ xx || L ∞ ((0,T )×R) and may vary from line to line. Moreover, using again midpoint quadrature rule, we have

(j+1/2)δx (j-1/2)δx Φ(t n , ξ) dξ = δx Φ(t n , jδx) + δx 3 24 Φ xx (t n , η j ), with η j ∈ [(j -1 2 )δx, (j + 1 2 )δx]. Hence, T 1,2 = A δx j=1 ξ -7/3 δx δx Φ n δx -δx Φ(t n , jδx) - δx 3 24 Φ xx (t n , η j ) .
But using again Taylor expansion, we get

φ x (t n , x i ) = φ(t n , x i+1 ) -φ(t n , x i-1 ) 2δx - δx 2 6 φ 3x (t n , x i ) + O(δx 3 ),
and so

Φ n δx -Φ(t n , jδx) = 4 9 φ n i-j -φ n i + φ n i+1 -φ n i-1 2 j -φ(t n , x i -jδx) + φ(t n , x i ) -φ x (t n , x i )jδx , = 4 9 ξ δx δx 2 6 φ 3x (t n , x i ) + jO(δx 4 ) .
Thus,

T 1,2 = C A δx j=1 ξ -7/3 δx ξ δx δx 3 6 φ 3x (t n , x i ) + jO(δx 5 ) - δx 3 24 Φ xx (t n , η j ) , = C A δx j=1 ξ -7/3 δx ξ δx δx 3 6 φ 3x (t n , x i ) + jO(δx 5 ) - 4 9 δx 3 24 φ xx (t n , x i -η j ) .
We finally get

|T 1,2 | ≤ C A δx j=1 j -4/3 δx 5/3 + jO(δx 5 ) + δx 3 , = O(δx 5/3 ) + O A 2 δx 3 + O A δx 2 ,
where C is a positive constant which depends on ||φ 3x || L ∞ ((0,T )×R) and ||φ xx || L ∞ ((0,T )×R) .

Let us now study T 2 . Using Taylor-Lagrange formula, we have

|T 2 | ≤ C||φ xx || L ∞ ((0,T )×R) δx 2 0 ξ 2 |ξ| 7/3 dξ, ≤ O(δx 2/3 ),
where C is a positive constant. Let us next consider T 3 :

T 3 = A+ δx 2 A |ξ| -7/3 Φ(t n , ξ) dξ, = 4 9
A+ δx 2 A [φ(t n , x i -ξ) -φ(t n , x i ) + φ x (t n , x i )ξ] |ξ| -7/3 dξ.
Then

|T 3 | ≤ C A+ δx 2 A (|ξ| -7/3 + |ξ| -4/3 ) dξ, ≤ C A -7/3 + A -4/3 δx, ≤ O δx A -4/3 ,
with C a positive constant which depends on ||φ|| L ∞ ((0,T )×R) and ||φ x || L ∞ ((0,T )×R) . And since

|T 4 | ≤ 4 9 +∞ A |φ(t n , x i -ξ) -φ(t n , x i )| |ξ| -7/3 dξ + 4 9 +∞ A |φ x (t n , x i )| |ξ| -4/3 dξ, ≤ CA -1/3 ,
where C is a positive constant which depends on ||φ|| L ∞ ((0,T )×R) and ||φ x || L ∞ ((0,T )×R) , we finally get

|E 2 δx,δt,A | := |P φ(t n , x i ) -P 2 δt,δx φ| ≤ O(δt) + O(δx 2/3 ) + O(A -1/3 ) + O A -4/3 δx + O A 2 δx 3 + O A δx 2 .
The proof of this proposition is now completed.

Remark 7.

From previous Proposition, we can see that the numerical scheme (14) with I 1 δx (resp. I 2 δx ) is consistent if δx << A -1/3 (resp. δx << A -2/3 ).

Convergence experiments.

In this section, we investigate the convergence using numerical simulations. Despite much effort we are unable to prove theoretically the convergence of the numerical solution towards the exact continuous solution. Indeed, the Lax procedure "stability + consistence = convergence" cannot be applied here due to the instability of low frequencies. In what follows, l 1 -norm is used to measure the accuracy of approximated solutions. Thus, we analyze the following error

E 1 = 1 N N n=0 (|u 1 j (T ) -u 2 j (T )|), (54)
where u 1 , u 2 are, respectively, computed for space steps δx/2 and δx/4, until a final time T . Figure 9 shows the numerical convergence rates obtained with the initial data displayed in Figure 8. These rates were obtained using δx = 10 -1 , 10 -2 , 10 -3 , 10 -4 . We plot the logarithm of the error E 1 versus the logarithm of δx.

Phase error

Numerical schemes produce, according to cases, results ahead of or delayed w.r.t exact solutions. In this section, we are interested in the error made on the velocity introduced by the discretization. Let us first note that, in addition to the anti-diffusive effect, the nonlocal term is also responsible of the motion of the initial data. Indeed, we saw in Section 2 that the continuous amplification factor has an imaginary part e

-i η √ 3 2 Γ( 2
3 ) k|k| 1/3 + vk δt . Therefore, the advection term is not the unique factor of displacement. It is the error on the argument -δt(v k + √ 3 2 Γ( 2 3 )η k|k| 1/3 ) that causes the phase error. To evaluate this error, we rewrite the discrete amplification factor g j introduced in Section 3 as

g j = |g j |e -iθ d j ,
for j = 1, 2, where θ d j is the argument of the discrete amplification factor g j . The phase lag during one time step is then given by

E j = (v k + √ 3 2 Γ( 2 3 )η k |k| 1/3 )δt -θ d j .
Thus, if E j is positive, the numerical wave goes slower than the physical wave and it goes faster if E j is negative. We have computed, in Tables 1, 2 and 3, the phase delay after one oscillation

∆ j = 1 - θ d j (v k + √ 3 2 Γ( 2 3 )η k |k| 1/3 )δt , for j = 1, 2,
for different values of Cr = v δt δx , D f = 2 ǫ δt δx 2 and F o = η δt δx 4/3 . We remark that scheme with the discretization I 2 δx involves a delay larger than the model with I 1 δx .

Discrete amplification vs. continuous amplification factors

Let us define

G 1 = |g 1 | |G cont | , G 2 = |g 2 | |G cont | , for different values of Cr = v δt δx , D f = 2 ǫ δt δx 2 and F o = η δt δx 4/3 .
Results are reported in Tables 1, 2 and 3. We have CF L

i mod = Cr + Df + λ i F o , where λ 1 = 2 -2 -1/3 , λ 2 = 4 9 ζ( 4
3 ) -1 . When the condition CF L i mod > 1 is violated because Cr or D f are close to one, we note that high frequencies θ > π/2 are more amplified by the discrete schemes than by the exact continuous problem. Whereas when the condition CF L i mod > 1 is violated because F o is close to one, high frequencies θ > π/2 are less amplified by the discrete schemes than by the continuous problem, as can be checked in Table 3. This is one unexpected benefit of the discretization scheme in the unfavourable case where nonlocal anti-diffusion is predominant. If we take a closer look at Table 3 we notice that |g 2 | may be greater than one even if CF L 2 mod < 1. This is due to the fact that η being big, the stability threshold frequency θ 0 is close to π, the aliasing limit frequency.

Concluding remarks

We have presented in this work a first investigation of finite differences schemes approximating the Fowler equation. We saw that the anti-diffusive behaviour of the nonlocal term does not enable to consider the classical notion of stability. Nevertheless, considering only the behaviour of the high frequencies (which should be quickly dampened), we exhibit numerical stability criteria which can be used to make simulations. Numerical computations have shown that numerical schemes dampened more than the continuous problem. Finally, consistency property has been proved and convergence of schemes has been investigated.

Cr CF L 1 mod CF L 2 mod θ ∆ 1 ∆ 2 G 1 G 2 0.

Figure 1 :

 1 Figure 1: Behaviour of Re (φ) for η, ǫ fixed. k 0 = 1 2 Γ 2 3 η ǫ 3/2 is the threshold frequency.

Figure 1

 1 Figure 1 shows that the modulus of the continuous amplification factor during one time step is controlled by e α * δt , with α * := -min Re (φ) = -Re φ(k * ) = 4 27

Lemma 1 .

 1 Let a, b ∈ R and d ∈ R + . Then we have ∀θ ∈ [0, 2π], |a + be -iθ | ≤ d if and only if a + |b| ≤ d and a -|b| ≥ -d.

Figure 2 :

 2 Figure 2: Dashed circle (resp. continuous circle) is centred at a (resp. 0) and of radius |b| (resp. d).

Figure 4

 4 Figure3shows the behaviour of amplification factors for I 1 δx and I 2 δx . We can see, for I 1 δx , that the maximal value of δt which ensures the numerical stability is δt max ≈ 0.042 and that for this value we have CF L 1 mod ≈ 0.99. Figure4displays the behaviour of the modulus of the amplification factor with

Figure 3 :

 3 Figure 3: Amplification factors for I 1 δx (blue line) and I 2 δx (dashed line).

Figure 4 :

 4 Figure 4: Amplification factor for I 1 δx with CF L 1 mod ≈ 1.22.

Figure 5 :

 5 Figure 5: Amplification factor I 1 δx with η = 8, v = 1, ǫ = 0.5 and δx = 0.05, δt = 0.001. For these coefficients, CF L 1 mod ≈ 0.94. 2. Conditions (24) and (25) are unusual and deserve some explanations. The term proportional to η δt δx 4/3

Figure 6 :

 6 Figure 6: Amplification factor g 1 for η = v = 1, ǫ = 0.1 and δx = 0.5, δt = 0.01. For these coefficients, we have CF L 1 mod ≈ 0.0584.

Figure 7 :

 7 Figure 7: Gain function for I 2 δx for η = v = ǫ = 1, δx = 0.1, δt = 0.05. For these coefficients, we have CF L 2 mod ≈ 1.3.

- 7 / 3 ,

 73 we have for all ξ ∈ [(j -1/2)δx; (j + 1/2)δx] 1/2) 10/3 .

Figure 8 :

 8 Figure 8: Initial data used for numerical experiments.

Figure 9 :

 9 Figure 9: Convergence in space for I 1 δx (red) and I 2 δx (blue). Dotted line has slope 2/3.

Table 1 :

 1 Dampening and phase error for D f = 0.2 and F o = 0.1

	2	0.5206	0.5156	π/6	0.0082	0.0333 0.9584 0.9788
				π/4	0.0024	0.0573 0.9102 0.9550
				π/2	-0.0715	0.1610 0.6824 0.8394
				3π/4 -0.2684	0.3911 0.3788 0.6626
	0.5	0.8206	0.8156	π/6	-0.0128	0.0104 0.9541 0.9736
				π/4	-0.0433	0.0091 0.9048 0.9452
				π/2 -0.02476 -0.0315 0.7439 0.8152
				3π/4 -0.4733 -0.1628 0.8284 0.6391
	0.9	1.2206	1.2156	π/6	-0.0103	0.0107 0.9870 1.0047
				π/4	-0.0306	0.0128 0.9869 1.0182
				π/2	-0.0781	0.0172 1.1776 1.1522
				3π/4 -0.0210	0.0371 1.8187 1.5053
	D f CF L 1 mod	CF L 2 mod	θ	∆ 1	∆ 2	G 1	G 2
	0.2	0.4206	0.4156	π/6	0.0213 0.0481 0.9654 0.9859
				π/4	0.0303 0.0868 0.9250 0.9707
				π/2	0.0525 0.2562 0.7340 0.9057
				3π/4 0.1721 0.5567 0.4834 0.8487
	0.4	0.6206	0.6156	π/6 -0.0066 0.0216 0.9649 0.9860
				π/4 -0.0358 0.0276 0.9216 0.9701
				π/2 -0.3470 0.0154 0.6750 0.8841
				3π/4 -2.0019 0.0277 0.4153 0.7059
	0.8	1.0206	1.0156	π/6 -0.0673 -0.0361 0.9614 0.9837
				π/4 -0.1989 -0.1169 0.9022 0.9568
				π/2 -3.1203 -1.4458 0.5305 0.6566
				3π/4 -3.8370 -3.6360 5.5254 3.5099

Table 2 :

 2 Dampening and phase error for Cr = 0.1 and F o = 0.1 .6556 0.5123 0.7449 1.2178 1.7707 2.3771 3π/4 0.2981 0.8366 0.3873 0.4083 1.5020 1.5836 3.8781 π 0.3924 1.0000 0.2696 0.2126 1.6583 1.3075 6.1519

	F o CF L 1 mod 0.2 0.5413	CF L 2 mod 0.5312	θ π/6	∆ 1 0.0307 0.0784 0.9455 0.9852 0.9741 1.0150 1.0302 ∆ 2 G 1 G 2 |g 1 | G cont |g 2 |
				π/4	0.0363 0.1360 0.8852 0.9707 0.9180 1.007 1.0371
				π/2 -0.0079 0.3495 0.6236 0.9052 0.6240 0.9057 1.0005
				3π/4 -0.1876 0.6374 0.3216 0.8182 0.2822 0.7180 0.8776
				π	-1.2548 1.0000 0.0574 0.6017 0.0574 0.6017 0.6950
	0.5	0.9031	0.8780	π/6	0.0591 0.1552 0.8992 0.9875 1.0092 1.1084 1.1224
				π/4	0.0650 0.2538 0.8043 0.9758 0.9664 1.1725 1.2016
				π/2 -0.0103 0.5264 0.5235 0.8684 0.7589 1.2590 1.4498
				3π/4 -0.0906 0.7619 0.4215 0.6524 0.6992 1.0824 1.6590
				π	-0.0430 1.0000 0.4202 0.5110 0.7435 0.9042 1.7694
	0.9	1.3857	1.3404	π/6	0.1064 0.2438 0.8602 0.9941 1.0822 1.2512 1.2583
				π/4	0.1361 0.3752 0.7523 0.9762 1.0999 1.4273 1.4621
				π/2	0.2002 0

Table 3 :

 3 Dampening and phase error for Cr = 0.1 and D f = 0.2. [2] ALVAREZ-SAMANIEGO B.; AZERAD P., Existence of travelling-wave and local well-posedness of the Fowler equation, Disc. Cont. Dyn. Syst., Ser. B, 12 (2009), pp. 671-692.

The authors are supported by the ANR MATHOCEAN ANR-08-BLAN-0301-02.


D Bjophthalmol

MD Maria Vadalà

PhD Alfredo Pece

MD Stefano Cipolla

MD Carla Monteleone

MD Giuseppe Fasolino

MD Alessandra Casuccio

BS Salvatore Cillino

Maria Vadalà
email: maria.vadala@unipa.it

Is Ranibizumab effective in stopping the loss of vision for Choroidal Neovascularization in Pathologic Myopia? A Long Term Follow-up Study

Keywords: pathological myopia, choroidal neovascularizaton, ranibizumab, intraocular injection, anti-VEGF drugs

Purpose: To assess the efficacy and safety of Ranibizumab in the treatment of choroidal neovascularization (CNV) caused by pathologic myopia (PM).

Design: Prospective, multicenter, interventional case series.

Methods: 40 of 39 consecutive patients with PM and CNV were treated with "on-demand" intravitreal injection of ranibizumab 0.5 mg. Final best-corrected VA (BCVA) and its change from baseline were the main outcome measures. Changes in Optical Coherence Tomography (OCT) central retinal thickness (CRT) were a secondary outcome.

Results: Mean age was 53 ± 13 years, mean refractive error -13.5 ± 6.5 diopters.

Median follow-up was 13.3 ± 2 months (range 12 -18). Fifteen eyes (37.5%) had previously been treated with photodynamic therapy (PDT). The mean baseline logMAR BCVA (ETDRS vision chart) was 0.68 ± 0.34 (Snellen equivalent 20/131) and 21 ± 16 letters. The final mean logMAR BCVA was 0.27 ± 0.2 (p 0.008) (20/42) and 40.5 ± 14 letters (p 0.01). Mean final VA improved in 82.5% of patients, in 60% by 3 or more lines,(median number of lines gained 2.9). Even six out of seven cases of low vision (≤ 1.1 LogMAR) at the final examination improved vision. Mean OCT CRT reduced from 218±70 μm to 175±46 μm (p 0.02). Age and previous PDT did not influence results (p›0,05). The mean number of injection was 2.8 ± 1.2 (range 1-6). No ocular or systemic side effects were observed.

Conclusion:

Ranibizumab was an effective treatment for stabilizing and improving vision in 92.5 % of myopic CNV in a long term follow-up with a low number of injections.

INTRODUCTION

Choroidal neovascularization (CNV) is the main cause of acute vision loss in pathologic myopia (PM). [START_REF] Yoshida | Myopic choroidal neovascularization[END_REF][START_REF] Soubrane | Choroiodal Neovascularization in pathologic myopia: recent developments in diagnosis and treatment[END_REF]. The natural course of the disease accounts for low visual acuity (20/200) in 44-60% of the patients after 24 months. Laser photocoagulation is no longer used in juxtafoveal CNV, because of large retinal scars increasing over time. [START_REF] Brancato | Photocoagulation scar expansion after laser theraphy for choroidal neovascularization in degenerative myopia[END_REF] Photodynamic therapy (PDT) with verteporfin may stabilize visual acuity (VA) and significantly reduce visual loss in eyes with subfoveal and juxtafoveal CNV PM. [START_REF]Photodynamic therapy of subfoveal choroidal neovascularisation in pathologic myopia with verteporfin. 1-year results of a randomized clinical trial-Vip report N° 1[END_REF][START_REF] Pece | Photodynamic Therapy With Verteporfin for subfoveal Choroidal Neovascularizzation secondary to Pathologic Myopia: A Longterm Study[END_REF][START_REF] Pece | Photodynamic Therapy With Verteporfin for Juxtafoveal Choroidal Neovascularizzation in Pathologic Myopia: A Long-term Follow-up Study[END_REF] Nevertheless the results show a decrease of VA in up to 1/3 of the eyes with subfoveal CNV and only 6% of the eyes gain more than 3 lines in juxtafoveal form. [START_REF] Pece | Photodynamic Therapy With Verteporfin for Juxtafoveal Choroidal Neovascularizzation in Pathologic Myopia: A Long-term Follow-up Study[END_REF] Meanwhile, the widespread use of anti-VEGF drugs in the treatment of CNV due to age-related macular degeneration (AMD) has opened new therapeutical strategies.

Bevacizumab off-label use in CNV related to PM has been evaluated in a few clinical nonrandomized clinical trials with a follow-up period of maximum 12 months follow-up. [START_REF] Ikuno | Intravitreal bevacizumab for choroidal neovascularization attributable to pathological myopia: one-year results[END_REF][START_REF] Gharbiya | Intravitreal bevacizumab treatment for choroidal neovascularization in pathologic myopia: 12-month results[END_REF] Few reports about ranibizumab show short-term efficacy in maintaining vision. [START_REF] Silva | Short-term efficacy of intravitreal ranibizumab for myopic choroidal neovascularization[END_REF][START_REF] Kostantidinis | Intravitreal ranibizumab (Lucentis) for the treatment of myopic choroidal neovascularization[END_REF][START_REF] Mones | Intravitreal ranibizumab for choroidal neovascularization secondary to pathologic myopia: 12-month results[END_REF][START_REF] Lai | Intravitreal ranibizumab for the primary treatment of choroidal neovascularization secondary to pathologic myopia[END_REF] We report the results of a 12-month prospective study in patients with pathologic myopia who were treated with intravitreally administered ranibizumab 0.5 mg for CNV.

MATERIALS AND METHODS

Forty eyes of 39 patients with CNV related to pathologic myopia were consecutively enrolled in this open-label, multicenter, non-randomized, prospective case series. They were treated with an off-label intravitreal administration of ranibizumab 0.5 mg, regardless of lesion location.

Inclusion criteria were: 1) highly myopic eyes, spherical equivalent greater than -6 diopters (D), 2) posterior pole myopic retinal changes (lacquer cracks, chorioretinal atrophy, papillary crescent, posterior staphyloma), 3) fluorescein angiography (FA) detection of the subfoveal or iuxtafoveal CNV, 4) best-corrected visual acuity >20/400 at baseline, 5) duration of symptoms no longer than 4 weeks before enrollment, 6) minimum follow up of 12 months, 7) clear ocular media.

Exclusion criteria included: 1) post inflammatory macular changes, 2) extrafoveal CNV, 3) age-related macular alterations, 4) refractive media opacities, 5) presence of other maculopathies as diabetic retinopathy or retinal vascular occlusion, 6) history of recent myocardial infarction or other thromboembolic events, 7) other previous intravitreal drugs injections, 8) ocular hypertension or glaucoma.

Previous photodynamic therapy was not considered as an exclusion criteria if administered more than three months before enrolling. However no additional PDT session was performed during the follow-up and treatment options included either "on demand" ranibizumab intravitreal injections or observation.

Best corrected visual acuity (BCVA) was recorded using the Early Treatment Diabetic

Retinopathy Study (ETDRS) chart at 4 meters distance using a logMAR scale: results were recorded as logMAR values and the number of letters read. At the baseline and at every examination, all patients underwent the standard ocular examination with ETDRS VA evaluation, biomicroscopy, tonometry, and OCT evaluation (Stratus OCT 4.01, Carl Zeiss Meditec, Inc.,Dublin, CA,USA). Retinal fluorescein angiography (TRC-50/IMAGEnet H 1024 system, TOPCON, Tokio. Japan, and/or Heidelberg Retina Angiograph 2, Heidelberg Engineering, Heidelberg, Germany) was performed at baseline within the 15 days prior to the initial visit, and subsequently, according to the judgment of the ophthalmologist, during the follow-up period.

The Retinal Thickness Map analysis program was applied with the extrapolation of data from six diagonal six-mm radial line scans. The scans were reviewed and manually caliper-assisted measurements were used in case of delineation errors. Standard deviation bigger than 10% of the macular thickness was considered out of range and scans were repeated or rejected. A decrease of 10% thickness in comparison with the baseline was considered as a reduction, while an increase of 10 % as an amount.

The off label use of the drug and its potential risks and benefits were discussed extensively with all the patients and written informed consent was obtained. The study was conducted in accordance with the tenets of the Declaration of Helsinki.

Patients received an initial dose of 0.5 mg (0.05 ml) of ranibizumab using the standard injection procedure in a surgical room. [START_REF] Aiello | Evolving guidelines for intravitreous injections[END_REF] Patients were reexamined the day after treatment. All patients were given detailed instructions and asked to call promptly if any significant changes in vision occurred. The patients were examined on a monthly basis, after the injection for the first three months, and monthly or quarterly thereafter, according to the judgement of the specialist.

Re-treatment was performed if any of the following conditions occurred: decrease of BCVA and/or increase of metamorphopsia, new retinal hemorrhage, increase of retinal/subretinal fluid or retinal thickness on OCT, CNV leakage on FA images. Visual acuity changes from baseline ware defined as improved if an increase of 1 ETDRS line or five letters was registered, or worsened if there was a decrease of 1 line; otherwise VA was considered unchanged, Data are presented as mean, standard deviation (SD), median, and percentage (%).

RESULTS

Thirty-nine patients were enrolled in this study and forty eyes were treated. Twelve patients were men and twenty-seven were women; mean age at the initial examination was 53 years old (SD 13, range 32-80 years). Mean refractive error was -13.5 D (SD 6.5 D, range -6 to -28 D). In twenty-nine eyes the CNV was subfoveal (72.5%) and in eleven eyes it was juxtafoveal (27.5%). Fifteen eyes (37.5%) were previously treated with PDT, twenty-three (57.5%) were not. The mean ETDRS BCVA at baseline was 0.67 logMAR (SD 0.34, range 1.1-0.1, Snellen equivalent 20/131), and 21 letters (SD 16, range 0-50).

The mean follow-up was 13.3 months (SD 2, range 12-18 months). All the patients completed a12-month follow-up, nine of them being followed-up for 18 months. Clinical data at baseline are reported in table 1. The mean BCVA was stable or improved during post-treatment follow-up with a mean final VA of 0.27 logMAR (SD 0.2, range 1.1-0.0, Snellen equivalent 20/41) and 40.5 letters (SD 14, range 0-55) at 12 months after enrolling (p = 0.008 for LogMAR, 0.01 for letters)

(Figure 1).

The final BCVA, compared with baseline, improved in 32 out of 40 eyes (80%), decreased in 3 eyes (7.5 %), and remained stable in 5 eyes (12.5 %). In the eyes showing improvement, the median number of lines gained from baseline was 2.9 (range 1-8), and 19 out of the 32 (59.4%) patients improved by three or more lines: nine eyes (28%) improved by one line, 4 eyes (12.5%) by 2 lines. Three eyes worsened in visual acuity despite therapy: one eye lost one ETDRS line, one eye two lines and the third case lost three lines. They however maintained BCVA between 0.2 and 0.6 LogMAR. Seven eyes had BCVA lower than 1 LogMAR (range 1.1-1.5) at the enrolling time but, at the end of the follow-up period, only one eye was lower (1.5 LogMAR): that was its exact BCVA at the start-up, so it was considered as stabilized.

Mean central macular thickness analyzed by OCT showed a reduction in all the cases, from 218.9 (SD 70) microns to 175 (SD 46) microns (p = 0.02) (Figure 2).

When analyzing visual acuity results in relation to refractive errors, age and lesion size, no statistical accordance was found. Eyes with an initially poor VA (lower than 1.1 LogMAR)

did not statistically differ in visual improvement from those with initially good VA.

Considering the entire group, the mean number of injections was 2.8 (SD 1.2): the number of treatments varied from 1 to a maximun of 6.

The eyes previously treated with PDT were compared with the naïve ones. No ocular or systemic side effects were observed: particularly, neither retinal treas nor detachment occurred during the observation.

Figure 3 shows fluorescein angiography and OCT changes from basal observation to the last examination.

DISCUSSION

Treatment of CNV in PM is still in debate and should be further defined. Currently laser photocoagulation is the standard treatment for extrafoveal CNV [START_REF] Secretan | Long-term visual outcome of choroidal neovascularization in pathologic myopia. Natural history and laser treatment[END_REF], while PDT is the only approved therapy for subfoveal lesions showing significant visual and anatomic benefits [START_REF]Photodynamic therapy of subfoveal choroidal neovascularisation in pathologic myopia with verteporfin. 1-year results of a randomized clinical trial-Vip report N° 1[END_REF][START_REF] Pece | Photodynamic Therapy With Verteporfin for subfoveal Choroidal Neovascularizzation secondary to Pathologic Myopia: A Longterm Study[END_REF][START_REF] Pece | Photodynamic Therapy With Verteporfin for Juxtafoveal Choroidal Neovascularizzation in Pathologic Myopia: A Long-term Follow-up Study[END_REF] . However two-and three -year observational reports [START_REF] Krebs | Choroidal neovascularization in pathologic myopia: three-year result after photodynamic therapy[END_REF][START_REF]Verteporfin therapy for subfoveal choroidal neovascularization in pathologic myopia: 3-years results of an open-label extension following completion of a 2-year randomized clinical trial-Vip report[END_REF] fail to demonstrate the effectiveness of PDT in stabilizing vision or reducing severe loss of vision versus a placebo. Moreover, in juxtafoveal form, about 25% of eyes had experience of a decrease in vision and only 6% of the eyes gained 3 or more lines. [START_REF] Pece | Photodynamic Therapy With Verteporfin for Juxtafoveal Choroidal Neovascularizzation in Pathologic Myopia: A Long-term Follow-up Study[END_REF] Vascular endothelial growth factor A (VEGF-A) has been implicated as a major mediator responsible for subretinal neovascularization in several ocular diseases and it is now the main target for therapy in exudative age-related macular degeneration. Ranibizumab and bevacizumab are the most diffused anti-VEGF drugs for AMD complicated by CNV, giving a pan-vascular endothelial growth factor blocking. [START_REF] Rosenfeld | Ranibizumab for neovascular age-related macular degeneration[END_REF][START_REF] Ciulla | Antivascular endothelial growth factor therapy for neovascular age-related macular degeneration[END_REF][START_REF] Mckay | VEGF and age-related macular degeneration[END_REF] The effect of anti-VEGFdrugs on aqueous VEGF and PEDF levels on other forms of CNV, as in PM, has been demonstrated [START_REF] Chan | Changes in aqueous vascular endothelial growth factor and pigment epithelial-derived factor levels following intravitreal bevacizumab injections for choroidal neovascularization secondary to age-related macular degeneration or pathologic myopia[END_REF], and such therapy is now having wide clinical support. [START_REF] Gunther | Bevacizumab (Avastin) for the treatment of ocular disease[END_REF] Previous papers about bevacizumab in myopic CNV indicate very promising results. [START_REF] Ikuno | Intravitreal bevacizumab for choroidal neovascularization attributable to pathological myopia: one-year results[END_REF][START_REF] Gharbiya | Intravitreal bevacizumab treatment for choroidal neovascularization in pathologic myopia: 12-month results[END_REF] Ikuno et coll. [START_REF] Ikuno | Intravitreal bevacizumab for choroidal neovascularization attributable to pathological myopia: one-year results[END_REF], in a retrospective analysis, enrolled 63 myopic eyes and treated them "on demand", disregarding about CNV clinical form, visual acuity or timing of symptoms at the enrolment. The mean number of injections was 2.4 ± 1.4: 40% of patients had improved 10 vision (27% from three to six lines,13% six or more lines), while 56% remained stable, with a gain or loss of less than three lines. Three patients (5%) lost more than three lines in twelve months, because of enlargement of macular atrophy.

Gharbiya et coll. [START_REF] Gharbiya | Intravitreal bevacizumab treatment for choroidal neovascularization in pathologic myopia: 12-month results[END_REF], in a one-year follow-up prospective study on 20 eyes, demonstrated that 1.25 mg intravitreal bevacizumab can improve vision significantly and, reduce macular thickness at the OCT and angiographic leakage in 95% of cases. In 90% of subjects vision improved by 10 or more letters and in 70% by 15 or more letters. None of the patients experienced a worsening of visual acuity with respect to the baseline condition.The mean number of injections was 4, ranging from 3 to 7, starting with a loading dose.

Very few papers report clinical experience with ranibizumab. [START_REF] Silva | Short-term efficacy of intravitreal ranibizumab for myopic choroidal neovascularization[END_REF][START_REF] Kostantidinis | Intravitreal ranibizumab (Lucentis) for the treatment of myopic choroidal neovascularization[END_REF][START_REF] Mones | Intravitreal ranibizumab for choroidal neovascularization secondary to pathologic myopia: 12-month results[END_REF][START_REF] Lai | Intravitreal ranibizumab for the primary treatment of choroidal neovascularization secondary to pathologic myopia[END_REF] Explanations could be that PDT has been stated as in-label treatment and the different costs between anti-VEGF drugs.

Silva et coll. [START_REF] Silva | Short-term efficacy of intravitreal ranibizumab for myopic choroidal neovascularization[END_REF] retrospectively evaluated the efficacy of ranibizumab in myopic CNV in 26 eyes followed up for three -six months, while Kostantidinis et al. [START_REF] Kostantidinis | Intravitreal ranibizumab (Lucentis) for the treatment of myopic choroidal neovascularization[END_REF] prospectively enrolled 14 eyes and followed them for up to an average of 8 months (range 3-16 months). They showed a positive trend in vision recovery, stable vision after the third month and a slower reduction in OCT macular thickness.

Monés [START_REF] Mones | Intravitreal ranibizumab for choroidal neovascularization secondary to pathologic myopia: 12-month results[END_REF] prospectively enrolled 23 patients with CNV related to PM, treated "on demand". Almost 70% of patients gained al least 1 line of visual acuity in a one-year follow-up, with a mean number of injections equal to 1.5. Intravitreal ranibizumab has also been suggested as first line treatment in 16 patients by Lai et al. [START_REF] Lai | Intravitreal ranibizumab for the primary treatment of choroidal neovascularization secondary to pathologic myopia[END_REF], who chose the loading dose as protocol: they observed a gain in vision in 75% of patients and only one patient needed retreatment during the follow-up of 12 months.

In our paper we confirm the efficacy of intravitreal injection in myopic CNV. The final mean VA improved in 32 out of the 40 eyes considered and was stable in five (92.5%). The median number of lines gained was 2.9 from the baseline. About 60 % of all patients improved by 3 or more lines. Ranibizumab stabilized BCVA to the range of legal blindness in one eye only, while it allowed six eyes to recover to mild low vision level (range 1-0.5 LogMAR). Even in those cases where BCVA worsened, it stabilized between 0.6 and 0.2 LogMAR.

We did not notice any difference in visual outcome between treated eyes previously treated with PDT and naïve eyes. Mechanisms of action of PDT and ranibizumab are different, and therefore where PDT is unsuccessful, anti-VEGF blocking is always a treatment opportunity. The role of PDT in promoting fibrosis should however be considered [START_REF] Ruiz-Moreno | Subretinal fibrosis after photodynamic therapy in subfoveal choroidal neovascularization in highly myopic eyes[END_REF] in evaluating visual prognosis.

Which is the most effective regimen to adopt in non-AMD CNV still remains uncertain.

Clinically myopic CNV differs from AMD-CNV, the former being slower in natural evolution with a variable visual outcome. [START_REF] Soubrane | Choroiodal Neovascularization in pathologic myopia: recent developments in diagnosis and treatment[END_REF] The role of pigment epithelium (RPE) and the release of different VEGF-related factors, is probably the key to the differences. [START_REF] Grossniklaus | Clinicopathologic correlation of surgically excised type 1 and type 2 submacular choroidal neovascular membranes[END_REF] Patients are often young and therapies can be more effective because of the health RPE that is able to react adequately, creating a barrier to the CNV development. The same results were observed after PDT treatments in the past. [START_REF] Pece | Photodynamic Therapy With Verteporfin for subfoveal Choroidal Neovascularizzation secondary to Pathologic Myopia: A Longterm Study[END_REF][START_REF] Pece | Photodynamic Therapy With Verteporfin for Juxtafoveal Choroidal Neovascularizzation in Pathologic Myopia: A Long-term Follow-up Study[END_REF][START_REF] Ergun | Prognostic factors influencig visual outcome of photodynamic therapy for subfoveal choroidal neovascularization in pathologic myopia[END_REF] For these reasons we think that, in myopic CNV, the "on demand" regimen may be preferable with respect to the loading phase; this would lead to a lower number of injections, lower risks and less expenses. By now we don't have any controlled trial assessing the best modality of treatment in myopic CNV, however, among the papers published, some indicate very good results with the on demand regimen, the final mean number of injection being very low. In our paper we observed a mean number of injections of 2.8 in the first year, indicating a better prognosis in myopic CNV than in AMD (5.6 reported in PrONTO study). [START_REF] Fung | An optical coherence tomographyguided, variable dosing regimen with intravitreal ranibizumab (Lucentis) for neovascular age-related macular degeneration[END_REF] In our series, the number of treatments were not related to poorer final visual acuity, in contrast with other reports about bevacizumab. [START_REF] Gharbiya | Intravitreal bevacizumab treatment for choroidal neovascularization in pathologic myopia: 12-month results[END_REF] reduction after treatment (p 0.02), the role of OCT needs to be clarified.

In myopic patients OCT may be considered useful to better define the diagnosis of CNV, excluding the presence of vitreoschisis with foveal alterations, tractions, and small macular holes. However, since PM-CNV can often have minimal leakage and choroid is usually very thin, OCT macular thickness can be irrelevant and the presence of fluid very scarce (see fig. 2 case 1). Furthermore the Stratus OCT may give us a lot of artifacts and would therefore need a very skilled examiner.

New spectral domain OCTs will give us more and better clinical information, but at the present moment, considering also that Stratus OCT is the gold standard in current scientific reports, in our opinion we should be very careful in carrying out a treatment or a retreatment of patients on the basis of OCT only. In our experience metamorphopsia is the leading symptom and should be highly considered in reaching a decision.

Several reports underline that refractive error, lesion dimension and age can influence treatment response, both for PDT and for anti-VEGF drugs. [START_REF] Pece | Photodynamic Therapy With Verteporfin for subfoveal Choroidal Neovascularizzation secondary to Pathologic Myopia: A Longterm Study[END_REF][START_REF] Gharbiya | Intravitreal bevacizumab treatment for choroidal neovascularization in pathologic myopia: 12-month results[END_REF][START_REF] Ergun | Prognostic factors influencig visual outcome of photodynamic therapy for subfoveal choroidal neovascularization in pathologic myopia[END_REF] However in our cohort of patients, statistical analysis did not show any correlation with age, refractive error or previous PDT. This issue surely deserves more investigation.

The safety of the drug was confirmed by our report. None of the patients developed any ocular or systemic side effects associated with intravitreal ranibizumab.

Efficacy and safety of a new drug should be tested on a large population as it happened in AMD trials; this seldom happened for myopic population. When evaluating literature about myopic CNV and anti-VEGF drugs, our series is large and can be accepted as representative of this population. We hope larger multicenter trial will offer stronger support to clinical practice. Further limitation of our paper is follow-up extension: we agree that clinical trial or observational studies should have at least 12 months of follow-up or 24-36 months in order to really compare results with natural history. The study lacks of a control group: however, nowdays we consider ethically unacceptable to perform comparison with randomized untreated patients, since therapies really changed prognosis of the disease.

CONCLUSION

International trials with ranibizumab negated any subgrouping of CNV and therefore treatment with anti-VEGF drugs is equally suggested for any form or localization of CNV in AMD. Anti-VEGF drugs give a high chance of stabilizing and improving vision that has otherwise deteriorated because of CNV. These drugs are therefore currently used in any form of CNV, also in off-label settings. [START_REF] Gunther | Bevacizumab (Avastin) for the treatment of ocular disease[END_REF] In our experience, intravitreal injection of ranibizumab appears to be very effective in the treatment of myopic CNV, with a high proportion of patients reporting an improvement in vision and a mean of 3 injections in the first year. Treatment is also able to avoid legal blindness. The correct regimen of treatment, as well as retinal and choroidal changes with the time, still need to be defined

  The primary outcome was considered to be the improvement of visual acuity with a minimum of 5 letters (one ETDRS line), or the stabilization of VA. Secondary outcome was the changes in OCT central retinal thickness. Alfa level (p) was fixed at 0.05, as a probability of 95%. Statistical analysis was performed by means of a non-parametric Wilcoxon signed-ranks test and two-tailed t test, provided by SAS (Statistical Analysis Software) for Windows 8.0 (Microsoft Corp. USA), in order to evaluate BCVA and OCT central retinal thickness changes. The influence of age, previous PDT and lesion size on VA changes was analyzed by means of multivariate linear regression.

Figure 1 :

 1 Figure 1: results are presented as scatterplot of basal and final BCVA. Scatterplot presents

Figure 2 :

 2 Figure 2: demonstrative fluorescein angiography and OCT of three patients (indicated as

Table 1 :

 1 baseline characteristics of patients at the study entry

	*central retinal thickness
	measured by optical
	coherence tomography.

Table 2 :

 2 outcomes at the end of the follow-up

	Findings	Values	p values §
	Mean LogMAR BCVA ± SD	0.27 (±0.2)	0.008
	Mean letter BCVA ± SD	40.5 (±14)	0.01
	Improvement (gain ≥ 1 line) (n.	32/40 (80)	=1 line 9/32 (28)
	eyes (%))		= 2 lines 4/32 (12.5)
			≥ 3 lines 19/32 (59.5)
	Stability (n. eyes (%))	5/40 (12.5)
	Worsening (loss ≤ 1 line)	3/40 (7.5)
	(n. eyes (%))
	Mean OCT CRT (± SD) (microns)	175 (±46)	0.02
	Mean n. injections (±SD)	2.8 (±1.2)
	§ p values when comparing with values at the baseline

  Also the mean number of injections (2.56 vs 2.6) was unremarkably different. The mean number of letters gained was 31 vs 34 (p 0.2).

No statistical accordance was found (p> 0.05) regarding age, refractive errors, lesion size, BCVA, central macular thickness and the response to a ranibizumab injection.

ACKNOWLEDGEMENTS

To Giuseppe Bianco, clever help in managing surgical times.

Licence for Publication:

The Corresponding Author has the right to grant on behalf of all authors and does grant on behalf of all authors, an exclusive licence (or non exclusive for government employees) on a worldwide basis to the BMJ Publishing Group Ltd and its Licensees to permit this article (if accepted) to be published in BJO editions and any other BMJPGL products to exploit all subsidiary rights, as set out in our licence (http://group.bmj.com/products/journals/instructions-forauthors/licence-forms/).

Competing Interest:

None to declare