Adnane Hamiaz

Michel Mehrenberger

Hocine Sellama

Eric Sonnendrücker

The semi-Lagrangian method on curvilinear grids

HAL is

Introduction

Plasmas, which are a collection of charged particles, can be described quite accurately by kinetic models like the Vlasov-Maxwell equations, or in some circumstances reduced models like the Vlasov-Poisson equation if low frequency phenomena are of interest or the gyrokinetic model when a strong background magnetic field is present. These models nonlinearly couple the Vlasov equation, which is a transport equation in phase space, with the Maxwell equations, which describe the evolution of the self-consistent electromagnetic field generated by the charged particles. The coupling is performed by solving alternatively the Vlasov and the Maxwell equations. In this paper, we are specifically interested by the Vlasov part and will not discuss the coupling, or the field equations. All these Vlasov or Vlasov type equations in standard cartesian coordinates in the d-dimensional phase space R d can be written equivalently in the two following abstract forms

∂ t f + ∇ x • (af ) = 0, (1.1)
or

∂ t f + a • ∇ x f = 0, (1.2)
where a : Ω ⊂ R d → R is a divergence free vector field, which means that it satisfies ∇ x • a = 0, and the solution f = f (t, x), t ∈ R + , x ∈ R d . Equation (1.1) is called the conservative form and is equivalent, as the vector field a is divergence free, to (1.2), which is called the advective form.

Because of the high dimension of phase space, up to six for the full physical problem, the Monte Carlo type particle in cell method (PIC) is one of the most used numerical method for its solution, see the book [START_REF] Birdsall | Plasma physics via computer simulation[END_REF] or the review article [START_REF] Verboncoeur | Particle simulation of plasmas: review and advances[END_REF] for more information on this method. Another method, which has been very successful for plasma simulations is the semi-Lagrangian method, which couples particle tracking with a projection on a phase space grid at each time step. The idea in the simple setting of the Vlasov-Poisson equations dates back to Cheng and Knorr [START_REF] Cheng | The integration of the vlasov equation in configuration space[END_REF] and to [START_REF] Sonnendrücker | The semilagrangian method for the numerical resolution of the vlasov equation[END_REF] in the general setting we are interested in here, based on the semi-Lagrangian method, which was already quite successful for climate applications [START_REF] Staniforth | Semi-lagrangian integration schemes for atmospheric models-a review[END_REF]. The big advantage of the semi-Lagrangian method, with respect to other grid based methods is that there is no strong CFL condition on the time step, which for many problems in plasma physics introduces a time step restriction, which is too severe due to the large velocity grid points compared to what is needed for physical accuracy. Since then, many variants of the semi-Lagrangian method have been introduced for plasma physics applications. A convergence proof of this method appears in [START_REF] Besse | Convergence of a semi-lagrangian scheme for the one-dimensional vlasov-poisson system[END_REF] for linear interpolation and in [START_REF] Besse | Convergence of classes of high-order semilagrangian schemes for the vlasov-poisson system[END_REF] for high order interpolation. The estimates have been improved recently in [START_REF] Charles | Enhanced convergence estimates for semi-lagrangian schemes application to the vlasov-poisson equation[END_REF]. A semi-Lagrangian method on unstructured meshes has been proposed in [START_REF] Besse | Semi-lagrangian schemes for the vlasov equation on an unstructured mesh of phase space[END_REF]. These methods, which solve the advective form (1.2) of the Vlasov equation, go in the standard class of semi-Lagrangian methods where characteristics are advected backward in time. We call this classical semi-Lagrangian method BSL for backward semi-Lagrangian method. This is opposed in particular to CSL, conservative semi-Lagrangian methods, which solve the conservative form (1.1) of the Vlasov equation still with backward characteristics as in [START_REF] Nakamura | Exactly conservative semi-lagrangian scheme for multi-dimensional hyperbolic equations with directional splitting technique[END_REF][START_REF] Crouseilles | Conservative semi-lagrangian schemes for vlasov equations[END_REF][START_REF] Umeda | A non-oscillatory and conservative semi-lagrangian scheme with fourth-degree polynomial interpolation for solving the vlasov equation[END_REF], and FSL, forward semi-Lagrangian methods introduced in [START_REF] Crouseilles | A forward semilagrangian method for the numerical solution of the vlasov equation[END_REF] for the Vlasov equation, which solves the characteristics forward in time. There have also been works on positivity preserving conservative Discontinuous Galerkin methods [START_REF] Qiu | Positivity preserving semi-lagrangian discontinuous galerkin formulation: theoretical analysis and application to the vlasov-poisson system[END_REF][START_REF] Rossmanith | A positivity-preserving high-order semilagrangian discontinuous galerkin scheme for the vlasov-poisson equations[END_REF] and also on semi-Lagrangian methods on adaptive grids based on wavelet interpolation [START_REF] Gutnic | Vlasov simulations on an adaptive phase-space grid[END_REF][START_REF] Besse | A wavelet-mra-based adaptive semi-lagrangian method for the relativistic vlasov-maxwell system[END_REF], see also [START_REF] Pinto | Convergence of an adaptive semilagrangian scheme for the vlasov-poisson system[END_REF] for a convergence proof. Note also a hybrid method using only a semi-Lagrangian method in the velocity space [START_REF] Guo | Hybrid semi-lagrangian finite element-finite difference methods for the vlasov equation[END_REF], and a conservative semi-Lagrangian method based on WENO reconstruction [START_REF] Qiu | Conservative semi-lagrangian finite difference weno formulations with applications to the vlasov equation[END_REF][START_REF] Qiu | Conservative high order semi-lagrangian finite difference weno methods for advection in incompressible flow[END_REF]. Now somewhat outdated comparisons of different types of methods for the Vlasov equation can be found in [START_REF] Arber | A critical comparison of eulerian-grid-based vlasov solvers[END_REF][START_REF] Filbet | Comparison of eulerian vlasov solvers[END_REF].

The present paper is dedicated to a difficulty appearing in applications of the semi-Lagrangian method to the gyrokinetic model used in turbulence simulations of magnetic fusion plasmas [START_REF] Grandgirard | A driftkinetic semi-lagrangian 4d code for ion turbulence simulation[END_REF][START_REF] Grandgirard | Computing itg turbulence with a full-f semi-lagrangian code[END_REF][START_REF] Latu | Gyrokinetic semi-lagrangian parallel simulation using a hybrid openmp/mpi programming[END_REF]. Due to the geometry of magnetic fusion devices and the strong magnetic field which creates a strong anisotropy on the solution and therefore imposes, if a reasonable cost is to be kept, that the mesh be aligned on magnetic flux surfaces, it is standard to use specific curvilin-ear grids, aligned on magnetic field lines or flux surfaces, for these simulations. It is well known that such grids introduce additional numerical difficulties, in particular the free stream preservation issue for conservative methods [START_REF] Colella | High-order, finite-volume methods in mapped coordinates[END_REF][START_REF] Persson | Discontinuous galerkin solution of the navier-stokes equations on deformable domains[END_REF]. The free stream preservation enforces that constant states are exactly preserved by the numerical method. This is enforced trivially by the classical BSL method independently of the accuracy of the advection as constants are preserved by interpolation. For this reason it appears that they are more robust in curvilinear coordinates. On the other hand, a first attempt of using a split conservative semi-Lagrangian method has been reported in [START_REF] Braeunig | Some numerical aspects of the conservative psm scheme in a 4d drift-kinetic code[END_REF], pinpointing some specific difficulties.

An important property of the advection field for the method to be stable in curvilinear coordinates is that it is divergence free at the discrete level in some sense. This is most easily realised, when using a potential formulation of the advection field. This is also most easily understood in two-dimensions, but can be generalised to arbitrary dimensions as shown in [START_REF] Colella | High-order, finite-volume methods in mapped coordinates[END_REF]. So in order to make our point it is enough to study the problem in a two-dimensional space. We suppose that we can write the advection field a

= (∂ x2 Φ, -∂ x1 Φ) = ∇ ⊥ Φ,
where stands for the transposition, and

∇ ⊥ = (∂ x2 , -(∂ x1 )) .
In applications a may depend on t and f (t, •), but we will omit this point for the moment.

The classical BSL methods automatically preserve constant states, but are not conservative. On the other hand the conservative CSL method is conservative, but does not preserve constant states. We will show in this paper that the BSL method can be made conservative up to at least first order in time and on the other some CSL methods can be made to preserve constant states up to at least first order in time. Moreover different interpolating methods will be compared for BSL to show its robustness. Some numerical investigations of the splitting issue raised in [START_REF] Huot | Instability of the time splitting scheme for the one-dimensional and relativistic vlasov-maxwell system[END_REF] will also be performed.

The outline of the paper is the following: First we are going to express the two abstract forms of the Vlasov equations in curvilinear coordinates. In the second part we recall the principles of the classical and conservative semi-Lagrangian methods. Then we are going to propose different versions of the BSL scheme and a new CSL scheme and prove their approximate conservation properties. Finally detailed numerical investigations will be performed to assess the strengths and weaknesses of the different schemes.

Obtention of the curvilinear equations

We consider here the curvilinear framework, that is we have a transformation (the mapping) F : Ω → Ω; η → x(η) from a logical domain Ω to the physical domain Ω ⊂ R 2 . Polar coordinates are a simple example of such a transformation, but more general, numerically defined transformations, are needed to align the grid with the magnetic flux surfaces in actual tokamak simulations. We want to define the equation for f (η) = f (x(η)). In the sequel we will omit the tilde on f and other quantities, as on each identity it should be clear form the context if we deal with functions of x (that is, from the physical domain) or functions of η (from the logical domain). We define the jacobian matrix [J] and the jacobian J by

[J] = ∂ η 1 x 1 ∂ η 2 x 1 ∂ η 1 x 2 ∂ η 2 x 2 , J = det([J]).
Using the chain rule we can relate the differential operators with respect to x to the differential operators with respect to η. Denoting by [J] -= ([J] -1 ) , we get

∇ x f = [J] -∇ η f, ∇ x • (af ) = 1 J ∇ η • ([J] -1 aJf ), (2.1)
so that the advective form reads

∂ t f + ([J] -1 a) • ∇ η f = 0,
and the conservative form becomes

∂ t (Jf ) + ∇ η • ([J] -1 a Jf ) = 0
and the divergence free condition is

∇ η • ([J] -1 a J) = 0.
Again from the chain rule, we find that

∇ ⊥ η Φ = J[J] -1 ∇ ⊥ x Φ.
So that if a = ∇ ⊥ x Φ, we get J[J] -1 a = ∇ ⊥ η Φ. The advective form then rewrites

∂ t f + 1 J ∇ ⊥ η Φ • ∇ η f = 0, (2.2)
or expanding the coordinates

∂ t f + 1 J ∂Φ ∂η 2 ∂f ∂η 1 - 1 J ∂Φ ∂η 1 ∂f ∂η 2 = 0,
and the conservative form rewrites

∂ t (Jf ) + ∇ η • ((∇ ⊥ η Φ) f ) = 0, (2.3)
which becomes, expanding the coordinates

∂ t (Jf ) + ∂ ∂η 1 ( 1 J ∂Φ ∂η 2 Jf ) - ∂ ∂η 2 ( 1 J ∂Φ ∂η 1 Jf ) = 0,
and the divergence free condition is automatically satisfied, as

1 J ∇ η • ([J] -1 a J) = 1 J ∇ η • ∇ ⊥ η Φ = 1 J ∂ ∂η 1 ( ∂Φ ∂η 2 ) - 1 J ∂ ∂η 2 ( ∂Φ ∂η 1 ) = 0.

The characteristics

The semi-Lagrangian method is based on conservation properties along the characteristics. Let us recall our Vlasov type equations in curvilinear coordinates that we are interested in, denoting by b = ∇ ⊥ η φ the divergence free advection field and J the jacobian of our curvilinear grid. The advective form (2.2) then writes

∂ t f + 1 J b • ∇ η f = 0, (3.1)
and introducing f = Jf the conserved variable, the conservative form (2.3) becomes

∂ t f + ∇ η • ( 1 J b f ) = 0. (3.2)
To both of these equations we associate the same characteristics, which are the integral curves of the differential equation

dH dt = 1 J b, H(s) = η. (3.3)
We shall denote classically by H(t; η, s), the solution at time t of the characteristic curve taking the value η at time s.

The solutions f of the advective form (3.1) and f = f J of the conservative form (3.2) obey the following conservation properties along the characteristics: Lemma 3.1 For given smooth vector field b and function J, the solution f of (3.1) is conserved along the characteristics. More precisely

d
+ f (t, H(t; η, s))(∇ η • b J )(t, H(t; η, s)) det(∇ η H(t; η, s)) = ∂ f ∂t (t, H(t; η, s)) + ∇ η • ( b J f ) (t, H(t; η, s)) det(∇ η H(t; η, s)) = 0,
as f is a solution of the conservation law at any point.

4

The Backward Semi-Lagrangian scheme (BSL2D)

The backward semi-Lagrangian scheme computes the time evolution of our equation in advective form (3.1) using the characteristics (3.3).

Formulation of the scheme

We consider a uniform cartesian grid of the logical space (η

1 , η 2 ) ∈ [0, 1] 2 defined by N 1 , N 2 ∈ N * and η 1 i = (i -1)∆η 1 , η 2 i = (i -1)∆η 2 , i ∈ R, using a Fortran like indexing, with ∆η k = 1/N k , k = 1, 2. We also define t n = n∆t, n ∈ R, with ∆t ∈ R + * .
The unknowns at time t n are

f n i,j f (t n , η 1 i , η 2 j ), i = 1, . . . , N 1 + 1, j = 1, . . . , N 2 + 1, n ∈ N.
The classical backward semi-Lagrangian scheme uses the conservation property (3.4) at each grid point to update the solution from one time step to the next. This can be expressed as

f n+1 i,j = f (t n , η * i,j , t n+1 )),
where

η i,j = (η 1 i , η 2 j
) is a grid point and η * i,j = H(t n ; η i,j , t n+1 ) the origin of the characteristic ending at that grid point. In order to get a numerical scheme from this equality, we need two steps:

1. Compute H(t n ; η i,j , t n+1 )) by numerically solving the characteristics equations (3.3) backward in time on one time step. If the advection field b/J is analytically known, any standard ODE solver can be used. When the problem is non linear and the advection field depends on f , a simple solution is to use the backward Euler method. Note however that iterations are needed as the advection field is not known at time t n+1 . Two are generally sufficient. The scheme reads

η * i,j = η i,j -∆t(b/J)(t n+1 , η i,j ), i = 1, . . . , N 1 + 1, j = 1, . . . , N 2 + 1.
2. Interpolate f n the solution at t n at the origin of the characteristics. From grid values f h = (f i,j ) i=1,...,N1+1, j=1,...,N2+1 , we need to define a reconstruction

Π[f h ] : R 2 → R.
This reconstruction (the interpolator ) should satisfy

Π[f h ](η 1 i , η 2 j ) = f i,j , i = 1, . . . , N 1 + 1, j = 1, . . . , N 2 + 1.
The algorithm then becomes:

1. Initialisation:

f 0 i,j = f (t 0 , η 1 i , η 2 j ), i = 1, . . . , N 1 + 1, j = 1, . . . , N 2 + 1,
2. From t n to t n+1 . Predict advection field b at time t n+1 for non linear problems and loop until convergence (a) Compute origin of characteristics

η * i,j = η i,j -∆t(b/J)(t n+1 , η i,j ), i = 1, . . . , N 1 + 1, j = 1, . . . , N 2 + 1. (b) Interpolate f n+1 i,j = Π[f n h ] η 1 i -∆ta 1 i,j , η 2 j -∆ta 2 i,j , i = 1, . . . , N 1 + 1, j = 1, . . . , N 2 + 1 with f n h = (f n i,j ) i=1,...,N1+1, j=1,...,N2+1 .

Preservation of constant states

As soon as the reconstruction Π reproduces the constants, we have preservation of constant states, which can be written in this form, as the operator is linear:

f n i,j = 1, i = 1, . . . , N 1 + 1, j = 1, . . . , N 2 + 1 ⇒ f n+1 i,j = 1, i = 1, . . . , N 1 + 1, j = 1, . . . , N 2 + 1 .

First order conservation of mass in BSL2D

A priori, conservation of mass is not satisfied for the BSL2D scheme. Using the jacobian at the grid points J i,j = J(η 1 i , η 2 j ) , this would read i,j

J i,j f n+1 i,j = i,j J i,j f n i,j .
However, assuming our reconstruction is smooth enough, we can write a Taylor expansion of first order in ∆t, from our scheme:

f n+1 i,j = Π[f n h ] η 1 i -∆tb 1 i,j /J i,j , η 2 j -∆tb 2 i,j /J i,j = Π[f n h ](η 1 i , η 2 j ) - ∆t J i,j b 1 i,j ∂Π[f h ] ∂η 1 (η 1 i , η 2 j ) + b 2 i,j ∂Π[f h ] ∂η 2 (η 1 i , η 2 j ) + O(∆t 2 ). So that i,j J i,j f n+1 i,j = i,j J i,j f n i,j -∆t b 1 i,j ∂Π[f h ] ∂η 1 (η 1 i , η 2 j ) + b 2 i,j ∂Π[f h ] ∂η 2 (η 1 i , η 2 j ) + O(∆t 2 ).
So a first order condition is the following

A = N1 i=1 N2 j=1 b 1 i,j ∂Π[f h ] ∂η 1 (η 1 i , η 2 j ) + b 2 i,j ∂Π[f h ] ∂η 2 (η 1 i , η 2 j ) = 0,
which should be valid for all periodic sequences f h . Now, considering that

b 1 i,j = ∂Φ h ∂η 2 (η 1 i , η 2 j ), b 2 i,j = - ∂Φ h ∂η 1 (η 1 i , η 2 j ),
with a reconstruction Φ h Φ, the first order condition rewrites

A = N1 i=1 N2 j=1 ∂Φ h ∂η 2 ∂Π[f h ] ∂η 1 - ∂Φ h ∂η 1 ∂Π[f h ] ∂η 2 (η 1 i , η 2 j ) = 0. (4.1)
Using a reconstruction of the derivatives of the form

∂F ∂η 1 (η 1 i , η 2 j ) = s1 =r1 α 1 F (η 1 i+ , η 2 j ), ∂F ∂η 2 (η 1 i , η 2 j ) = s2 =r2 α 2 F (η 1 i , η 2 j+ ), (4.2)
for both F = Π[f h ] and F = Φ h leads to

A = N1 i=1 N2 j=1 s1 k=r1 s2 =r2 α 1 k α 2 (Π[f h ](η 1 i+k , η 2 j )Φ h (η 1 i , η 2 j+ ) -Φ h (η 1 i+k , η 2 j )Π[f h ](η 1 i , η 2 j+ )) = N1 i=1 N2 j=1 Π[f h ](η 1 i , η 2 j ) s1 k=r1 s2 =r2 α 1 k α 2 (Φ h (η 1 i-k , η 2 j+ ) -Φ h (η 1 i+k , η 2 j-)) = N1 i=1 N2 j=1 Π[f h ](η 1 i , η 2 j ) k (α 1 i-k α 2 -j -α 1 k-i α 2 j-)Φ h (η 1 k , η 2 ),
prolongating α 1 i = 0, i ∈ {r 1 , . . . , s 1 } and α 2 j = 0, j ∈ {r 2 , . . . , s 2 }. The condition A = 0 is thus true as soon as

α k -i = -α k i , i ∈ Z, k = 1, 2
, which is a reasonable assumption for the derivative formula (4.2). We emphasize that such condition holds when the reconstructions Π[f h ] and Φ h use the same formula for the derivative. We also remark that the condition (4.1) is the same for cartesian and curvilinear case, as the jacobian is simplified and does no more appear. For splines, the derivative is generally not expressed like in (4.2): derivatives are obtained, by inverting a linear system. However, in the periodic setting, the matrix for computing the derivatives is circulant and its inverse also, which leads to a formula of type (4.2).

Examples of reconstructions: the 1D advective case

In this section, we review some choices for Π h . We focus here on the 1D interpolation in the advective form, as the 2D case can be obtained through tensor product.

Spline interpolation: SPL(d) (see e.g. [START_REF] Sonnendrücker | The semilagrangian method for the numerical resolution of the vlasov equation[END_REF] (for d = 3), [START_REF] Besse | Convergence of classes of high-order semilagrangian schemes for the vlasov-poisson system[END_REF], [START_REF] Afeyan | Simulations of kinetic electrostatic electron nonlinear (keen) waves with two-grid, variable velocity resolution and high-order time-splitting[END_REF] (for d = 3 and non-uniform grid)) We can consider here a non uniform mesh. Let η min < η max ∈ R. We consider a discretization

η 1 = η min ≤ η 2 ≤ • • • ≤ η N +1 = η max of the interval [η min , η max ]. Let d ∈ N *
, an odd number. We define a knot sequence τ i , i = 1, . . . , N + 2d + 1. The knot sequence satisfies

τ j+d = η j , j = 1, . . . , N + 1.
In the case of open boundary conditions, we take

τ j = τ d+1 , j = 1, . . . , d, τ N +j = τ N +d+1 , j = d + 2, . . . , 2d + 1,
and for periodic boundary conditions, we take

τ j = τ j+N -L, j = 1, . . . , d, τ N +j = τ j + L, j = d + 2, . . . , 2d + 1,
where L = η max -η min . We can consider other possibilities for the 2d knots τ 1 , . . . , τ d and τ N +d+2 , . . . , τ N +2d+1 ; we consider that

τ 1 ≤ • • • ≤ τ d ≤ τ d+1 , τ N +d+1 ≤ τ N +d+2 ≤ • • • ≤ τ N +2d+1 .
We then can define, for j = 1, . .

b 1 ← 1 for = 1, . . . , d do α ← x-τ m+1- τm+1-τ m+1-b 1 b 1 ← b 1 -α for k = 2, . . . , do β = x-τ m+k- τ m+k -τ m+k-b k b k ← b k + α -β α ← β end for b +1 ← α end for
Now, we look for a function f h of the form

f h (η) = N +d j=1 c j B d+1 j (η), η min ≤ η < η max .
For i = 1, . . . , N , we have

f i = f h (η i ) = d k=1 b k (τ i+d )c i-1+k .
This is enough for periodic boundary conditions. For open boundary conditions, we can add the d conditions

f ( ) 1 = f ( ) h (η 1 ) = d k=1 b ( ) k (τ d+1 )c k , = 1, . . . , (d -1)/2, f ( ) N +1 = f ( ) h (η N +1 ) = d+1 k=2 b ( ) k (τ - N +d+1 )c N -1+k , = 0, . . . , (d -1)/2.
Typically, when periodic boundary conditions are not used, we take 0 for the derivatives and extend the solution outside [η min , η max [ to be

f h (η) = f 1 , η ≤ η min , f h (η) = f N +1 , η ≥ η min .
Note that, from [START_REF] Volkov | On complete interpolation spline finding via b-splines, Sibirskie Èlektronnye Matematicheskie Izvestiya[END_REF], we can then obtain explicitly the d+1 coefficients c k , c N -1+k , k = 1, . . . , (d + 1)/2, and the N -1 remaining coefficients are then solution of a ddiagonal system, with matrix

M =       b (d+1)/2 (τ 2+d ) . . . b d (τ 2+d ) 0 . . . . . . 0 0 . . . 0 b 1 (τ i+d ) . . . b (d+1)/2 (τ i+d ) . . . b d (τ i+d ) 0 . . . 0 0 . . . . . . 0 b 1 (τ N +d ) . . . b (d+1)/2 (τ N +d )      
.

Once the coefficients are computed, we obtain for i = 1, . . . , N, 0 ≤ α < 1

f h (η) = d+1 k=1 c i-1+k b k (τ i+d + α(τ i+d+1 -τ i+d )), η = η i + α(η i+1 -η i ).
Lagrange interpolation: LAG(2d+1) (see e.g. [START_REF] Filbet | Comparison of eulerian vlasov solvers[END_REF][START_REF] Charles | Enhanced convergence estimates for semi-lagrangian schemes application to the vlasov-poisson equation[END_REF]33])

We consider here a uniform mesh. Lagrange interpolation of degree 2d + 1 reads

f h (η) = d+1 k=-d w d k (α)f i+k , η = η i + α(η i+1 -η i ), i = 1, . . . , N, 0 ≤ α < 1, with w d k (α) = d+1 j=-d, j =k (α-j)
d+1 j=-d, j =k (k-j) . As boundary conditions, apart from periodic, we often use

f j = f N +1 , j ∈ N, j ≥ N + 1; f j = f 1 , j ∈ Z, j ≤ 1. (4.3)
Cubic Hermite interpolation: H(p) (see e.g. [START_REF] Filbet | Comparison of eulerian vlasov solvers[END_REF] (for p = 4), [START_REF] Steiner | Numerical computation of the gyroaverage operator and coupling with the Vlasov gyrokinetic equations[END_REF]) Let p ∈ N. We consider here again a uniform mesh, and write

f h (η) = w 0,0 (α)f i + w 0,1 (α)f i + + w 1,0 (α)f i+1 + w 1,1 (α)f (i+1) -, η = η i + α(η i+1 -η i ), i = 1, . . . , N, 0 ≤ α < 1, with w 1,0 (α) = α 2 (3-2α), w 0,0 (α) = 1-w 1,0 (α), w 0,1 (α) = α(1-α) 2 , w 1,1 (α) = α 2 (α -1).
The derivatives are reconstructed in the following way

f i ± = s ± =r ± b ± f i+ .
For p ∈ N * , we take r

+ = -p 2 , s + = p+1 2 , r -= -s + , s -= -r + , -b - -= b + = s +
j=r + , j ∈{0, } (-j)

s + j=r + , j = ( -j) , = r + , . . . , s + , = 0, -b - 0 = b + 0 = - s + j=r + , j =0 b + j .
As examples, we have (negative indices are before ; in the brackets) Boundary conditions are treated like for the Lagrange case.

p = 1 b + = (-1, 1), p = 2 b + = (-1/2; 0, 1/2), p = 3 b + = (-1/3; -1/2, 1, -1/6), p = 4 b + = (1/

Extension to 2D for BSL2D

This is achieved using a tensor product formulation. For the splines, for each j = 1, . . . , N 2 + 1, we compute first coefficients c 1 ij , i = 1, . . . , N 1 + d 1 , from the values f ij , i = 1, . . . , N 1 + 1. Then, for each i = 1, . . . , N 1 + d 1 , we compute coefficients c ij , j = 1, . . . , N 2 + d 2 , from the values c 1 ij , i = 1, . . . , N 1 + d 1 , j = 1, . . . , N 1 + 1. The formula for interpolation is then

f h (η 1 , η 2 ) = d1+1 k=1 d2+1 =1 c i-1+k,j-1+ b k (τ i+d1 + α(τ i+d1+1 -τ i+d1 )) b (τ j+d + β(τ j+d2+1 -τ j+d2 )), for η 1 i = η 1 i + α(η 1 i+1 -η 1 i ), i = 1, . . . , N 1 , 0 ≤ α < 1, η 2 j = η 2 j + β(η 2 j+1 -η 2 j ), j = 1, . . . , N 2 , 0 ≤ β < 1.
Note that for open boundary conditions (in both directions), in addition to

f ij , i = 1, . . . , N 1 + 1, j = 1, . . . , N 2 + 1, we have to specify f ( ,0) ij = ∂ f h ∂η 1 (η 1 i , η 2 j ), j = 1, . . . , N 2 + 1, i ∈ {1, N 1 + 1}, = 1, . . . (d -1)/2,
and also

f (0, ) ij = ∂ f h ∂η 2 (η 1 i , η 2 j ), i = 1, . . . , N 1 + 1, j ∈ {1, N 2 + 1}, = 1, . . . (d -1)/2, together with f ( 1 , 2 ) ij = ∂ 1 ∂ 2 f h ∂η 1 ∂η 2 (η 1 , η 2 ), η k ∈ {η k 1 , η k N k +1 }, k = 1, . . . (d k -1)/2, k = 1, 2.
From f

( ,0) ij and f ij , we compute c 1 ij , i = 1, . . . , N 1 + d, j = 1, . . . , N 2 + 1. From f ( 1 , 2) ij and f (0, ) ij , we compute missing coefficients called c 1,( ) ij , i = 1, . . . , N 1 + d, j ∈ {1, N 2 + 1}, = 1, . . . , (d 2 -1)/2,
in order to compute finally, together with c 1 ij the coefficients c ij . For Lagrange interpolation, the formula is

f h (η 1 , η 2 ) = d1+1 k=1 d2+1 =1 w d1 k (α)w d2 (β)f i+k,j+ ,
and for Hermite, we can write

f h (η 1 , η 2 ) = 1 k1,k2, 1, 2=0 w k1,k2 (α)w 1, 2 (β)f (k1, 1 ) i(k2),j( 2) , defining i(0) = i + , i (1)
= (i + 1) -and similarly, j(0) = j + , j(1) = (j + 1) -. We take here f

(0,•) i ± ,• = f (0,•) i,•
and f

(•,0) •,j ± = f (•,0)
•,j . We have to compute the values

f (1,0) i±,j , i = 1, . . . , N 1 + 1, j = 1, . . . , N 2 + 1, f (0,1) i,j± , i = 1, . . . , N 1 + 1, j = 1, . . . , N 2 + 1, f (1,1)
i±,j± , i = 1, . . . , N 1 + 1, j = 1, . . . , N 2 + 1.

The values f

(1,0) i±,j are obtained from f i,j using the operator described in the 1D case in the first direction. The values f (0,1) i,j± are obtained from f i,j using the same operator in the second direction. Finally, the values f (1,1) i±,j± are obtained from f (1,0) i,j± using the operator in the second direction (or equivalently from f (0,1) i±,j using the operator in the first direction).

A split conservative semi-Lagrangian method

On curvilinear grids the free stream preservation property is essential for conservative finite difference, finite volume or finite element solvers. This means that for our problem a constant function must be preserved by a constant coefficient advection.

Although the natural form of the conservative semi-Lagrangian method, as proposed in [START_REF] Crouseilles | Conservative semi-lagrangian schemes for vlasov equations[END_REF], is based on a finite volume formulation where the unknown is a cell average, we find it simpler to use a point based conservative method with a point based unknown as in a Finite Difference scheme. This facilitates a symmetric reconstruction between the stream function Φ of the advection field and the unknown, which for the conservative form of the equation (3.2) is f = f J as we saw before. And, even more importantly this enables to define the Jacobian J only at the grid points, as well for the distribution function as for the computation of the characteristics. This avoids interpolation of the Jacobian, which is known to be bad for free stream preservation. In the classical Finite Volume type method of [START_REF] Crouseilles | Conservative semi-lagrangian schemes for vlasov equations[END_REF], the advected points are the cell edges and the Jacobian needs to be interpolated or defined in another manner there, which prevents a direct cancellation as we have for the point based method used in BSL and proposed here for the conservative scheme. A conservative semi-Lagrangian method of this type was introduced by Qiu and Shu [START_REF] Qiu | Conservative high order semi-lagrangian finite difference weno methods for advection in incompressible flow[END_REF]. We will consider here a specific example, that permits to have the first order free stream preservation; a further study may be pursued for designing a class of high order stable schemes, sharing this property, in the same spirit of the BSL method; one difficulty is to get the stability of the scheme, as already mentioned in [START_REF] Qiu | Conservative high order semi-lagrangian finite difference weno methods for advection in incompressible flow[END_REF]. A general discussion on the respective advantages and drawbacks of Finite Difference, Finite Volume and Discontinuous Galerkin discretisations for conservation laws is given in [START_REF] Shu | High-order finite difference and finite volume weno schemes and discontinuous galerkin methods for cfd[END_REF].

Formulation of the 1D split conservative semi-Lagrangian method (CSL)

Let us start from the model equation in conservative form (3.2) written in logical coordinates, with f = Jf , which we recall

∂ f ∂t + ∇ η • ( 1 J (∇ ⊥ η φ) f ) = 0 (5.1)
with (∇ ⊥ η φ) = (∂ η2 φ, -∂ η1 φ) . Using an operator splitting method, which can be of arbitrary order by alternating the 1D solves with adequate coefficients, we need to solve the following two 1D advections:

∂ f ∂t + ∂ ∂η 1 1 J ∂φ ∂η 2 f = 0, ∂ f ∂t - ∂ ∂η 2 1 J ∂φ ∂η 1 f = 0. (5.2)
This leads us to introducing the conservative algorithm for a 1D conservative problem of the form

∂ f ∂t + ∂ ∂η a J f = 0. ( 5

.3)

A key ingredient in our algorithm is the sliding average reconstruction, which is classical in conservative Finite Difference methods for conservations laws. This reconstruction aims at defining for a function G defined at grid points a high order approximation of its derivative ∂ η G at the same grid points. For this, we introduce a function G related to G by the formula

G(η) = 1 ∆η η+ ∆η 2 η-∆η 2 G(η) dη, so that ∂G ∂η (η i ) = G(η i + ∆η 2 ) -G(η i -∆η 2 ) ∆η . (5.4)
When the averages of G on each cell, which are also the grid values G(η i ), are known, the reconstruction by primitive, standard in Finite Volume methods can be used to approximate G at the cell extremities. Formulas for third and fifth order reconstructions are given in [START_REF] Qiu | Conservative high order semi-lagrangian finite difference weno methods for advection in incompressible flow[END_REF], we shall use here a fourth order reconstruction, for which it is straightforward to compute

G(η i + ∆η 2 ) 7 12 (G(η i ) + G(η i+1 ) - 1 12 (G(η i-1 ) + G(η i+2 )). (5.5)
Now, a conservative semi-Lagrangian algorithm is based on the conservation property (3.5), which becomes for the 1D case

d dt f (t, H(t; η, s)) ∂H ∂η (t; η, s) = 0.
Introducing F a primitive in η of f ,

F (t, η) := η η1 f (t, η) dη,
where η 1 is chosen arbitrarily for fixing the constant, this can be written equivalently ∂ ∂t ∂F * ∂η (t, η, s) = 0, with F * (t, η, s) := F (t, H(t; η, s)).

(5.6)

Let us introduce a uniform grid of our logical domain [0, 1] characterised by the number of cells N . The grid points are then defined by η i = (i -1)∆η, 1 ≤ i ≤ N + 1, with ∆η = 1/N . The unknowns are defined at the grid points f n i f (t n , η i ) and we denote by f n h = ( f n i ) 1≤i≤N +1 the collection of the grid values. Integrating first in t between t n and t n+1 in (5.6) and evaluating at

s = t n+1 , it yields f (t n+1 , η) = ∂F * ∂η (t n , η, t n+1 ),
as H(t n+1 ; η, t n+1 ) = η and ∂H ∂η (t n+1 ; η, t n+1 ) = 1. We suppose for the moment that η * i = H(t n ; η i , t n+1 ), the origin of the characteristic ending at η i satisfies

η i ≤ η * i < η i+1 . (5.7) Defining G(η) = F * (t n , η, t n+1 ) -F (t n , η), we get f (t n+1 , η i ) = f (t n , η i ) + ∂G ∂η (η i ).
In order to get the scheme, we approximate ∂G ∂η (η i ) from values

G(η i ) = η * i ηi f (t, η) dη (5.8)
as explained before, that is using (5.4) together with (5.5).

To complete our algorithm we need to approximate G(η i ). For this, we use Lagrange interpolation of order 3 on the interval [η i , η i+1 [, using the values (5.8). Note that other choices may be possible, but we warn the reader that, as noted in [START_REF] Qiu | Conservative high order semi-lagrangian finite difference weno methods for advection in incompressible flow[END_REF], for stability reasons, such interpolation has to be combined carefully with the quadrature formula which is here (5.5).

f n i+ , = -1, . . . , 2, in order to approximate f (t, η), η ∈ [η i , η i+1 [, in
To summarize, our conservative semi-Lagrangian algorithm, generalized to arbitrary displacement, as in [START_REF] Qiu | Conservative high order semi-lagrangian finite difference weno methods for advection in incompressible flow[END_REF], is based on the following steps:

1. The backward solution of the characteristics ending at the grid points

η i solution of dH dt = a J , H(t n+1 ) = η i , denoted by η * i = H(t n ; η i , t n+1
). These are the same characteristics as for the split 1D BSL method and the algorithm to compute them is the same. We shall consider here the backward Euler method.

2. The computation of the flux F i+1/2 . We denote by

η * i = η i * + α i ∆η, 0 ≤ α i < 1. Let P i be the polynomial of degree ≤ 3 satisfying P i (η i * +k ) = f n i * +k , k = -1, 0, 1, 2. We compute f i,k = 1 ∆η η * i+k η i * +k P i (η)dη, k = -1, 0, 1, 2,
by using Simpson rule for example. We then define

F i+1/2 = 7 12 (f i,0 +f i,1 )- 1 12 (f i,-1 +f i,2 )+ i+1≤k≤i * f n k - i * +1≤k≤i
f n k . (5.9)

3. The update of the solution. The new value of f at the grid points is then given by

f n+1 i = f n i + F i+1/2 -F i-1/2 . (5.10)
Note that in the derivation of the scheme, the assumption (5.7) corresponds to the special case i * = i. Formula (5.9) is chosen by considering that the local displacement is the decomposition of an exact displacement of a given number of cells plus a positive displacement less than one cell. We can check, as in [START_REF] Qiu | Conservative high order semi-lagrangian finite difference weno methods for advection in incompressible flow[END_REF] that the formula is continuous when α i approachs 1. We clearly see from the flux form (5.10), that the scheme is mass conservative. Note that in [START_REF] Qiu | Conservative high order semi-lagrangian finite difference weno methods for advection in incompressible flow[END_REF] explicit formulae are given for third and fifth order in the context of WENO. We consider here the example of fourth order and do not apply a WENO procedure. For constant displacement (with Jacobian equal to one), this scheme corresponds to use the BSL scheme with Lagrange interpolation of degree 4, where the stencil is chosen according to the sign of the displacement.

Free stream preservation

As for the BSL method, where the mass is only preserved up to first order in ∆t, we shall prove here that fi,j is preserved by the update formula up to first order in ∆t if f = 1 or equivalently fi,j = J i,j , which is our version of free stream preservation. For this, taking f n i = J i , we consider for simplifying the exposition of the proof that (5.7) occurs, but we could adapt the proof, removing the hypothesis. We use again the formulation of the update formula from the last section:

f n+1 i = f n i + (F * -F ) (η i ) = J i + (F * -F ) (η i )
and compute an approximation of (F * -F ) (η i ) in the limit of small ∆t. As we saw in the previous subsection F * -F is obtained by reconstruction of the grid function with sliding average

η * i ηi f h (η) dη = -∆t a(η i ) J(η i ) f h (η i ) + O(∆t 2 ) = -∆t a(η i ) + O(∆t 2 )
as we assume that f h (η i ) = J(η i ) = J i . Then denoting by S[a h ] the sliding average reconstruction of a h = (a(η i )) 1≤η≤N +1 , the update formula becomes for

f n i = J i f n+1 i = J i -∆t S[a h ](η i ) + O(∆t 2 ).
Let us now apply this formula to the two split steps of our method assuming a first order Lie splitting, but the same argument holds for higher order splittings. For the advection in the first direction we have for a given j, a = ∂ η2 Φ. Numerically Φ is approximated by its grid values Φ h = (Φ i,j ) 1≤i≤N1,1≤j≤N2 , we denote also by Φ h,j = (Φ i,j ) 1≤i≤N1 and Φ h,i = (Φ i,j ) 1≤j≤N2 one grid line at constant j and i respectively. The derivative in the η 2 direction is approximated by the sliding average reconstruction in the

η 2 direction S 2 [Φ h,i ]. Then a h,j = (S 2 [Φ h,i ](η j )) 1≤i≤N1 . So that our first update reads f * i,j = J i,j -∆t S 1 [S 2 [Φ h,i ](η j )](η i ) + O(∆t 2 ).
For the second split step, we first note that, as the reconstruction is linear the two terms above can be treated separately (and the higher order terms also). Then as there is a ∆t in factor of the second term, this will be a O(∆t 2 ) after reconstruction. So, as now a = -∂ η1 Φ, the formula is the symmetric of the other split step, which yields

f n+1 i,j = J i,j -∆t S 1 [S 2 [Φ h,i ](η j )](η i ) -S 2 [S 1 [Φ h,j ](η i )](η j ) + O(∆t 2 ) = J i,j + O(∆t 2 )
because of the linearity and the symmetry of the sliding reconstruction in the two directions. This proves that the constant states are preserved up to first order, which is needed for a consistent conservative approximation on curvilinear grids.

6 Numerical results

Rotation

The domain is Ω = [-π, π] 2 . We take

Φ(x 1 , x 2 ) = -((x 1 ) 2 + (x 2 ) 2 )/2
and solve (2.2). The mapping, referred to as deformed mesh (see e.g. [START_REF] Colella | High-order, finite-volume methods in mapped coordinates[END_REF]) is given by

x 1 (η 1 , η 2 ) = η 1 + α sin(η 1 2π/L x 1 ) sin(η 2 2π/L x 2 ), x 2 (η 1 , η 2 ) = η 2 + α sin(η 1 2π/L x 1 ) sin(η 2 2π/L x 2 ) (6.1)
for (η 1 , η 2 ) ∈ Ω = Ω and for 0 ≤ α < 1. Here, L x 1 = L x 2 = 2π. We use the cos-bell initial function defined by

f 0 (x 1 , x 2 ) = cos(r) 6 , for r < π/2, 0, else, with r = (x 1 -x 1 c ) 2 + (x 2 -x 2 c ) 2 and x 1 c = 1, x 2 c = -0.2.
We use a Verlet scheme for computing the characteristics backward in BSL2D, which leads to a second order in time scheme; for the evaluation of the fields, we use cubic splines in the fixed point algorithm (tolerance is put to 10 -12 an a maximum of 1000 iterations is allowed; when the convergence is not reached, a warning message is given). Periodic boundary conditions are used for this test, which can lead to small errors as the solution is not periodic but vanishes in the vicinity of the boundary. For the choice of the interpolation, we consider here only piecewise cubic polynomial reconstructions: SPL(3) and H(p). Note that LAG3 corresponds to H(3).

When not specified, we use the natural derivative for the potential, that is: the same derivative, as the derivative of the interpolation scheme. For H(p), when p is odd, the derivative is discontinuous; we choose the derivative of order p + 1 as natural derivative for the potential (we do not consider a discontinuous electric field here); this corresponds to taking the mean between the right and left derivative of order p.

We will refer to FD4 (resp. FD6) for computing the derivative using (4.4), with p = 4 (resp p = 6). We consider different discretizations in space: N = N 1 = N 2 ∈ {32, 64, 128, 256} and time: ∆t = 2 -i , i ∈ {0, 1, . . . , 10}. The schemes under consideration are here H2, H3, H4, H5, H6, H9, H10, H17 and SPL3. Numerical results are given on Figure 1 (solution and error at time T = 20), Figure 2 (error vs time step) and Figure 3 (evolution of mass). On Figure 1, we see the solution and error of H17, using ∆t = 2 -6 , N = 32 (top) and N = 64 (bottom) at time T = 20 that is after 10/π turns. We see that the error is reduced, taking a finer mesh and that the error is not located on the boundary, that is periodic boundary conditions do not lead to problems, in this test. The mesh is superposed to the solution, and we see that it is quite -For fixed value of N , by diminishing the time step, mass conservation is generally improved.

-For N = 64, we only see 3 types of curves: H2, H3 and the others which are almost the same. In order to see more differences on the different methods, we have to take a smaller value of N , like N = 32.

-We also can check that the centered reconstruction H2 (resp. H4) is worse for mass conservation than H3 (resp H5) when ∆t is big, but they become better, when ∆t is smaller, as expected from the analysis.

-We can also observe worse mass conservation using non natural reconstructions for the derivatives (not shown here), when N = 32, but the differences are slight, and we will see more difference in the guiding center case.

Kelvin-Helmholtz instability in a periodic box

We refer for example to [START_REF] Crouseilles | Conservative semi-lagrangian schemes for vlasov equations[END_REF] for this test case. This is a non linear test case, where the stream function Φ depends on f via the solution of the Poisson equation -∆Φ = f . The initial distribution f 0 is given by the formula :

f 0 (x 1 , x 2 ) = sin(x 2 ) + β cos(σx 1 )
where β = 0.015 and σ = 0.5. Periodic conditions are considered both in x 1 and x 2 direction. The domain is [0, L

x 1 ] × [0, L x 2 ], with L x 1 = 2π
σ , L x 2 = 2π. We use the same deformed mesh as in (6.1), with α = 0.9. We use mudpack [START_REF]Mudpack[END_REF] for the curvilinear Poisson solver. For BSL2D, we use a predictor-corrector scheme for the time loop as in [START_REF] Hamiaz | Guiding center simulations on curvilinear grids[END_REF], together with Verlet for the characteristics (as in the case of the rotation). For BSL1D, we replace the 2D advection, by 1D advection, using Strang splitting, as in [START_REF] Crouseilles | Conservative semi-lagrangian schemes for vlasov equations[END_REF]; a trapezoidal rule is used for the 1D characteristics, and again cubic splines for the evaluation of the fields, in the fixed point algorithm. For CSL, we compute the characteristics as for BSL1D and use the same Strang splitting. The schemes under consideration are here SPL3, SPL3-FD6 (SPL3 with FD6 instead of cubic splines for the derivatives of the potential), SPL3-1d (BSL1d with SPL3), H6, H6-SPL3 (H6 with cubic splines instead of FD6 for the derivatives of the potential), H17, H17-1d (BSL1d with H17) and CSL.

Numerical results for BSL are given on Figure 4 (solution at time T = 50) and Figure 5 (evolution of mass). Numerical results for CSL are given on Figure 6 (solution at time T = 50) and Figure 7 (evolution of free streaming error). On Figure 7, L 2 norm et energy evolution are also shown for CSL and BSL.

On Figure 4, we represent f at time T = 50 using H17 (left) and SPL3 (right), with N = 256, ∆t = 2 -10 (top), N = 128, ∆t = 2 -10 (middle) and N = 128, ∆t = 2 -5 (bottom). Using N = 128, we have less details than with N = 256 for H17. This scheme leads to similar results, using different small enough time steps for a given N . Note that if we would take bigger time steps very bad mass conservation (here, there is no picture; the result is similar to α = 0.5, but with values around 1 for ∆t = 2 -4 and between 6 and 7 for ∆t = 2 -4 at final time T = 100). We see on Figure 10 that for α = 0.9, the quantities are changed, especially for ∆t = 2 -3 , and also for the CSL method. Note that for ∆t = 2 -3 , the simulation stops in the CSL method, as we have added a test that prevents from having intersection of characteristics. We see that the energy is less accurately conserved for ∆t = 2 -4 , for the CSL method, whereas it was not changing much from α = 0 to α = 0.5. For the BSL method, the energy conservation is really worse for ∆t = 2 -3 . On the other hand, we remark that H17 performs the best for ∆t = 2 -4 in comparison to the other methods. A picture of the meshes for α = 0.25 and α = 0.5 is depicted on Figure 11, using N = 32 instead of N = 256 in order to better see the mesh.

Conclusion

We have studied the semi-Lagrangian method on curvilinear grids and looked at mass and constant states preservation issues. A (mass) conservative method that preserves constant states up to first order in time has been exhibited and compared with the classical advective backward semi-lagrangian method, which in turn always preserves constant states and for which we have shown that it conserves mass up to first order in time. When the mesh becomes distorted, the time step has to be small enough; otherwise important errors appear, as in particular characteristics are badly solved. In this setting, Hermite interpolation with high order odd reconstruction of the derivatives exhibits more favorable behavior than cubic splines, but the Hermite method can be also more expensive; so we may stick to cubic splines when the mesh is not too distorted. The next step is to adapt such methods to gyrokinetic simulations. The classical BSL method may be preferred, as we can make it work for several reconstructions and it is more robust, even though there is a loss of mass conservation. The CSL method could also be tried, but we fear for stability and robustness issues, as the design of such a scheme was already not easy in this simplified context. We have tested the methods for a quite severely deformed mesh; we have seen that, when the mesh is less deformed, the time constraint is relaxed. The grid in a gyrokinetic simulation is imposed by the background magnetic field, which does not vary in time. Hence the methodology described in this paper can be directly applied as we can define a transformation, based on a B-spline representation, which maps a logical domain onto the desired physical domain. Typical gyrokinetic simulations have magnetic configurations not far from equilibrium, thus including the fluctuations of the magnetic field in the grid information is not necessary for most applications. Among the new difficulties that we can encounter in 5D gyrokinetic simulations, one is due to the dimension, which leads to heavy computations; the model can also lead to numerical instabilities, but collisions can be added to prevent from this fact; the whole system can also become quite complex with all the physical terms; boundary conditions and change of geometry have also to be treated correctly. A possible further study is also to add slope limiters or other tools in order to prevent from sporadic out of bounds values which can destroy the solution when they become too high.

  . , N + d, the B-spline B d+1 j , of degree d and order d + 1, whose support is [τ j , τ j+d+1 ]. For all m ∈ N * such that τ m < τ m+1 and m ≥ d + 1 and m ≤ N + d, the B-splines that are non-zero in the interval [τ m , τ m+1 [ are B d+1 , = m -d, . . . , m, with the exception that B d+1 m (τ m ) = 0. For x ∈ [τ m , τ m+1 [, we can compute b j = b j (x), j = 1, . . . , d + 1, where b j (x) = B d+1 m-d-1+j (x), by the following algorithm:

  12, -2/3; 0, 2/3, -1/12), p = 5 b + = (1/20, -1/2; -1/3, 1, -1/4, 1/30), p = 6 b + = (-1/60, 3/20, -3/4; 0, 3/4, -3/20, 1/60).(4.4)

Figure 1 :

 1 Figure 1: Rotation case. Solution at time T = 20, on deformed mesh: H17 α = 0.9, ∆t = 1/2 6 (left: solution, right: error), with N = 32 (top) and N = 64 (bottom).

Figure 2 :

 2 Figure 2: Rotation case. Average of the error in maximum norm between t = 10 and t = 20 divided by t versus i, where ∆t = 1/2 i .

Figure 3 :

 3 Figure 3: Rotation case. Evolution of mass.

Figure 4 :

 4 Figure 4: Guiding center case. Solution at time T = 50, on deformed mesh, with α = 0.9 (left: H17, right: SPL3), N = 256, ∆t = 2 -10 (top), N = 128, ∆t = 2 -10 (middle) and N = 128, ∆t = 2 -5 (bottom). Mesh for N = 256 (top right) and N = 128 (middle right) are also shown.

Figure 5 :

 5 Figure 5: Guiding center case. Evolution of mass.

Figure 6 :

 6 Figure 6: Guiding center case. Solution at time T = 50, on deformed mesh for CSL, with α = 0.9 (left: N = 128, right: N = 256), ∆t = 2 -5 (top), ∆t = 2 -6 (middle) and ∆t = 2 -7 (bottom).

Figure 8 :

 8 Figure Guiding center case. Evolution of free stream error: FD4 for field (top), FD6 (middle left) SPL3 (middle right). Evolution of L 2 norm (bottom left) and energy (bottom right).

Figure 9 :

 9 Figure 9: Guiding center case. Evolution of mass for α = 0, 0.1, 0.25, 0.5; N = 256 (from left to right, top to bottom).

Figure 10 :

 10 Figure 10: Guiding center case. Evolution of energy (left) and L 2 norm (right) for α = 0, 0.5, 0.9 (from top to bottom).

Figure 11 :

 11 Figure 11: Guiding center case. Representation of mesh, on 32 × 32 grid, for α = 0.25 and α = 0.5.

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme

2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

deformed. Note that the Jacobian becomes singular for α = 1. The higher α is, the more difficult it is to get accurate solution (see [START_REF] Hamiaz | Guiding center simulations on curvilinear grids[END_REF]). Here we always take α = 0.9 to see the robustness of the different methods. On Figure 2, we give the error in function of i, with ∆t = 2 -i using the different reconstructions and N = 32 (top left), N = 64 (top right), N = 128 (bottom left) and N = 256 (bottom right). The error is computed in L ∞ norm taking the best constant that approaches in least square sense the error between T = 10 and T = 20 (we use the fit command of gnuplot). The error is multiplied by (N/32) 3 as we expect third order convergence in space for the different schemes (except H2). So, we observe how the error behaves with respect to time for different discretizations in space. One can expect to have an error of the form

(see [START_REF] Charles | Enhanced convergence estimates for semi-lagrangian schemes application to the vlasov-poisson equation[END_REF] for example). We distinguish basically three zones. A first zone (i = 0, 1, 2, 3) when ∆t is too big to observe convergence. Typically the time step is too big for solving properly the characteristics which are quite stiff to solve, as the mesh is very distorted. Then we observe a second zone, where the error decreases like ∆t 2 . This zone becomes bigger as N increases. Finally, there is a third zone, where the error saturates. Multiplying by (N/32) 3 permits to have the error for N = 32 and to see how the "constant" in the space error behaves (which is C 2 if the error is of the form (6.2)). For N = 32 (top left), we observe that increasing the value of p in H(p) increases the accuracy and SPL3 is little better that H6. Moreover, when ∆t becomes smaller, we see bad effect for the schemes that use a centered reconstruction for the derivatives: SPL3 and H10 (looking at the solution, we see that the error comes from the boundary; using other boundary conditions may lead to improvements); for H6, which is a little more diffusive, it is not observed here. When i <= 3 the error is big and does not depend on the reconstruction. H2 gets worser increasing N , as it is not a third order scheme; we also observe that taking smaller time steps leads to bigger error, when time error is not predominant. For H3, the error saturates or little increase with ∆t. We remark that is of order 3, comparing N = 128 and N = 256. For the other schemes, increasing N , the difference between the schemes is diminished, as the reconstruction for the derivatives is higher than third order (for cubic splines it is fifth order, and the constant seems to be smaller than for H5).

On Figure 3, we represent the evolution of mass, for the different reconstructions, with N = 32 (left), N = 64 (right), and ∆t = 2 -5 (top), ∆t = 2 -6 (middle) and ∆t = 2 -7 (bottom).

-We observe an increase of error of mass, when the bell crosses the region where the mesh is the most distorted.

-Increasing N improves the mass for the lower accuracy methods (essentially H2, H3, H4 and H5); it improves also the mass for all the methods, when the bell is not touching the most distorted region.

(like ∆t = 2 -3 which can be taken in the case of uniform mesh with α = 0), the error in time would increase and we would get different (wrong) results. We see that this scheme is particularly robust with respect to ∆t. For example, for the fixed value of N = 128 (the exercise could be repeated for N = 256), using ∆t = 2 -5 or ∆t = 2 -10 leads to similar result. On the other hand, SPL3, works very badly when the time step is small; but, using a bigger time step as ∆t = 2 -5 , for N = 128, gives an acceptable result almost as accurate as H17 for the same value of N . Indeed, we see some more details in H17, which are also present on the "reference solution" (H17, N = 256 top left). On Figure 5, we represent the evolution of mass for N = 256, ∆t = 2 -5 (top left), ∆t = 2 -6 (top right), ∆t = 2 -7 (middle left), ∆t = 2 -8 (middle right), ∆t = 2 -9 (bottom left) and ∆t = 2 -10 (bottom right).

-We remark that BSL1d (here SPL3-1d and H17-1d) leads to huge mass error, when the time step is big (see Figure top left); this confirms previous analysis on uniform grid (see [START_REF] Crouseilles | Conservative semi-lagrangian schemes for vlasov equations[END_REF] for example); on the other hand, diminishing the time step (i = 5 to i = 9) BSL1d and BSL2d become similar: SPL3-1d converges to SPL3, and similarly H17-1d to H17.

-Diminishing the time step (i = 5 to i = 9), we get better conservation of mass, except for the schemes that do not use the natural reconstruction of the derivatives (SPL3-FD6 and H6-SPL3). This confirms the theoretical analysis of Section 4.3, which was more difficult to point out in the rotation case.

-For i = 10, we have first better conservation of mass for SPL3 than for H17 (which was the best for i ≤ 9) this is coherent with the theoretical analysis of Section 4.3. On the other hand, the solution develops a lot of oscillations (as shown on Figure 4 top right) and a numerical instability appears. Note that instability appears before for SPL3-1d than for SPL3. For H6, we also get an instability (which was not in the rotation case).

On Figure 6, we represent f at time T = 50 using CSL, with N = 128 (left) and N = 256 (right), for ∆t = 2 -5 (top), ∆t = 2 -6 (middle) and ∆t = 2 -7 (bottom). The solution seems to converge to the "reference solution" (H17, N = 256 top left). Sporadic out of bounds value of f appear and become bigger by taking a smaller time step and these values are attenuated by refining the space grid. We remark some slight grid effect, taking the biggest ∆t (that is ∆t = 2 -5 ); this may be due to the fact that the characteristics are not computed accurately enough. Some numerical oscillations are present and in the mean time some details are diffused; this may be explained from the order of the scheme and the behaviour as ∆t tends to 0 (here fourth order).

On Figure 7, we represent the free stream error. To compute it, we consider a function initialized to 1 and let it evolve as f with the same advection field and the same scheme. As the field is divergence free, the function should remain 1, as it is the case for BSL. The free stream error is then the L ∞ error between the solution and 1. It is represented versus time. On top left, we multiply the error by (2 -5 /∆t) 2 and otherwise not. It is represented on top for CSL with the natural reconstruction of the field, that is, FD4. On Figure 7 middle, we change the reconstruction to FD6 (left) and SPL3 (right). On Figure 7, bottom, we represent the evolution of L 2 -norm (left) and energy (right) for CSL, and BSL (H3,H6,H17 and SPL3), using ∆t = 2 -5 , or ∆t = 2 -6 and N = 256. As seen on Figure 7 top left, we check the order 2 accuracy of the free stream preservation. We can notice, that is holds both in the linear phase, when everything is smooth and remains in the non linear phase, where we are far from convergence. On the other hand, the other reconstructions of the fields do not lead to this property. Taking ∆t = 2 -7 , N = 128, we even observe that the scheme can become unstable (Figure 7 middle left). So the natural reconstruction of the field seems here to prevent from the appearance of instability; the sporadic out of bounds values can become higher and higher; so that at a certain level, this destroys the solution. Other ingredients may be needed for garantying the stability of the scheme and to permit to go to higher order schemes. On Figure 7, bottom, it is confirmed that diffusion is not as strong as H3, but more than other usual BSL schemes (H6, SP L3, H17). We also see that energy is quite well preserved for CSL, certainly because the scheme is mass conservative; on the other hand, diminishing the time step also improves energy conservation for the BSL method (except for H3); this is coherent with the analysis about mass conservation.

On Figure 8, we represent the free streaming error, the mass, the L 2 norm and the energy, in the cartesian case (i.e. α = 0). We observe that free streaming error is less important for the conservative method (Figure 8, top left). We remark that it is more important for N = 256 than N = 128, when ∆t = 2 -5 . This may be due to the fact that the solution is more complex when N = 256. On the other hand, diminishing the time step, we observe that the free streaming error diminishes, as expected, because we use FD4 for the reconstruction of the derivatives. Mass conservation is really improved, using α = 0 (Figure 8 top right) for the BSL methods. The effect of diminishing the time step leads to small improvement of the mass conservation, in constrast to the curvilinear case (α = 0.9), but the conservation is already a lot better. We remark no real difference for the L 2 norm (Figure 8, bottom left) and the energy is strangely not conserved as well; the curve is however less oscillating. On Figure 9, we represent the mass evolution, for bigger time steps: ∆t = 2 -3 or ∆t = 2 -4 , which corresponds to more standard time steps that are used for this test case in cartesian geometry (see [START_REF] Crouseilles | Conservative semi-lagrangian schemes for vlasov equations[END_REF] e.g.) and for α ∈ {0, 0.1, 0.25, 0.5}. We remark that the mass is really better conserved on the cartesian mesh (α = 0), and as α increases, mass conservation is degraded (for non conservative methods, of course). This can be seen even for α = 0.1 which is almost uniform. As already shown, diminishing the time step has a beneficial effect, which is not so clear for α = 0 or α = 0.1 but which becomes clearer increasing the value of α. L 2 norm and energy evolutions are plotted on Figure 10 for α ∈ {0, 0.5, 0.9}. On these quantities which are less accurately conserved in the cartesian case, we do not see differences between the cartesian case α = 0 and the case α = 0.5. This is good news, as it tells us that the time step restriction is not so severe; diminishing the time step mainly helps in the conservation of the mass. In the case of α = 0.9, for which reasonable simulation are shown for ∆t ≤ 2 -5 , exhibit


Rafik Dey

Laurent Cavalié

Christine Vernet

Jacques Bodennec

Pierre Pernin
email: pernin@sante.univ-lyon1.fr

The free-living amoeba Willaertia magna, is particularly resistant to infection by the pathogenic bacteria Legionella pneumophila

Keywords: Willaertia magna, Legionella pneumophila, free-living amoebas, protozoabacterial interactions

Legionella pneumophila, the causative agent of Legionnaire's disease, is well characterized as a bacteria surviving and developing, almost exclusively, as intracellular parasite within freshwater protozoa. Several species of protozoa and ciliae have been shown to support the growth of the pathogenic bacteria. In the present study, we report for the first time the behaviour of the protozoan Willaertia magna towards L. pneumophila and compared it with Acanthamoeba castellanii and Hartmannella vermiformis, two known L. pneumophila permissive protozoa. The results show that Willaertia amoebas displayed the ability to internalize the pathogenic bacteria but at a different extent when compared to other species. Surprisingly, and at the opposite of the other two amoebic species, coculture experiments showed that Willaertia magna inhibited the growth and the development of L. pneumophila and that it was resistant to bacterial induced cytotoxicity. Electron microscopy demonstrated that the formation of replicative phagosomes observed within Acanthamoeba sp. and Hartmannella sp. is deficient in Willaertia sp. These observations describing for the first time the occurrence of a protozoan species that display resistance towards the bacteria L. pneumophila, shed new light on a cellular model that may provide useful information on the mechanisms underlying bacterial infection and show that among free-living amoebas, the different species are not equally permissive to L. pneumophila.

INTRODUCTION

The gram negative bacteria Legionella pneumophila, the causative agent of the Legionnaire's disease, is characterized by a facultative intracellular replication inside macrophages, monocytes or epithelial cells [START_REF] Horwitz | The Legionnaires' disease bacterium (Legionella pneumophila) inhibits phagosome-lysosome fusion in human monocytes[END_REF]). In the environment, L. pneumophila presents a ubiquitous aquatic repartition and replicates inside protozoa as free-living amoebas [START_REF] Molmeret | Cell biology of the intracellular infection by Legionella pneumophila[END_REF][START_REF] Molmeret | Amoebae as training grounds for intracellular bacterial pathogens[END_REF]. In all these phagocytic cells, the genes of the dot/icm system of L. pneumophila are essential to warrant a successful infection with evasion of the endocytic maturation and inhibition of the phagosome lysosome-fusion inside the diverse cellular hosts [START_REF] Roy | The dot/Icm transporter of Legionella pneumophila: a bacterial conductor of vesicle trafficking that orchestrates the establishment of a replicative organelle in eukaryotic hosts[END_REF]. Rowbotham was the first to demonstrate an intracellular replication of the bacteria in free-living amoebae [START_REF] Rowbotham | Preliminary report on the pathogenicity of Legionella pneumophila for freshwater and soil amoebae[END_REF]. Since this previous work, numerous studies have been performed that strengthen the existence of a close relationship between Legionella and free-living amoebas. Hence, co-culture experiments in liquid media have shown an important growth promoting effect of amoebas towards L. pneumophila. During outbreaks of Legionnaire's disease, water sources are usually found to contain significant amount of free-living amoebas, and conversely, high numbers of Legionella species are found when amoebic concentrations are high [START_REF] Barbaree | Isolation of protozoa from water associated with a legionellosis outbreak and demonstration of intracellular multiplication of Legionella pneumophila[END_REF][START_REF] Fields | Intracellular multiplication of Legionella pneumophila in amoebae isolated from hospital hot water tanks[END_REF][START_REF] Breiman | Association of shower use with Legionnaires' disease. Possible role of amoebae[END_REF]. In addition, the isolation of L. pneumophila from water is strongly improved by addition of amoebas [START_REF] Sanden | Incubation of water samples containing amoebae improves detection of legionellae by the culture method[END_REF]. Moreover, Steinert [START_REF] Steinert | Resuscitation of viable but nonculturable Legionella pneumophila Philadelphia JR32 by Acanthamoeba castellanii[END_REF] showed that the addition of amoebas can trigger the revival of viable Legionella strains that became non cultivable (VBNC state) due to a prolonged stage in poor medium such as distilled water.

Although Fields reported 13 species of amoebas and 2 species of ciliated protozoa that display the ability to sustain L. pneumophila replication [START_REF] Fields | The molecular ecology of legionellae[END_REF], most of in vitro experiments have been performed using two amoebic genera: Acanthamoeba [START_REF] Anand | Interaction of L. pneumophila and a free living amoeba (Acanthamoeba palestinensis)[END_REF][START_REF] Holden | Intracellular growth of Legionella pneumophila within Acanthamoeba castellanii Neff[END_REF][START_REF] Vandenesch | Relationship between free amoeba and Legionella: studies in vitro and in vivo[END_REF][START_REF] Moffat | A quantitative model of intracellular growth of Legionella pneumophila in Acanthamoeba castellanii[END_REF][START_REF] Bozue | Interaction of Legionella pneumophila with Acanthamoeba castellanii: uptake by coiling phagocytosis and inhibition of phagosome-lysosome fusion[END_REF][START_REF] Gao | Utilization of similar mechanisms by Legionella pneumophila to parasitize two evolutionarily distant host cells, mammalian macrophages and protozoa[END_REF][START_REF] Neumeister | Multiplication of different Legionella species in Mono Mac 6 cells and in Acanthamoeba castellanii[END_REF][START_REF] Harb | Heterogeneity in the attachment and uptake mechanisms of the Legionnaires' disease bacterium, Legionella pneumophila, by protozoan hosts[END_REF] and Hartmannella [START_REF] King | Effects of cytochalasin D and methylamine on intracellular growth of Legionella pneumophila in amoebae and human monocyte-like cells[END_REF][START_REF] Kwaik | Protein expression by the protozoan Hartmannella vermiformis upon contact with its bacterial parasite Legionella pneumophila[END_REF]Abu Kwaik, 1996.;Abu Kwaik et al., 1998a;[START_REF] Harb | Heterogeneity in the attachment and uptake mechanisms of the Legionnaires' disease bacterium, Legionella pneumophila, by protozoan hosts[END_REF] and more recently Dictyostelium, which is not really a freeliving amoeba [START_REF] Hagele | Dictyostelium discoideum: a new host model system for intracellular pathogens of the genus Legionella[END_REF][START_REF] Solomon | Intracellular growth of Legionella pneumophila in Dictyostelium discoideum, a system for genetic analysis of host-pathogen interactions[END_REF]. But, few studies have been carried out with other amoebic genera. For these reasons, cocultures of L. pneumophila with Willaertia magna, an amoebic species never tested until now, were performed and compared with two bacteria permissive amoebas : A. castellanii and H. vermiformis.

MATERIAL AND METHODS

Strains

L. pneumophila (Paris serogroup 1 CIP 107 629T), was cultured at 37°C for 4 to 5 days on buffered charcoal yeast extract (BCYE) agar before coculture experiments in order to be in post-exponential phase.

Willaertia magna (deposit of this strain was performed at ATCC under the number PTA-7824), Hartmannella vermiformis (Ax.5.2e4b) and Acanthamoeba castellanii (By 02.2.4) amoebas are environmental strains isolated from a French thermal spa by filtration of water, and grown at 30°C on non nutrient agar (NNA) overlaid with a thin film of Escherichia coli. These three strains were established in axenic culture at 37°C in SCGYEM liquid medium containing 10% of foetal calf serum.

Coculture of L. pneumophila with the different amoebas

Tubes (FALCON ® 3033) containing 3 ml of SCGYEM medium were seeded with 5.5×10 4 trophozoites/ml of the different amoebic strains maintained in exponential growth phase by subcultures every 3-4 days. L. pneumophila, grown on BCYE medium, was suspended in sterile distilled water at 10 9 /ml (= 1 OD at 550 nm) and, after appropriate dilutions, was inoculated at a multiplicity of infection (MOI) of 50, i.e around 2.8 x 10 6 bacteria/ml. The tubes were immediately centrifuged at low speed (5 min at 760 g) to initiate interactions between bacteria and amoebas. After 10 min, the pellet was resuspended and the tubes incubated in a slanting position at 37°C. The cocultures were analyzed for up to 5 days (day 0 to day 4) as follows. At the indicated time, the tubes were placed on ice for 6 min to reduce cell adherence and, after a vigorous agitation by vortex, amoebas were numbered using a hematocytometer. Controls of amoebic growth were also performed in the same conditions in absence of bacteria. The number of total L. pneumophila (intra-and extracellular), expressed in CFU/ml, was determined after serial 10 fold dilutions of the coculture medium with sterile H 2 O that were spread in triplicate on BCYE plates and were further incubated at 37°C for at least 6 days.

Statistical differences were addressed using Student t-test.

Cytotoxicity of L. pneumophila

Firstly, the cytotoxicity of L. pneumophila towards the different strains of amoebas was studied qualitatively by phase contrast microscopy in 24 well plates containing 5×10 4 amoebas/well infected at an MOI of 50 and amoebic monolayer formation was observed after 72 hr infection with L. pneumophila.

Bacterial cytotoxicity was also determined on Acanthamoeba and Willaertia species by trypan blue exclusion test as follows: at 48 and 72 hr of infection (MOI 50), after agitation of the coculture medium, the amoebic suspension was pelleted by low centrifugation and suspended in 200 μl of trypan blue-SCGYEM medium (4:1 v/v) mix.

The extent of cell death was determined as the percentage of trypan blue positive cells. This assay could not be performed with H. vermiformis due to the punctiform, shrivelled appearance of the remaining cells resulting from the high cytotoxicity of L. pneumophila 72 hours after infection. Trypan blue exclusion tests were also performed on control amoebic cultures grown in the same conditions but without L. pneumophila. Statistical differences were addressed using Student t-test.

Transmission electron microscopy

Axenic cultures of H. vermiformis, A. castellanii and W. magna were infected for 36 hr with L. pneumophila at a MOI of 50 and 100. After decantation of the medium, the amoebas were fixed for 20 min at room temperature using 2% glutaraldehyde in serum free SCGYEM medium. After removal of the medium, the fixation of cells was further completed for 30 min with 2% glutaraldehyde in 0.1 M sodium cacodylate buffer (pH 7.4).

After scrapping, amoebas were transferred in 1.5 ml Eppendorf microtubes, washed three times (15 min each) by repeated cycles of suspension-centrifugation, and post-fixed for 40 min with 1% OsO 4 in 0.15 M sodium cacodylate buffer. After dehydration of the samples in ethanol solutions of increasing concentration (50°, 70°, 90°, 95°, 2 x 100°), the cells were embedded in EPON resin. Ultrathin sections were stained with uranyl acetate followed by lead citrate and examined using a JEOL 1200 CX electron microscope at 80 kV. A minimum of more than 100 cells was observed and the percentage of amoebas displaying a replicative phagosome (i.e vacuoles filled with several bacterias) determined.

RESULTS

We first studied the ability of L. pneumophila to infect and growth within the different strains of amoebas. As expected, and in accordance to what was previously published [START_REF] Anand | Interaction of L. pneumophila and a free living amoeba (Acanthamoeba palestinensis)[END_REF][START_REF] Holden | Intracellular growth of Legionella pneumophila within Acanthamoeba castellanii Neff[END_REF][START_REF] Vandenesch | Relationship between free amoeba and Legionella: studies in vitro and in vivo[END_REF][START_REF] King | Effects of cytochalasin D and methylamine on intracellular growth of Legionella pneumophila in amoebae and human monocyte-like cells[END_REF][START_REF] Moffat | A quantitative model of intracellular growth of Legionella pneumophila in Acanthamoeba castellanii[END_REF][START_REF] Kwaik | Protein expression by the protozoan Hartmannella vermiformis upon contact with its bacterial parasite Legionella pneumophila[END_REF][START_REF] Bozue | Interaction of Legionella pneumophila with Acanthamoeba castellanii: uptake by coiling phagocytosis and inhibition of phagosome-lysosome fusion[END_REF]Abu Kwaik et al., 1996.;[START_REF] Gao | Utilization of similar mechanisms by Legionella pneumophila to parasitize two evolutionarily distant host cells, mammalian macrophages and protozoa[END_REF][START_REF] Neumeister | Multiplication of different Legionella species in Mono Mac 6 cells and in Acanthamoeba castellanii[END_REF]Abu Kwaik et al., 1998a), the bacteria displayed the ability to grow within H. vermiformis and A. castellanii with a respective increase of 1.3 and 1.4 log in CFU/ml after 3 days of infection (Fig. 1).

Surprisingly, and at the opposite of these two amoebas, W. magna did not favour bacterial yield and even reduced it significantly within the coculture medium (Fig. 1). A 2.3 log difference in bacterial yield was measured at day 2 post infection between the Hartmannella-Acanthamoeba species and W. magna.

The observation that amoebic species display a different behaviour towards L. pneumophila was further strengthened by the fact that Hartmannella and Acanthamoeba growths were particularly affected by the bacteria when compared to Willaertia amoebas (Table 1). Firstly, a significant decrease in Hartmannella and Acanthamoeba number was detected on day 3 of coculture with L. pneumophila when compared to day 1, a phenomenon that could not be observed with Willaertia. Secondly, a significant drop was observed in Hartmannella and Acanthamoeba growth cocultured with Legionella when compared to their counterparts cultured in absence of the bacteria (control amoebas).

Hence, after 3 days of coculture, the number of Hartmannella and Acanthamoeba was reduced by 82 and 93 % (P< 0.001) when compared to the respective amoebas cultured without L. pneumophila; a decrease was also observed in the number of W. magna cocultured with the bacteria when compared to bacterial-free control, but it was limited to 25 % and less significant (P<0.05).

We then examined the putative cytotoxicity of L. pneumophila towards the different species of amoebas. Phase contrast microscopic observations of microplates clearly shown alterations occurring within Hartmannella and Acanthamoeba populations while Willaertia magna was not affected (Fig. 2). Hence, after 72 hr of bacterial infection, few adherent Hartmannella and Acanthamoeba amoebas effectively remained when compared to controls, while Willaertia still formed a dense monolayer of cells similar to control.

Following infection, Hartmannella and Acanthamoeba trophozoites became progressively rounded, shrivelled, and lost their adherence ability after 2 days in presence of the bacteria.

Moreover, further experiments shown that 28.4 ± 8% of the very few Acanthamoeba amoebas remaining at 72 hr versus only 3.8 ± 0.9% of the numerous W. magna were trypan blue positive (Fig. 3), showing again a particular resistance of this latter species to bacterial induced cytotoxicity when compared to A. castellanii.

We next determined whether the defective growth of L. pneumophila in presence of W. magna could be due to a defect in intracellular replication. The number of amoebas displaying a replicative phagosome was determined by electron microscopy. In Acanthamoeba and Hartmannella amoebas, the percentage of cells displaying a replicative phagosome at 36 hr was very high, i.e., 36 % to 51% according to the MOI (Table 2). In contrast, replicative phagosomes were not expressed in Willaertia amoebas infected at a MOI of 50, although some rare bacteria were sometimes internalized as demonstrated in figures 5A and 5B. When L. pneumophila was cocultured with W. magna at a higher ratio (MOI of 100), few amoebas displayed a replicative phagosome but at a very low level (1.5%) when compared to Acanthamoeba and Hartmannella (Table 2). Hence, electron microscopy showed that in A. castellanii and H. vermiformis, all the stages of infection could be observed. Shortly after infection, L. pneumophila appeared in these two species as isolated bacteria in a vacuole-like structure (Fig. 4A) surrounded by rough endoplasmic reticulum (RER) according to previous descriptions (Abu Kwaik, 1996.). More pronounced invasion leads to the formation of clear distended replicative phagosomes containing numerous bacteria juxtaposed with cytoplasmic amorphous elements (Fig. 4B) and ended by the necrosis of these heavily infected amoebas. Nevertheless, even in these heavily infected amoebas, the contractile vacuole was still active (Fig. 4C).

In contrast, out of the few W. magna displaying intracellular bacteria (only 2/133

Willaertia amoebas examined at a MOI of 100), the replicative phagosome showed morphological differences when compared to those observed in Acanthamoeba or Hartmannella amoebas. These phagosomes were not clearly delimited by a membrane as in A. castellanii (compare Fig 5C and 5D with Fig. 4B) and the bacteria inside contain numerous vacuoles usually described as polyhydroxybutyrate granules (Fig. 5D). In Acanthamoeba and Hartmannella amoebas, the bacteria displayed the characteristic morphology of the typical gram-negative cell envelope composed of a wavy outer membrane around a clear periplasmic space (Fig. 6A), while the L. pneumophila observed within the two infested Willaertia amoebas appeared as short stubby rods with a dense cytoplasm delimited by a thick laminar outer layer (Fig. 6B). These observations suggest that the replication of L. pneumophila cannot occur properly within W. magna and may explain why bacterial growth is inhibited in presence of this amoeba (Fig. 1).

DISCUSSION

The main finding of these experiments is that W. magna displays a particular resistance to infection by L. pneumophila when compared to the amoebas, A. castellanii and H. vermiformis. To our knowledge, this is the first study that systematically addresses the ability of Willaertia amoebas to resist against L. pneumophila. It also confirms previous studies showing that the Acanthamoeba and Hartmannella species support the proliferation of L. pneumophila. Numerous studies effectively shown that multiplication of L. pneumophila occurs within various protozoa such as amoebas or ciliates which make that free-living amoebas are now considered as a natural reservoir for the pathogenic bacteria (Abu Kwaik et al., 1998b;[START_REF] Molmeret | Cell biology of the intracellular infection by Legionella pneumophila[END_REF][START_REF] Molmeret | Amoebae as training grounds for intracellular bacterial pathogens[END_REF]. Fields reported that five amoebic genera are able to support the intracellular growth of L. pneumophila [START_REF] Fields | The molecular ecology of legionellae[END_REF]. However, out of these 5 genera, only two (Acanthamoeba and Hartmannella) have been repeatedly tested in coculture with L. pneumophila [START_REF] Anand | Interaction of L. pneumophila and a free living amoeba (Acanthamoeba palestinensis)[END_REF][START_REF] Holden | Intracellular growth of Legionella pneumophila within Acanthamoeba castellanii Neff[END_REF][START_REF] Vandenesch | Relationship between free amoeba and Legionella: studies in vitro and in vivo[END_REF][START_REF] King | Effects of cytochalasin D and methylamine on intracellular growth of Legionella pneumophila in amoebae and human monocyte-like cells[END_REF][START_REF] Moffat | A quantitative model of intracellular growth of Legionella pneumophila in Acanthamoeba castellanii[END_REF][START_REF] Kwaik | Protein expression by the protozoan Hartmannella vermiformis upon contact with its bacterial parasite Legionella pneumophila[END_REF]Abu Kwaik, 1996;[START_REF] Bozue | Interaction of Legionella pneumophila with Acanthamoeba castellanii: uptake by coiling phagocytosis and inhibition of phagosome-lysosome fusion[END_REF][START_REF] Gao | Utilization of similar mechanisms by Legionella pneumophila to parasitize two evolutionarily distant host cells, mammalian macrophages and protozoa[END_REF][START_REF] Neumeister | Multiplication of different Legionella species in Mono Mac 6 cells and in Acanthamoeba castellanii[END_REF]Abu Kwaik et al., 1998a;[START_REF] Harb | Heterogeneity in the attachment and uptake mechanisms of the Legionnaires' disease bacterium, Legionella pneumophila, by protozoan hosts[END_REF]. Recently, it was shown that Dictyostelium discoideum is also able to support the bacterial growth and proliferation [START_REF] Hagele | Dictyostelium discoideum: a new host model system for intracellular pathogens of the genus Legionella[END_REF][START_REF] Solomon | Intracellular growth of Legionella pneumophila in Dictyostelium discoideum, a system for genetic analysis of host-pathogen interactions[END_REF]. Our results with A. castellanii and H. vermiformis are in agreement with these previous studies although the bacteria grow to a lesser extent in our experiments when compared to available data in the literature. This may be explained by the fact that bacterial growth determination was performed upon spontaneous lysis of amoebas instead of the artificial lysis that was often used in other studies. However, Acanthamoeba and Hartmannella species displayed a 2.3 log higher ability to support bacterial growth when compared to W. magna tested exactly in the same conditions. Concomitantly to the lack of L. pneumophila replication, the bacteria did not display any obvious cytotoxicity towards W. magna, demonstrating again the particular resistance of this latter species when compared to the other two permissive genera.

Different possibilities, related to either the bacteria or the amoeba itself, may explain the particular resistance of W. magna. Until now, in all studies that reported a defect of L. pneumophila growth in coculture with macrophages and/or amoebas, the underlying mechanism results from the use of mutant strains of the bacteria that split in two great categories. Several bacterial mutants of the dot/icm system are defective in evasion of phagolysosomal fusion and consequently, have lost their capacity of intracellular replication [START_REF] Cianciotto | Legionella pneumophila mip gene potentiates intracellular infection of protozoa and human macrophages[END_REF][START_REF] Swanson | Identification of Legionella pneumophila mutants that have aberrant intracellular fates[END_REF][START_REF] Gao | Utilization of similar mechanisms by Legionella pneumophila to parasitize two evolutionarily distant host cells, mammalian macrophages and protozoa[END_REF][START_REF] Coers | Identification of icm protein complexes that play distinct roles in the biogenesis of an organelle permissive for Legionella pneumophila intracellular growth[END_REF]. In contrast, the rib mutants are defective in their cytotoxic effect and are unable to egress from the phagocytic cell despite an active intracellular replication [START_REF] Alli | Temporal pore formation-mediated egress from macrophages and alveolar epithelial cells by Legionella pneumophila[END_REF]. However, these two possibilities can be excluded in our study since the bacterial strain of L. pneumophila used in coculture experiments with W. magna was the same that displayed an obvious multiplication within A. castellanii and H. vermiformis associated to a pronounced cytolytic effect in these two species. Consequently, the mechanism of the defective replication and the absence of cytopathogenicity in W. magna lie in the amoeba itself rather than in the bacterial strain. Willaertia may reduce L. pneumophila growth in coculture either by inhibiting its intracellular replication and/or by blocking bacterial induced lysis. Another explanation may be the occurrence of differences in the membrane composition of amoebas that may result in a lower efficiency of bacterial internalization within W. magna. Some studies have shown heterogeneity in the attachment and uptake mechanisms of L. pneumophila by H. vermiformis and A. polyphaga [START_REF] Harb | Heterogeneity in the attachment and uptake mechanisms of the Legionnaires' disease bacterium, Legionella pneumophila, by protozoan hosts[END_REF].

Although none of the hypothesis may be excluded, our observations support the possibility that L. pneumophila replication is inhibited by W. magna.

Hence, a defective replication of the bacteria in W. magna amoebas clearly occurs as shown by the very low occurrence of replicative phagosomes in this species when compared to Acanthamoeba and Hartmannella amoebas. The microscopic observations suggest that the fusion of phagosome with lysosomes is not inhibited by L. pneumophila within Willaertia amoebas. Moreover, in the very few W. magna that displayed a replicative phagosome, this latter subcellular organelle appeared not as clearly delimited as the one in other amoebic species. To test the hypothesis of an efficient phago-lysosomal formation within Willaertia amoebas, colocalisation studies of L. pneumophila with late endosomal and lysosomal markers such as acid phosphatase or LAMP-1 glycoprotein could be performed.

The expression of replicative phagosomes in W. magna only occurred for the highest MOI (i.e 100), suggesting that the cell defence mechanisms of Willaertia amoebas may be partially overwhelmed at this bacterial concentration. Concomitantly, this observation demonstrates that ingestion of L. pneumophila occurs in W. magna. Furthermore, even when L. pneumophila exceptionally succeeded in this limited intracellular replication within W. magna, it displayed a particular aspect recalling the "cyst-like" or mature intracellular form described by Garduno in HELA cells [START_REF] Garduno | Intracellular growth of Legionella pneumophila gives rise to a differentiated form dissimilar to stationaryphase forms[END_REF]. The authors suggested that these "cyst-like" bacteria are resilient and metabolically dormant bacteria that showed to be more resistant to diverse chemical and physical stresses. It is thus tempting to speculate that the bacteria observed within Willaertia are individuals that display a particular ability to resist to the intracellular environment of this amoeba. Concomitantly, the observations showing that the L. pneumophila internalized within Willaertia displays a dormant like morphology similar to the one described by Garduno [START_REF] Garduno | Intracellular growth of Legionella pneumophila gives rise to a differentiated form dissimilar to stationaryphase forms[END_REF] strengthen our assumption that W. magna is an amoebic species having a particular resistance towards the bacteria, at least in our experimental conditions. Globally, our results show that at least one amoebic species may not be a good vector for the pathogenic bacteria. Our observations are in agreement with a recent publication of Declerck [START_REF] Declerck | Impact of non-Legionella bacteria on the uptake and intracellular replication of Legionella pneumophila in Acanthamoeba castellanii and Naegleria lovaniensis[END_REF] in which the authors shown that Naegleria lovaniensis is less permissive than A. castellanii to invasion and replication of L. pneumophila. This observation, reported for the first time for a Naegleria species, seems perfectly in line with our results since Naegleria and Willaertia are two related genera, which belong to the same family of Vahlkampfiidae (De [START_REF] Jonckheere | Willaertia magna gen. nov., sp. nov. (Vahlkampfiidae), a thermophilic amoeba found in different habitats[END_REF][START_REF] Page | The classification of "naked" amoebae (Phylum Rhizopoda)[END_REF][START_REF] Jonckheere | The phylogenetic position of the amoeboflagellate Willaertia deduced from SSUrDNA sequences[END_REF].

Indeed, these two genera share some characteristics that clearly differentiate them from the group of Hartmannella spp. and Acanthamoeba spp. An important difference between these amoebas is their mode of locomotion: Hartmannella and Acanthamoeba amoebas are moving by slow and progressive pseudopodial deformations while Willaertia and Naegleria amoebas do it by fast deformations with emission of eruptive lobopods.

Hartmannella and Acanthamoeba mitosis is similar to that of higher eucaryotic cells, while Willaertia and Naegleria division occurs through promitosis characterized by the persistence of the nuclear membrane and of the nucleolus.

The results reported herein demonstrate marked differences between the three tested amoebic species in their ability to serve as host cells for the replication of L. pneumophila.

Consequently, the distribution and the density of the pathogenic bacteria in the environment may be subject to variations according to the composition and the evolution of the amoebic population in the biotope with potential implications on the occurrence of outbreaks of Legionnaire's disease. Our results also suggest that W. magna may be a useful model to study the mechanisms of host resistance to L. pneumophila.

Table. 1. Effect of bacterial infection on the growth of amoebas

The different amoebic species were cocultured with L. pneumophila at an MOI of 50 and cells number determined as described in methods. The results are expressed as the number of amoebic cells/ml of medium and are the mean ± standard deviation of 6 (Hartmannella)

to 8 (

 8 Acanthamoeba and Willaertia) independent experiments. Statistical differences between the growth of Willaertia and other amoebic species are indicated (*: P<0.05; **: P<0.001). The number of Acanthamoeba and Hartmanella determined on day 3 of coculture represents a drop of ~ 93 and 82 % in the number of respective control amoebas grown in the same conditions for 3 days without L. pneumophilla. This drop was limited to ~ 25 % in the case of coculture experiments using W. magna.

Fig. 1 .

 1 Fig. 1. Growth of Legionella pneumophila in coculture with the three amoebic genera.

Fig. 2 .

 2 Fig. 2. Effect of L. pneumophila on amoebic monolayer formation.

Fig. 3 .

 3 Fig. 3. Cytotoxicity of L. pneumophila towards Acanthamoeba and Willaertia amoebas.

Fig. 4 .

 4 Fig. 4. Transmission electron micrographs of Acanthamoeba infected by L.

Fig. 5 .

 5 Fig. 5. Transmission electron micrographs of the rare Willaertia infected by L.

Fig. 6 .Fig

 6 Fig. 6. Sections of L. pneumophila in Acanthamoeba and Willaertia phagosomes.

Fig

  Fig. 2.

Table . 2 : Occurrence of replicative phagosomes within amoebas

 .

		Acanthamoeba	Hartmannella	Willaertia
		castellanii	vermiformis	magna
	MOI*= 50	36 %	42%	0 %
	MOI = 100	51 %	ND	1,5 %

* MOI : Multiplicity of infection

Amoebas were cocultured for 36 hours with L. pneumophila at a MOI of 50 or 100, fixed and processed for electron microscopy. The results are expressed as the percentage of cells displaying a replicative phagosome and were obtained out of the observation of more than 100 cells. The determination of phagosome occurrence at a MOI of 100 could not be performed in H. vermiformis on account of too many bacterial-induced lysed cells (ND: not determined).

Acknowledgments: This work was supported by a grant of the French "Direction Générale de la Santé" (ECOMICTH Association) and by the "Service des Etudes Médicales d'Electricité de France"