Francis Filbet

Luis Miguel Miguel

Miguel Luis

Rodrigues

Asymptotically stable particle-in-cell

Keywords: High order time discretization, Vlasov-Poisson system, Guiding-centre model, Particle methods. 1

Introduction

Magnetized plasmas are encountered in a wide variety of astrophysical situations, but also in magnetic fusion devices such as tokamaks, where a large external magnetic field needs to be applied in order to keep the particles on the desired tracks.

In Particle-In-Cell (PIC) simulations of such devices, this large external magnetic field obviously needs to be taken into account when pushing the particles. However, due to the magnitude of the concerned field this often adds a new time scale to the simulation and thus a stringent restriction on the time step. In order to get rid of this additional time scale, we would like to find approximate equations, where only the gross behavior implied by the external field would be retained and which could be used in a numerical simulation. In the simplest situation, the trajectory of a particle in a constant magnetic field B is a helicoid along the magnetic field lines with a radius proportional to the inverse of the magnitude of B. Hence, when this field becomes very large the particle gets trapped along the magnetic field lines. However due to the fast oscillations around the apparent trajectory, its apparent velocity is smaller than the actual one. This result has been known for some time as the guiding center approximation, and the link between the real and the apparent velocity is well known in terms of B.

Here, we consider a plasma constituted of a large number of charged particles, which is described by the Vlasov equation coupled with the Maxwell or Poisson equations to compute the self-consistent fields. It describes the evolution of a system of particles under the effects of external and self-consistent fields. The unknown f (t, x, v), depending on the time t, the position x, and the velocity v, represents the distribution of particles in phase space for each species with (x, v) ∈ R d × R d , d = 1, .., 3. Its behaviour is given by the Vlasov equation, (1.1) ∂f ∂t

+ v • ∇ x f + F(t, x, v) • ∇ v f = 0,
where the force field F (t, x, v) is coupled with the distribution function f giving a nonlinear system. We first define ρ(t, x) the charge density and J(t, x) the current density which are given by

ρ(t, x) = q R d f (t, x, v)dv, J(t, x) = q R d v f (t, x, v)dv,
where q is the elementary charge. For the Vlasov-Poisson model

(1.2) F(t, x, v) = q m E(t, x), E(t, x) = -∇ x φ(t, x), -∆ x φ = ρ ε 0 ,
where m represents the mass of one particle. On the other hand for the Vlasov-Maxwell model, we have

F(t, x, v) = q m (E(t, x) + v ∧ B(t, x) ),
and E, B are solutions of the Maxwell equations

                   ∂E ∂t -c 2 ∇ × B = - J ε 0 , ∂B ∂t + ∇ × E = 0, ∇ • E = ρ ε 0 , ∇ • B = 0,
with the compatibility condition ∂ρ ∂t + div x J = 0, which is verified by the solutions of the Vlasov equation.

Here we will consider an intermediate model where the magnetic field is given,

B(t, x) = 1 ε B ext (t, x ⊥ ),
with ε > 0 and we focus on the long time behavior of the plasma in the orthogonal plane to the external magnetic field, that is the two dimensional Vlasov-Poisson system with an external strong magnetic field

(1.3)        ε ∂f ∂t + v ⊥ • ∇ x ⊥ f + E(t, x ⊥ ) + 1 ε v ⊥ ∧ B ext (t, x ⊥ ) • ∇ v ⊥ f = 0, (x ⊥ , v ⊥ ) ∈ R 4 , E = -∇ x ⊥ φ, -∆ x ⊥ φ = ρ, x ⊥ ∈ R 2 .
Here, for simplicity we set all physical constants to one and consider that ε > 0 is a small parameter related to the ratio between the reciprocal Larmor frequency and the advection time scale. The term ε in front of the time derivative of f stands for the fact that we want to approximate the solution for large time. We want to construct numerical solutions to the Vlasov-Poisson system (1.3) by particle methods (see [START_REF] Birdsall | Plasma Physics via Computer Simulation Institute of Physics[END_REF]), which consist in approximating the distribution function by a finite number of macro-particles. The trajectories of these particles are computed from the characteristic curves corresponding to the Vlasov equation

(1.4)                ε dX dt = V, ε dV dt = 1 ε V ∧ B ext (t, X) + E(t, X), X(t 0 ) = x 0 ⊥ , V(t 0 ) = v 0 ⊥ ,
where the electric field is computed from a discretization of the Poisson equation in (1.3) on a mesh of the physical space.

The main purpose of this work is the construction of efficient numerical methods for stiff transport equations of type (1.3) in the limit ε → 0. Indeed, setting

Z = V ε - E ∧ B ext B ext 2 ,
the system (1.4) can be re-written for (X, Z) as

(1.5)                      dX dt = E ∧ B ext B ext 2 + Z, dZ dt = - 1 ε 2 B ext ∧ Z - d dt E ∧ B ext B ext 2 , X(t 0 ) = x 0 ⊥ , Z(t 0 ) = 1 ε v 0 ⊥ - E ∧ B ext B ext 2 (t 0 , x 0 ⊥ ).
Therefore, we denote by (X ε , Z ε ) the solution to (1.5), and under some classical smoothness assumptions on the electromagnetic fields (E, B ext ), it is well-known at least when v 0 ⊥ = 0 or B ext is homogeneous that (Z ε ) ε>0 converges weakly to zero when ε → 0, and X ε Y, where Y corresponds to the guiding center approximation

(1.6)        dY dt = E ∧ B ext B ext 2 (t, Y), Y(t 0 ) = x 0 ⊥ .
Here, we are of course interested in the behavior of the sequence solution (f ε ) ε>0 to the Vlasov-Poisson system (1.3) when ε → 0, which corresponds to the gyro-kinetic approximation of the Vlasov-Poissons system. Following the work of L. Saint-Raymond [START_REF] Saint-Raymond | The gyro-kinetic approximation for the Vlasov-Poisson system[END_REF], it can be proved -at least when B ext is homogeneous -that the charge density (ρ ε ) ε>0 converges to the solution to the guiding center approximation

(1.7)    ∂ρ ∂t + U • ∇ρ = 0, -∆ x ⊥ φ = ρ,
where the velocity U is

U = E ∧ B ext B ext 2 , E = -∇ x ⊥ φ.
We observe that the limit system (1.6) corresponds to the characteristic curves to the limit equation (1.7).

We seek a method that is able to capture these properties, while the numerical parameters may be kept independent of the stiffness degree of these scales. This concept is known and widely studied for dissipative systems in the framework of asymptotic preserving schemes [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multi scale kinetic equations[END_REF][START_REF] Klar | An asymptotic-induced scheme for non-stationary transport equations in the diffusive limit[END_REF]. Contrary to collisional kinetic equations in hydrodynamic or diffusion asymptotic, collisionless equations like the Vlasov-Poisson system (1.3) involve time oscillations. In this context, the situation is more complicated than the one encountered in collisional regimes since we cannot expect any dissipative phenomenon. Therefore, the notion of two-scale convergence has been introduced both at the theoretical and numerical level [START_REF] Crouseilles | Two-Scale Macro-Micro decomposition of the Vlasov equation with a strong magnetic field[END_REF][START_REF] Frénod | Long time simulation of a beam in a periodic focusing channel via a two-scale PIC-method[END_REF][START_REF] Frénod | Long time behavior of the Vlasov equation with strong external magnetic field[END_REF] in order to derive asymptotic models. However, these asymptotic models, obtained after removing the fast scales, are valid only when ε is small. We refer to E. Frénod and E. Sonnendrücker [START_REF] Frénod | Long time behavior of the Vlasov equation with strong external magnetic field[END_REF], and F. Golse and L. Saint-Raymond [START_REF] Golse | The Vlasov-Poisson system with strong magnetic field[END_REF][START_REF] Saint-Raymond | Control of large velocities in the two-dimensional gyro-kinetic approximation[END_REF] for a theoretical point of view on these questions, and E. Frénod, F. Salvarani and E. Sonnendrücker [START_REF] Frénod | Long time simulation of a beam in a periodic focusing channel via a two-scale PIC-method[END_REF] for numerical applications of such techniques.

Another approach is to combine both disparate scales into one and single model. Such a decomposition can be done using a micro-macro approach (see [START_REF] Crouseilles | Two-Scale Macro-Micro decomposition of the Vlasov equation with a strong magnetic field[END_REF] and the references therein). Such a model may be used when the small parameter of the equation is not everywhere small. Hence, a scheme for a micro-macro model can switch from one regime to another without any treatment of the transition between the different regimes. A different method consists in separating fast and slow time scales when such a structure can be identified [START_REF] Crouseilles | Asymptotic preserving schemes for highly oscillatory VlasovPoisson equations[END_REF] or [START_REF] Frénod | Long Time Behaviour of an Exponential Integrator for a Vlasov-Poisson System with Strong Magnetic Field[END_REF].

Theses techniques work well when the magnetic field is uniform since fast scales can be computed from a formal asymptotic analysis, but for more complicated problems, that is, when the external magnetic field depends on time and position x, the generalization of this approach is an open problem.

In this paper, we propose an alternative to such methods allowing to make direct simulations of systems (1.3) with large time steps with respect to ε. We develop numerical schemes that are able to deal with a wide range of values for ε, so-called Asymptotic Preserving (AP) class [START_REF] Klar | An asymptotic-induced scheme for non-stationary transport equations in the diffusive limit[END_REF][START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multi scale kinetic equations[END_REF], such schemes are consistent with the kinetic model for all positive value of ε, and degenerate into consistent schemes with the asymptotic model when ε → 0.

Before presenting our time discretization technique, let us first briefly review the basic tools of particle-in-cell methods which are widely used for plasma physics simulations [START_REF] Birdsall | Plasma Physics via Computer Simulation Institute of Physics[END_REF].

A brief review of particle methods

The numerical resolution of the Vlasov equation and related models is usually performed by Particle-In-Cell (PIC) methods which approximate the plasma by a finite number of particles. Trajectories of these particles are computed from characteristic curves (1.4) corresponding to the the Vlasov equation (1.3), whereas self-consistent fields are computed on a mesh of the physical space. This method yields satisfying results with a relatively small number of particles but it is sometimes subject to fluctuations, due to the numerical noise, which are difficult to control. To improve the accuracy, direct numerical simulation techniques have been developed. The Vlasov equation is discretized in phase space using either semi-Lagrangian [START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF][START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF][START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: theoretical analysis and application to the Vlasov-Poisson system[END_REF][START_REF] Sonnendrücker | The semi-Lagrangian method for the numerical resolution of Vlasov equation[END_REF], finite difference [START_REF] Yang | Conservative and non-conservative methods based on Hermite weighted essentially non-oscillatory reconstruction for Vlasov equations[END_REF] or discontinuous Galerkin [START_REF] De Dios | Discontinuous Galerkin methods for the multi-dimensional Vlasov-Poisson problem[END_REF][START_REF] Heath | A discontinuous Galerkin method for the Vlasov-Poisson system[END_REF] schemes. But these direct methods are very costly, hence several variants of particle methods have been developed over the past decades. In the Complex Particle Kinetic scheme introduced by Bateson and Hewett [START_REF] Bateson | Grid and Particle Hydrodynamics[END_REF][START_REF] Hewett | Fragmentation, merging, and internal dynamics for PIC simulation with finite size particles[END_REF], particles have a Gaussian shape that is transformed by the local shearing of the flow. Moreover they can be fragmented to probe for emerging features, and merged where fine particles are no longer needed. In the Cloud in Mesh (CM) scheme of Alard and Colombi [START_REF] Alard | A cloudy Vlasov solution[END_REF] particles also have Gaussian shapes, and they are deformed by local linearization of the force field. More recently in [START_REF] Pinto | Noiseless VlasovPoisson simulations with linearly transformed particles[END_REF], the authors proposed a Linearly-Transformed Particle-In-Cell method, that employs linear deformations of the particles.

Here we focus on the time discretization technique, hence we will only consider standard particle method even if our approach is completely independant on the choice of the particle method.

The particles method consists in approximating the initial condition f 0 in (1.3) by the following Dirac mass sum

f 0 N (x, v) := 1≤k≤N ω k δ(x -x 0 k ) δ(v -v 0 k ) ,
where (x 0 k , v 0 k ) 1≤k≤N is a beam of N particles distributed in the four dimensional phase space according to the density function f 0 . Afterwards, one approximates the solution of (1.3), by

f N (t, x, v) := 1≤k≤N ω k δ(x -X k (t)) δ(v -V k (t)) ,
where (X k , V k ) 1≤k≤N is the position in phase space of particle k moving along the characteristic curves (1.4) with the initial data (x 0 k , v 0 k ), for 1 ≤ k ≤ N . However when the Vlasov equation is coupled with the Poisson equation for the computation of the electric field, the Dirac mass has to be replaced by a smooth function ϕ α

f 0 N,α (x, v) := 1≤k≤N ω k ϕ α (x -x 0 k ) ϕ α (v -v 0 k ) ,
where ϕ ε = ε -d ϕ(•/ε) is a particle shape function with radius proportional to ε, usually seen as a smooth approximation of the Dirac measure obtained by scaling a compactly supported "cut-off" function ϕ for which common choices include B-splines and smoothing kernels with vanishing moments, see e.g. [START_REF] Koumoutsakos | Inviscid Axisymmetrization of an Elliptical Vortex[END_REF][START_REF] Cottet | Vortex Methods: Theory and Practice[END_REF]. Particle centers are then pushed forward at each time t by following a numerical approximation of the flow (1.4), leading to

f N,α (t, x, v) := 1≤k≤N ω k ϕ α (x -X k (t)) ϕ α (v -V k (t)) .
In the classical error analysis [START_REF] Beale | Vortex methods. II. Higher order accuracy in two and three dimensions[END_REF][START_REF] Raviart | An analysis of particle methods[END_REF], the above process is seen as

• An approximation (in the distribution sense) of the initial data by a collection of weighted Dirac measures ; • The exact transport of the Dirac particles along the flow ;

• The smoothing of the resulting distribution

1≤k≤N ω k δ(. -X k (t)) δ(. -V k (t))
with the convolution kernel ϕ ε . The classical error estimate reads then as follows [START_REF] Cohen | Optimal Approximations of Transport Equations by Particle and Pseudoparticle Methods[END_REF]: Proposition 2.1. Consider the Vlasov equation with a given electromagnetic field (E, B ext ) and a smooth initial datum f 0 ∈ C s c (R d ), with s ≥ 1. If for some prescribed integers m > 0 and r > 0, the cut-off ϕ ≥ 0 has m-th order smoothness and satisfies a moment condition of order r, namely,

R d ϕ(y) dy = 1, R d |y| r ϕ(y) dy < ∞,
and

R d y s 1 1 . . . y s d d ϕ(y) dy = 0, for s = (s 1 , . . . , s d ) ∈ IN d with 1 ≤ s 1 + • • • + s d ≤ r -1.
Then there exists a constant C independent of f 0 , N or α, such that we have for all 1 ≤ p ≤ +∞,

f (t) -f N,α (t) L p → 0,
when N → ∞ and α → 0 where the ratio N 1/d α 1.

Note that following [START_REF] Cohen | Optimal Approximations of Transport Equations by Particle and Pseudoparticle Methods[END_REF], it is also possible to get explicit order of convergence for the linear Vlasov equation. Let us also mention related papers where the convergence of a numerical scheme for the Vlasov-Poisson system is investigated. G.-H. Cottet and P.-A. Raviart [START_REF] Cottet | Particle methods for the one-dimensional Vlasov-Poisson equations[END_REF] present a precise mathematical analysis of the particle method for solving the one-dimensional Vlasov-Poisson system. We also mention the papers of S. Wollman and E. Ozizmir [START_REF] Wollman | Numerical approximation of the one-dimensional Vlasov-Poisson system with periodic boundary conditions[END_REF] and S. Wollman [START_REF] Wollman | On the approximation of the Vlasov-Poisson system by particle methods[END_REF] on the topic. K. Ganguly and H.D. Victory give a convergence result for the Vlasov-Maxwell system [START_REF] Ganguly | On the convergence of particle methods for multidimensional Vlasov-Poisson systems[END_REF].

The rest of the paper is organized as follows. In Section 3 we present several time discretization techniques based on high-order semi-implicit schemes [START_REF] Boscarino | High order semi-Implicit schemes for time dependent partial differential equations[END_REF] for the Vlasov-Poisson system with a strong external magnetic field, and we prove uniform consistency of the schemes in the limit ε → 0 with preservation of the order of accuracy (from first to third order accuracy). In Section 4 we perform a rigorous analysis of the first order scheme for smooth electromagnetic fields.

Section 5 is then devoted to numerical simulations for one single particle motion and for the Vlasov-Poisson model for various asymptotics ε ≈ 1 and ε 1, which illustrate the advantage of high order schemes.

A particle method for Vlasov-Poisson system with a strong magnetic field

Let us now consider the system (1.3) and apply a particle method, where the key issue is to design a uniformly stable scheme with respect to the parameter ε > 0, which is related to the magnitude of the external magnetic field. Assume that at time t n = n ∆t, the set of particles are located in (x n k , v n k ) 1≤k≤N , we want to solve the following system on the time interval [t n , t n+1 ],

(3.1)

               ε dX k dt = V k , ε dV k dt = 1 ε V k ∧ B ext (t, X k ) + E(t, X k ), X(t n ) = x n k , V(t n ) = v n k ,
where the electric field is computed from a discretization of the Poisson equation (1.3) on a mesh of the physical space.

The numerical scheme that we describe here is proposed in the framework of Particle-In-Cell method, where the solution f is discretized as follows

f n+1 N,α (x, v) := 1≤k≤N ω k ϕ α (x -x n+1 k ) ϕ α (v -v n+1 k ),
where

(x n+1 k , v n+1 k
) represents an approximation of the solution X k (t n+1 ) and V k (t n+1 ) to (3.1). When the Vlasov equation (1.3) is coupled with the Poisson equation (1.2), the electric field is computed in a macro-particle position x n+1 k at time t n+1 as follows

• Compute the density ρ ρ n h,ε (x) = k∈Z Z d w k ϕ ε (x -x n k ).
• Solve a discrete approximation to (1.2)

-∆ h φ n (x) = ρ n h,ε (x).
• Interpolate the electric field with the same order of accuracy on the points (x n k ) k∈Z Z d . To discretize the system (3.1), we apply the strategy developed in [START_REF] Boscarino | High order semi-Implicit schemes for time dependent partial differential equations[END_REF] based on semi-implicit solver for stiff problems. In the rest of this section, we propose several numerical schemes to the system (3.1) for which the index k ∈ {1, . . . , N } will be omitted.

3.1.

A first order semi-implicit scheme. We start with the simplest semi-implicit scheme for (3.1), which is a combination of backward and forward Euler scheme. It gives for a fixed time step ∆t > 0 and a given electric field E and an external magnetic field B ext ,

(3.2)          x n+1 -x n ∆t = v n+1 ε , v n+1 -v n ∆t = 1 ε v n+1 ε ∧ B ext (t n , x n ) + E(t n , x n ) .
Note that only the second equation on v n+1 is fully implicit and requires the inversion of a linear operator. Then, from v n+1 the first equation gives the value for the position x n+1 .

Proposition 3.1 (Consistency in the limit ε → 0 for a fixed ∆t). Let us consider a time step ∆t > 0, a final time T > 0 and set

N T = [T /∆t]. Assume that the sequence (x n ε , v n ε ) 0≤n≤N T given by (3.2) is such that for all 1 ≤ n ≤ N T , (x n ε , εv n ε ) ε>0
is uniformly bounded with respect to ε > 0 and (x 0 ε , εv 0 ε ) ε>0 converges in the limit ε → 0 to some (y 0 , 0). Then, for 1 ≤ n ≤ N T , x n ε → y n , as ε → 0 and the limit (y n ) 1≤n≤N T is a consistent first order approximation with respect to ∆t of the guiding center equation provided by the scheme

(3.3) y n+1 -y n ∆t = E(t n , y n ) ∧ B ext (t n , y n ) B 2 .
Proof. For all 1 ≤ n ≤ N T , we consider (x n ε , v n ε ) the solution to (3.2) now labeled with respect to ε > 0. Since, the sequence (x n ε ) ε>0 is uniformly bounded with respect to ε > 0, we can extract a subsequence still abusively labeled by ε and find some (y n ) such that x n ε → y n as ε goes to zero. Then, we observe that the second equation of (3.2) can be written as

ε 2 v n+1 ε -v n ε ε∆t = v n+1 ε ε ∧ B ext (t n , x n ε ) + E(t n , x n ε ) .
and that, for any 0 ≤ n ≤ N T , (εv n ε ) ε>0 is uniformly bounded. From this we conclude first that for any 1 ≤ n ≤ N T , (ε -1 v n ε ) ε>0 is uniformly bounded then that εv n ε → 0 for any 0 ≤ n ≤ N T . Therefore, taking the limit ε → 0, it yields that for 0

≤ n ≤ N T -1, ε -1 v n+1 ε → E(t n , y n ) ∧ B ext (t n , y n ) B ext 2 , when ε → 0.
Substituting the limit of ε -1 v n+1 in the first equation of (3.2) we prove that the limit y n satisfies (3.3). Since the limit point y n is uniquely determined, actually all the sequence (x n ε ) ε>0 converges.

Remark 3.2. The consistency provided by the latter result is far from being uniform with respect to the time step. However we do prove in the next section that the solution to (3.2) is both uniformly stable and consistent with respect to ∆t and ε > 0.

Of course, such a first order scheme is not accurate enough to describe correctly the long time behavior of the numerical solution, but it has the advantage of the simplicity and we will prove in the next section that it is uniformly stable with respect to the parameter ε and the sequence (x n ε ) converges to a consistent approximation of the guiding center model when ε → 0. Now, let us see how to generalize such an approach to second and third order schemes.

3.2.

Second order semi-implicit Runge-Kutta schemes. Now, we consider second order schemes with two stages.

3.2.1.

A second order A-stable scheme. A first example of scheme satisfying the second order conditions is given by a combination of Heun method (explicit part) and an A-stable second order singly diagonal implicit Runge-Kutta SDIRK method (implicit part) [START_REF] Hairer | Solving ordinary differential equations. I. Non-stiff problems[END_REF][START_REF] Boscarino | High order semi-Implicit schemes for time dependent partial differential equations[END_REF]. The first stage corresponds to 1) ,

(3.4)          x (1) = x n + ∆t 2ε v (
v (1) = v n + ∆t 2ε v (1) ε ∧ B ext (t n , x n ) + E(t n , x n ) .
Then the second stage is given by

(3.5)          x (2) = x n + ∆t 2ε v (2) , v (2) = v n + ∆t 2ε v (2) ε ∧ B ext (t n+1 , 2x (1) -x n ) + E(t n+1 , 2x (1) -x n ) .
Finally, the numerical solution at the new time step is

(3.6)    x n+1 = x (1) + x (2) -x n , v n+1 = v (1) + v (2) -v n .
A similar numerical scheme has been proposed in the framework of δf simulation of the Vlasov-Poisson system [START_REF] Cheng | A second-order semi-implicit df method for hybrid simulation[END_REF].

Under stability assumptions on the numerical solution to (3.4)-(3.6), we get the following consistency result in the limit ε → 0. Proposition 3.3 (Consistency in the limit ε → 0 for a fixed ∆t). Let us consider a time step ∆t > 0, a final time T > 0 and set N T = [T /∆t]. Assume that the sequence (x n ε , v n ε ) 0≤n≤N T given by (3.4)-(3.6) is such that for all 1 ≤ n ≤ N T , (x n ε , εv n ε ) ε>0 is uniformly bounded with respect to ε > 0 and (x 0 ε , εv 0 ε ) ε>0 converges in the limit ε → 0 to some (y 0 , 0). Then, for 1 ≤ n ≤ N T , x n ε → y n , as ε → 0 and the limit (y n ) n≥1 is a consistent and second order approximation with respect to ∆t of the guiding center equation given by the scheme

(3.7)            y (1) = y n + ∆t 2 E(t n , y n ) ∧ B ext (t n , y n ) B ext 2 , y (2) = y n + ∆t 2 E(t n+1 , 2y (1) -y n ) ∧ B ext (t n+1 , 2y (1) -y n ) B ext 2 ,
and y n+1 = y (1) + y (2) -y n .

Proof. We follow the lines of the proof of Proposition 3.1 and first choose a subsequence of (x n ε ) converging to some y n . The second equation of (3.4) implies that ε -1 v

(1) ε is bounded. From the first equation of (3.4) it follows that so is x ε is also bounded. Now from the second equation of (3.6) stems that εv n ε converges to zero as ε → 0 for any 0 ≤ n ≤ N T . Coming back to the second equation of (3.4) we show that for all n, ε -1 v

(1)

ε converges when ε → 0 and v (1) ε ε → E(t n , y n ) ∧ B ext (t n , y n ) B ext 2 , when ε → 0.
Using the first equation of (3.4) we conclude that x (1) also converges and that its limit y (1) is given by the first equation of (3.7). Going on with the same arguments shows that

v (2) ε ε → E(t n+1 , 2y (1) -y n ) ∧ B ext (t n+1 , 2y (1) -y n ) B ext 2
, when ε → 0 and that x (2) converges to a y (2) given by the second equation of (3.7). One may then take a limit in the first equation of (3.6) and conclude that indeed the limit y n satisfies (3.7). Again uniqueness supplies the convergence of the whole sequence.

3.2.2.

A second order L-stable scheme. Another choice is a combination of Runge-Kutta method (explicit part) and an L-stable second order SDIRK method in the implicit part. This implicit scheme should give better stability properties on the numerical solution with respect to the stiffness parameter ε > 0.

We first choose γ > 0 as the smallest root of the polynomial γ 2 -2γ + 1/2 = 0, i.e. γ = 1 -1/ √ 2, then the scheme is given by the following two stages. First, we have 1) , 1) , with

(3.8)          x (1) = x n + γ∆t ε v (
v (1) = v n + γ∆t ε F (
F (1) := v (1) ε ∧ B ext (t n , x n ) + E(t n , x n ).
For the second stage, we first define 1) , then the solution (x (2) , v (2) ) is given by (3.10)

(3.9) t(1) := t n + ∆t 2γ , x(1) := x n + ∆t 2γε v (
         x (2) = x n + (1 -γ)∆t ε v (1) + γ∆t ε v (2) , v (2) = v n + (1 -γ)∆t ε F (1) + γ∆t ε F (2) , with
F (2) := v (2) ε ∧ B ext t(1) , x(1) + E t(1) , x (1)
.

Finally, the numerical solution at the new time step is (2) .

(3.11)    x n+1 = x (2) , v n+1 = v
Under stability assumptions on the numerical solution to (3.4)-(3.6), we get the following consistency result in the limit ε → 0. Proposition 3.4 (Consistency in the limit ε → 0 for a fixed ∆t). Let us consider a time step ∆t > 0, a final time T > 0 and set N T = [T /∆t]. Assume that the sequence (x n ε , v n ε ) 0≤n≤N T given by (3.8)-(3.11) is such that for all 1 ≤ n ≤ N T , (x n ε , εv n ε ) ε>0 is uniformly bounded with respect to ε > 0 and (x 0 ε , εv 0 ε ) ε>0 converges in the limit ε → 0 to some (y 0 , 0). Then, for 1 ≤ n ≤ N T , x n ε → y n , as ε → 0 and the limit (y n ) n≥1 is a consistent second order approximation of the guiding center equation, given by (3.12) 1) , where

       U n = E(t n , y n ) ∧ B ext (t n , y n ) B ext 2 , y n+1 = y n + (1 -γ)∆t U n + γ∆t U (
ŷ(1) := y n + ∆t 2γ U n , U (1) := E( t(1) , ŷ(1) ) ∧ B ext ( t(1) , ŷ(1) ) B ext 2 .
We omit the proof of Proposition 3.4 as almost identical to the one of Proposition 3.3. The present scheme is L-stable, which means uniformly linearly stable with respect to ∆t.

3.3.

Third order semi-implicit Runge-Kutta schemes. A third order semi-implicit scheme is given by a four stages Runge-Kutta method introduced in the framework of hyperbolic systems with stiff source terms [START_REF] Boscarino | High order semi-Implicit schemes for time dependent partial differential equations[END_REF]. First, we set α = 0.24169426078821, β = α/4 and η = 0.12915286960590 and γ = 1/2 -α -β -η. Then we construct the first stage as 1) , with

(3.13)          x (1) = x n + α∆t ε v (1) , v (1) = v n + α∆t ε F (
F (1) := v (1) ε ∧ B ext (t n , x n ) + E(t n , x n ).
For the second stage, we have

(3.14)          x (2) = x n - α∆t ε v (1) + α∆t ε v (2) , v (2) = v n - α∆t ε F (1) + α∆t ε F (2) , with
F (2) := v (2) ε ∧ B ext (t n , x n ) + E (t n , x n ) .
Then, for the third stage we set 2) .

(3.15)          x (3) = x n + (1 -α)∆t ε v (2) + α∆t ε v (3) , v (3) = v n + (1 -α)∆t ε F (2) + α∆t ε F (3) , with          F (3) := v (3) ε ∧ B ext t n+1 , x (2) + E t n+1 , x (2) , x (2)
:= x n + ∆t ε v (
Finally, for the fourth stage we set (3.16)

         x (4) = x n + β∆t ε v (1) + η∆t ε v (2) + γ∆t ε v (3) + α∆t ε v (4) , v (4) = v n + β∆t ε F (1) + η∆t ε F (2) + γ∆t ε F (3) + α∆t ε F (4) , with        F (4) := v (4) ε ∧ B ext t n+1/2 , x (3) + E t n+1/2 , x (3) , x (3)
:= x n + ∆t 4ε v (2) + v (3) , and the numerical solution at the new time step is (3.17) 4) , 4) . As for the previous schemes, under uniform stability assumptions with respect to ε > 0, we prove the following Proposition Proposition 3.5 (Consistency in the limit ε → 0 for a fixed ∆t). Let us consider a time step ∆t > 0, a final time T > 0 and set N T = [T /∆t]. Assume that the sequence (x n ε , v n ε ) 0≤n≤N T given by (3.13)-(3.17) is such that for all 1 ≤ n ≤ N T , (x n ε , εv n ε ) ε>0 is uniformly bounded with respect to ε > 0 and (x 0 ε , εv 0 ε ) ε>0 converges in the limit ε → 0 to some (y 0 , 0). Then, for 1 ≤ n ≤ N T , x n ε → y n , as ε → 0 and the limit (y n ) n≥1 is a consistent third order approximation of the guiding center equation provided by the scheme 2) + U (3) , and 4) .

   x n+1 = x n + ∆t 6ε v (2) + v (3) + 4 v (
v n+1 = v n + ∆t 6ε F (2) + F (3) + 4 F (
(3.18)                        U (2) = E(t n , y n ) ∧ B ext (t n , y n ) B ext 2 , U (3) = E(t n+1 , y (2) ) ∧ B ext (t n+1 , y (2) ) B ext 2 , U (4) = E(t n+1/2 , y (3) ) ∧ B ext (t n+1/2 , y (3) ) B ext 2 , where      y (2) = y n + ∆t U (2) , y (3) = y n + ∆t 4 U (
y n+1 = y n + ∆t 6 U (2) + U (3) + 4 U (
Proof. We follow the lines of the proof of Proposition 3.3. In particular we prove successively the boundedness of ε -1 v

(1)

ε , x (1)
ε , ε -1 v (2) ε , x (2)
ε , ε -1 v (3) ε , x (3)
ε , ε -1 v (4) ε , x (4)
ε . Arguing inductively on n one then shows that εv n ε , F (1) , F (2) , F (3) , F (4) all converge to zero when ε → 0. Then for any converging subsequence of x n ε we may successively identify limits for ε -1 v

(1)

ε , x (1)
ε , ε -1 v (2) ε , x (2) ε , ε -1 v (3) ε , x (3) ε , ε -1 v (4) ε , x (4)
ε and prove that any accumulating point of x n ε solves the limiting scheme involving (3.18). Hence the result.

Analysis of the first-order semi-implicit scheme

Consider the first order Euler semi-implicit scheme (4.1)

                 ε x n+1 k -x n k ∆t = v n+1 k , ε v n+1 k -v n k ∆t = 1 ε v n+1 k ∧ B ext (t n , x n k ) + E(t n , x n k ), x 0 k = x 0 k , v 0 k = v 0 k .
We focus here on the case where both E and B ext are given external fields. Since, in this case, there is no coupling between particles we may focus safely on one of them and drop the k suffix. Then we prove Theorem 4.1 (Uniform consistency with respect to ε). We set

λ := ∆t ε 2 , assume that the (E, B ext ) ∈ W 1,∞ loc (R + × R 2
) and consider the solution (x n , v n ) to (4.1). Then there exist positive constants C and λ 0 , only depending on E, B ext and t n such that when 1 0 ≤ ∆t ≤ 1 and λ ≥ λ 0

x n -y n ≤ C∆t λ 1 + ε -1 v 0 -R[E(t 0 , x 0 )] ,
1 Obviously one may replace 1 by any upper bound on ∆t but λ 0 and C would depend on this upper bound.

where (y n ) corresponds to the numerical solution of the guiding center model

(4.2)        y n+1 -y n ∆t = E(t n , y n ) ∧ B ext B ext 2 , y 0 = x 0 .
Proof. To start with we analyze the case where E is a given bounded Lipschitz field and B ext is constant and of norm 1. Our estimates shall be expressed in terms of

K 0 = E L ∞ K t = ∂E ∂t L ∞ K x = ∇ x E L ∞ .
For concision's sake we also introduce the operator to denote

R[W] = W ∧ B ext .
By introducing the key quantity for n ≥ 1,

z n = ε -1 v n -R[E(t n-1 , x n-1 )],
the first equation of (4.1) reads

x n = x n-1 + ∆t R[E(t n-1 , x n-1 )] + ∆t z n , n ≥ 1 while z 1 = [Id -λR] -1 (ε -1 v 0 -R[E(t 0 , x 0 )]
) and

z n+1 = [Id -λR] -1 (z n -R[E(t n , x n ) -E(t n-1 , x n-1 )]) , n ≥ 1 .
We first observe that R ≤ 1 and [

Id -λR] -1 = (Id + λR)/(1 + λ 2 ) since R 2 = -Id.
This leads to

z n+1 ≤ 1 + λ 1 + λ 2 [(1 + K x ∆t) z n + (K t + K x K 0 )∆t] , n ≥ 1 . Assuming 1 + λ 1 + λ 2 (1 + K x ∆t) < 1 and introducing a = 1 + λ 1 + λ 2 K t + K x K 0 1 -1+λ 1+λ 2 (1 + K x ∆t) ∆t we infer z n ≤ a + 1 + λ 1 + λ 2 (1 + K x ∆t) n-1 b , n ≥ 1 , where b = a + 1 + λ 1 + λ 2 1 ε v 0 -R[E(t 0 , x 0 )] .
For comparison we define y n solution to (4.2). Then,

x n -y n ≤ (1 + K x ∆t) x n-1 -y n-1 + ∆t a + 1 + λ 1 + λ 2 (1 + K x ∆t) n-1 b , n ≥ 1 .
Hence assuming moreover K x > 0 (or replacing K x with some arbitrary positive number if

∇ x E L ∞ = 0) x n -y n ≤ a K x [ 1 + (1 + K x ∆t) n ] + b ∆t (1 + K x ∆t)[1 -1+λ 1+λ 2 ] 1 + 1 + λ 1 + λ 2 n (1 + K x ∆t) n .
As a conclusion, for any 0 ≤ θ < 1, there exists C θ such that if

1 + λ 1 + λ 2 (1 + K x ∆t) ≤ θ then, it yields (4.3) x n -y n ≤ C θ ∆t λ K 0 + K t K x + 1 ε v 0 -R[E(t 0 , x 0 )] e Kx n∆t ,
which concludes the proof of Theorem 4.1 when the magnetic field is uniform. We now relax the assumption that B ext is constant and of norm 1. We set b ext = B ext and let R be dependent on t and x. We observe that now

R 2 = -b 2 ext Id so that [Id -λR] -1 = (Id + λR)/(1 + λ 2 b 2 ext ).
Then introducing the drift force

F(t, x) = 1 B ext (t, x) 2 E(t, x) ∧ B ext (t, x)
we essentially obtain the same estimates with 1+λ 1+λ 2 replaced by 1

Λ := 1 + λb ext 1 + λ 2 b 2 ext L ∞ and K 0 = F L ∞ K t = ∂F ∂t L ∞ K x = ∇ x F L ∞ .
Indeed introducing for n ≥ 1,

z n = v n ε -F(t n-1 , x n-1 ),
the scheme is written

z 1 = [Id -λR(t 0 , x 0 )] -1 v 0 ε -F(t 0 , x 0 ) ,
and then

z n+1 = [Id -λR(t n , x n )] -1 (z n -(F(t n , x n ) -F(t n-1 , x n-1 ))) , n ≥ 1 together with x n = x n-1 + ∆t F(t n-1 , x n-1 ) + ∆t z n , n ≥ 1 . Mark that if λ × (inf B ext ) > √ 2 -1 then Λ ≥ 1 + λ 2 × (inf B ext ) 2 1 + λ × (inf B ext ) .
Note that this result can be slightly improved when we modify the initial condition of the asymptotic discrete model. Indeed, consider y n solving (4.4)

     y n+1 -y n ∆t = E(t n , y n ) ∧ B ext ,
y 0 = x 0 + ε (v 0 ∧ B ext + ε E(t 0 , x 0 )) .
The gain is that now

x 1 -y 1 ≤ ∆t λ K x ∆t + 1 + λ 1 + λ 2 1 ε v 0 -R[E(t 0 , x 0 )] since z 1 = 1 λ R 1 ε v 0 -R[E(t 0 , x 0 )] - 1 λ R [Id -λR] -1 1 ε v 0 -R[E(t 0 , x 0 )] .
This leads, for n ≥ 1 to

x n -y n ≤ a K x 1 + (1 + K x ∆t) n-1 + b ∆t 1 -1+λ 1+λ 2 1 + λ 1 + λ 2 1 + 1 + λ 1 + λ 2 n-1 (1 + K x ∆t) n-1 + ∆t λ K x ∆t + 1 + λ 1 + λ 2 1 ε v 0 -R[E(t 0 , x 0 )] (1 + K x ∆t) n-1
then with notation as above there exists a constant

C = C θ (K 0 + Kt Kx + K x + 1)e Kxt n > 0 such that (4.5) x n -y n ≤ C ∆t λ 1 + 1 λ + ∆t 1 ε v 0 -R[E(t 0 , x 0 )] .
Concerning the analysis of high-order schemes, it is not straightforward to adapt directly the strategy of Theorem 4.1. Indeed, the use of a semi-implicit scheme does not necessarily guarantee that the particle trajectories are under control. Remark 4.2. Consider the scheme (3.4)- (3.6) in the simplest situation where the electric field is zero and the magnetic field is B ext = (0, 0, 1), we show that the scheme preserves the kinetic energy and we have

v n 2 = v 0 2 , ∀ n ∈ IN,
hence as ε goes to zero, the velocity v n /ε cannot converge to the null guiding center velocity except if the initial velocity does. Fortunately, for the other high order schemes (3.8)-(3.11) and (3.13)-(3.17), kinetic energy is dissipated and converges to 0.

Numerical simulations

In this section, we discuss some examples to validate and to compare the different time discretization schemes. We first consider the single motion of a particle under the effect of a given electromagnetic field. It allows us to illustrate the ability of the semi-implicit schemes to capture the guiding center velocity with large time step ∆t in the limit ε → 0.

Then we consider the Vlasov-Poisson system with an external magnetic field. A classical particle-in-cell method is applied with different time discretization techniques to compute the particle trajectories. Hence this collection of charged particles move and give rise to a selfconsistent electric field, obtained by solving numerically the Poisson equation in a space grid. 5.1. One single particle motion. Before going to the statistical descriptions, let us investigate the accuracy and stability properties of the semi-implicit algorithms presented in Section 3 on the motion of individual particles in a given electromagnetic field.

Here we consider an electric field E = -∇φ, where

φ(x) = 1 2 x 2 + α cos 2 (2π y) , x = (x, y) ∈ R 2 ,
with α = 0.02 and a magnetic field B(x) = 1 + 10 -1 sin(2πx) with x = (x, y) ∈ R 2 . We choose for all simulations ∆t = 0.1 and the initial data as x 0 = (1, 1.4) and ε -1 v 0 = (3, 5), such that the initial data z 0 = ε -1 v 0 + F(0, x 0 ) is bounded with respect to ε > 0. Thus, we apply the schemes proposed in Section 3 to compute a numerical solution to (3.1) and present the particle trajectory and velocity in Figures 1 and2. These results are compared with those obtained with a fourth order Runge-Kutta scheme using a small time step.

In Figure 1, we first investigate the case where ε is of order one (ε = 1 and 0.5), which is the non stiff regime. On the left column, we clearly observe that the space trajectory obtained from high order schemes agree very well with the reference trajectory, whereas after few time steps the first order scheme does not give satisfying results. On the right hand side, we present the evolution of the velocity at each time step and compare it with the reference velocity and the guiding center velocity. In this non stiff regime, the velocity obtained from high order schemes coincides with the reference velocity and the guiding center velocity is meaningless. Therefore, this first test illustrates the ability of high order schemes to describe accurately the particle motion in phase space when ε ∼ 1.

In Figure 2, we now propose the numerical results when ε 1, that is, ε = 0.1 and ε = 0.01, which corresponds to the high field regime. In that case, the space trajectory in the orthogonal plane to the magnetic field can be decomposed into a relatively slow motion due to the guiding center velocity

F(t, x) = 1 B ext (t, x) 2 E(t, x) ∧ B ext (t, x)
along the field line and a fast circular motion with a frequency of order 1/ε. Our aim here is to capture the slow motion using a fixed time step ∆t = 0.1 independently of the value ε 1. On the one hand, the space trajectory (left column of Figure 2) of the numerical solution remains stable for various ε > 0 even if we do not solve the fast scales. Moreover, when ε → 0, the numerical solution approaches the correct trajectory and fast fluctuations are somehow filtered thanks to the implicit treatment of the velocity v n . As in the previous simulations, we observe a discrepancy between the first order scheme and other high order schemes for large time.

Furthermore, we focus on the velocity variable (v n ) 0≤n≤N T and compare its time evolution with the reference velocity and the guiding center velocity (right column of Figure 2). Since we use a large time step ∆t = 0.1, we cannot expect to follow all the details of the velocity variable, but only the slow motion corresponding to the guiding center velocity. Now we clearly observe a different behavior of the numerical solutions. On the one hand, the scheme (3.4)- (3.6), which preserves the kinetic energy 1 2 |v n | 2 when E = 0 and B ext = (0, 0, 1) (see Remark 4.2), still oscillates with a large amplitude and the velocity does not coincide with the guiding center velocity. On the other hand, the schemes (3.2), (3.8)-(3.11) and (3.13)-(3.17) are more dissipative and after few time steps the velocity follows the line of the guiding center velocity (see the zoom on the right column of Figure 2). Hence the amplitude of oscillations in the physical space are diminishing and the particle follows the trajectory corresponding to the guiding center model (1.7).

As a conclusion, these elementary numerical simulations confirm the ability of the semiimplicit discretization to capture the slow motion corresponding to the guiding center model (1.7) uniformly with respect to ε 1 and the interest of high order time discretization for the long time behavior of the solution. 5.2. Diocotron instability. We now consider the diocotron instability [START_REF] Davidson | Physics of Nonneutral Plasmas[END_REF] for an annular electron layer usually described by the guiding center model (1.7). This instability is well studied numerically as in [START_REF] Yang | Conservative and non-conservative methods based on Hermite weighted essentially non-oscillatory reconstruction for Vlasov equations[END_REF][START_REF] Pétri | Non-linear evolution of the diocotron instability in a pulsar electrosphere: 2D PIC simulations[END_REF]. It may give rise to electron vortices, which is the analog of the Kelvin-Helmholtz fluid dynamics and may occur when charge neutrality is not locally maintained.

Here we want to investigate the development of such instability when we consider the Vlasov-Poisson system with an external magnetic field (1.3). It is expected that such instability holds true for large magnetic fields, whereas dissipative effects dominate when the magnetic field is not large enough to confine particles. Therefore, we perform numerical simulation using a particle-in-cell method [START_REF] Birdsall | Plasma Physics via Computer Simulation Institute of Physics[END_REF] for the Vlasov equation (1.3), where the particle trajectories are approximated thanks to the third order semiimplicit scheme (3.13)-(3.17). On the other hand, we compute an approximation of the guiding center model using a finite difference method developed in [START_REF] Yang | Conservative and non-conservative methods based on Hermite weighted essentially non-oscillatory reconstruction for Vlasov equations[END_REF]. This reference solution will be used to compare our results obtained from the Vlasov-Poisson system with a large magnetic fields.

The initial density is given by

ρ 0 (x) = (1 + α cos( θ)) exp (-4( x -6.5) 2 ), if r -≤ x ≤ r + , 0, otherwise,
where α is a small parameter, θ = atan(y/x), with x = (x, y) ∈ R 2 . In the following tests, we take α = 0.01, r -= 5, r + = 8, = 7. Furthermore, we will also consider the Vlasov-Poisson system (1.3) with an external magnetic fiel with the initial datum f 0

f 0 (x, v) = ρ 0 (x) 2π exp - v 2 2 , (x, v) ∈ R 4 .
A particle method with a third order semi-implicit solver (3.13))-(3.17) will be applied for different values of ε = 1, 10 -1 and 10 -2 . For both systems a high order finite difference scheme in Cartesian coordinates will be applied [START_REF] Yang | Conservative and non-conservative methods based on Hermite weighted essentially non-oscillatory reconstruction for Vlasov equations[END_REF] for the numerical approximation of the Poisson equation in a disc domain. We choose a grid with 100 × 100 points in the physical space and 400 particles per cell for the discretization of the velocity space.

We define the total energy E(t) as the sum of the potential energy E pot (t) and the kinetic energy E kin (t) with

E pot (t) = 1 2 R 2 |∇φ(t, x)| 2 dx, and E kin (t) = 1 2 R 4 f (t, x, v) |v| 2 dxdv.
For the Vlasov-Poisson system (1.3), the total energy is exactly conserved with respect to time. First, we consider the case where ε = 1, that is, the particle trajectories do not coincide with the trajectory corresponding to the guiding center model (1.7). On the one hand, we present in Figures 3, the evolution of relative energy with respect to the initial data t → E(t) -E(0) and the L ∞ norm of the self-consistent electric field t → E(t) ∞ . The discrete total energy is not exactly preserved but it oscillates around zero and these variations remain relatively small compared to the variations of the kinetic and potential energy which also oscillate with a frequency around 1/2π. The same phenomenon can be observed on the time evolution of the L ∞ norm of the electric field. On the other hand in Figure 4, we plot the time evolution of the charged density and observe that when the amplitude of the external magnetic field B is of order one, the plasma is not well confined and does not develop any instability. The density seems to oscillate around the steady state. Then, we take a small value ε = 0.01, and for this case we present a comparison between the finite difference approximation to the guiding center model (1.7) and a particle method for the Vlasov-Poisson system (1.3), using the third order semi-implicit scheme (3.13)-(3.17). The numerical results are presented in Figures 5 and6.

Let us emphasize that for the guiding center model (1.7), that is in the limit ε → 0, the potential energy E pot is conserved with respect to time. Therefore, for ε 1, it is expected that both the variations of E pot and E kin are small. In Figure 5, we can indeed observe that the variation of both quantities E pot and E kin are varying with a small amplitude of order 10 -3 . In that case, we do not see any oscillations since a large time step is used. Moreover, the second picture in Figure 5 represents the time evolution of E(•) ∞ for the Vlasov-Poisson system (1.3) and the guiding center model (1.7). Both results agree well which illustrates the accuracy of the particle method in the limit ε 1. Finally in Figure 6, we present a comparison between the density ρ obtained from the Vlasov-Poisson system discretized with the particle method and the one corresponding to the guiding center model discretized with a finite difference method. These figures show the development of the diocotron instability on the density ρ for both models. Indeed in this regime (ε 1), it is expected that the density ρ ε computed from the Vlasov-Poisson system obeys to the same evolution as the one of the guiding center model. Once again, the results agree very well and it is remarkable that the particle-in cell method does not suffer from too many fluctuations. The vortices are well described even for large time t ∼ 120.

Conclusion and perspective

In this paper we proposed a class of semi-implicit time discretization techniques for particlein cell simulations. The main feature of this approach is to guarantee the accuracy and stability when the amplitude of the magnetic field becomes large and to get the correct long time behavior (guiding center approximation). We formally showed that the present schemes preserve the initial order of accuracy when ε → 0. Furthermore, we performed a complete analysis of the first order semi-implicit scheme when we consider a given and smooth electromagnetic field (E, B ext ).

The time discretization techniques proposed in this paper seem to be a very simple and efficient tool to filter fast oscillations and have nice stability and consistency properties in the limit ε → 0. However, a complete analysis of high order semi-implicit schemes is still missing. The main issue is to control the space trajectory (x n ) n uniformly with respect to ε and as we have shown in Remark 4.2, the use of a semi-implicit scheme does not necessarily guarantee that the particle trajectories are under control. A complete analysis of high order schemes is currently under study.

On the other hand, the present techniques will be applied to more advanced problems as the three dimensional Vlasov-Poisson system when the magnetic field is non uniform and the particle trajectories become more complicated.

Luis Miguel Rodrigues

Université de Rennes 1, UMR6625, IRMAR, 263 avenue du General Leclerc, F-35042 Rennes Cedex, FRANCE e-mail: luis-miguel.rodrigues@univ-rennes1.fr

( 1 )

 1 ε . Then the second equation of (3.5) yields boundedness of ε-1 v (2)ε and the first that x[START_REF] Bateson | Grid and Particle Hydrodynamics[END_REF]

Figure 1 .

 1 Figure 1. One single particle motion. Numerical solution obtained with a large time step ∆t = 0.1 with the 1st order scheme (3.2), the 2nd order schemes (3.4)-(3.6) and (3.8)-(3.11) and the third order scheme (3.13)-(3.17) for ε = 1 (top) and ε = 5. 10 -1 (bottom) : (a) particle trajectory in physical space (x n ) 0≤n≤N T (b) particle velocity (v n ) 0≤n≤N T .

Figure 2 .

 2 Figure 2. One single particle motion. Numerical solution obtained with a large time step ∆t = 0.1 with the 1st order scheme (3.2), the 2nd order schemes (3.4)-(3.6) and (3.8)-(3.11) and the third order scheme (3.13)-(3.17) for ε = 10 -1 (top) and ε = 10 -2 (bottom) : (a) particle trajectory in physical space (x n ) 0≤n≤N T (b) particle velocity (v n ) 0≤n≤N T .

Figure 3 .

 3 Figure 3. Diocotron instability ε = 1. Time evolution of (a) relative energy E(t) -E(0) with respect to the initial one and (b) E(t) ∞ . For ε ∼ 1 the instability does not occur.

Figure 4 .

 4 Figure 4. Diocotron instability ε = 1. Time evolution of the density ρ for time t = 36, t = 38, t = 40 and t = 42 units. For ε ∼ 1 the instability does not occur.

Figure 5 .

 5 Figure 5. Diocotron instability ε = 10 -2 . Time evolution of the relative energy with respect to the initial one (a) for the Vlasov-Poisson system (1.3) (b) for the guiding center model (1.7).

Figure 6 .

 6 Figure 6. Diocotron instability ε = 10 -2 . Time evolution of the density ρ for time t = 60, t = 90, and t = 120 units for (a) the Vlasov-Poisson system (1.3) and (b) the guiding center model (1.7).

Acknowledgements

FF was supported by the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. views and opinions expressed herein do not necessarily reflect those of the European Commission.

LMR was supported in part by the ANR project BoND (ANR-13-BS01-0009-01).


Kunal Kumar
email: kkumar13@etudiant.mines-nantes.fr

Christian Javier Clavijo Lopez

Oscar Augusto Tellez Sanchez

Amit Gupta
email: agupta13@etudiant.mines-nantes.fr

Olivier Péton
email: olivier.peton@mines-nantes.fr

Thomas G Yeung
email: thomas.yeung@mines-nantes.fr

Adrien Vanuxem
email: vanuxem@sol3d.com

Integrated strategic and tactical optimization of animal-waste sourced biopower supply chains

Many models have been recently developed for the optimization of biomass related supply chains. However, models for biopower supply chains powered by animal waste have not received much attention yet. In this paper, we propose a mixed integer linear programming model for supplier selection and procurement planning for a biopower plant. The model integrates time window constraints for the collection of animal waste as well as inventory constraints. We show that the model is intractable with a state-of-the art commercial solver and propose a heuristic approach based on the Adaptive Large Neighbourhood Search (ALNS) framework. We show the efficiency of this approach on a case study in central France.

I. INTRODUCTION

The declining reserves of fossil fuels and an ever-increasing demand for energy has forced many economies to gradually shift their dependency for energy derived from non-renewable to renewable sources. The shift can be seen from the fact that between 2003 and 2012, the gross energy consumption of energy derived from renewable sources in France has increased from 12.1% to 15.6% [START_REF]Energy -Country factsheets version 3.0[END_REF].

Among the various sources of renewable energy produced in France, biomass holds a major share with 52% [START_REF]Market opportunities for biomass in France[END_REF]. Owing to a large agricultural area and immense resources of animal waste, the production of energy from biomass holds a lot of potential in France [START_REF] Nielsen | The future of biogas in Europe: Visions and targets until 2020[END_REF]. Moreover, the potential of biomass to produce energy in different forms -as biofuel (bioethanol) and, biopower (electricity), is really beneficial as well. At present, there exist more than 230 active biomass processing plants worth an installed capacity of 110MWe and these figures are expected to triple by the year 2020 [START_REF]Biogas plants in France market and new site potentials[END_REF]. This boost in the number of biomass processing plants is due to the increased feed-in tariffs and subsidies provided to the producers. However, despite the abundance and ease of availability of biomass (agricultural and animal waste), producing energy from biomass is costlier than other sources [START_REF] Kavkler | A Comparison of Electricity Generation Reference Costs for Different Technologies of Renewable Energy[END_REF]. In particular, transportation cost of animal waste is an issue because of its high water content [START_REF] De Figueiredo | A systemic approach for dimensioning and designing anaerobic bio-digestion/energy generation biomass supply networks[END_REF].

Thus, it is imperative to produce energy at the lowest possible costs for biomass to remain an attractive source of energy in the future.

Biomass collected from agricultural farms can be broadly divided into two categories waste coming from plants (agricultural biomass) and animal waste. Agricultural waste is generally produced during the harvesting season in high quantities. While on the other hand, animal waste is characterized by irregular production, and in small quantities. This imposes time-windows on the collection of animal waste, and thus differentiates it with the agricultural waste. In this paper, we consider the strategic selection of animal waste suppliers for a biopower plant and the tactical optimization of the collection schedule at suppliers. In must be noted here that the biomass considered in this study is restricted to animal waste, due to its characteristics which make it different from other biomass supply chains. In section II, we review the related literature to better position our work in the corresponding field of research. In section III, we propose a mixed integer linear programming (MILP) formulation. The MILP model minimizes the total costs, by taking into account the transportation costs incurred while procuring the biomass. We take into account the timewindows (minimum and maximum gap between consecutive deliveries) to align this problem closer to practice. The model also takes into account several practical constraints such as the collection of a minimum percentage of biomass from a selected supplier. This constraint does not exist in all situations, but is still common in many real-life situations. In section IV, we present a case study in central France and show that the CPLEX solver cannot solve it within acceptable computation time. In section V, we propose an Adaptive Large Neighborhood Search (ALNS) [START_REF] Pisinger | A general heuristic for vehicle routing problems[END_REF], [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF] meta-heuristic for the problem. We also run scenarios with different plant capacities and time discretizations. Section VI concludes this paper, and maps issues for further research.

II. LITERATURE SURVEY

Extensive literature reviews on the biomass supply chains have been proposed by [START_REF] Sharma | Biomass supply chain design and analysis: basis, overview, modeling, challenges, and future[END_REF], [START_REF] Mafakheri | Modeling of biomass-to-energy supply chain operations: applications, challenges and research directions[END_REF] and [START_REF] Gold | Supply chain and logistics issues of bio-energy production[END_REF] and can be referred for the current state of the art of the problem. In their literature review, Eskandarpour et al. [START_REF] Eskandarpour | Sustainable supply chain network design: an optimization-oriented review[END_REF] identify biomass to bioenergy supply chains as the leading application area of sustainable supply chain network design.

Studies have been made to reduce costs by the optimization of the biomass to bioethanol supply chains which take into account the various aspects of the production from harvesting, transportation, and handling from the supply to the end-customer (see [START_REF] Sharma | Scenario optimization modeling approach for design and management of biomass-to-biorefinery supply chain system[END_REF], [START_REF] Yue | Biomass-to-bioenergy and biofuel supply chain optimization: Overview, key issues and challenges[END_REF], [START_REF] Eksioglu | Analyzing the design and management of biomass-to-biorefinery supply chain[END_REF], [START_REF] Lin | Integrated strategic and tactical biomass-biofuel supply chain optimization[END_REF] or [START_REF] Zhang | An integrated optimization model for switchgrass-based bioethanol supply chain[END_REF]). However, these papers mainly focus on bioethanol production, which have different characteristics compared to biopower production systems. In bioethanol production systems, the end-product is usually supplied to several customers and the production is guided through the stochastic demand for the product. However, in case of biopower production systems, the output is injected into the power grids which serve the local populations. It is thus independent of the customers demand and boundedabove by the treatment capacity of the facility at each period. Further, there have been studies which take into consideration the different types of biomass like switch-grass [START_REF] Zhang | An integrated optimization model for switchgrass-based bioethanol supply chain[END_REF], cotton-stalk [START_REF] Tatsiopoulos | Economic aspects of the cotton-stalk biomass logistics and comparison of supply chain methods[END_REF] and miscanthus [START_REF] Shastri | Optimization of miscanthus harvesting and handling as an energy crop: Biofeed model application[END_REF] for the production of energy products. Most of these types of biomass are agricultural products and offer high seasonality due to limited periods of harvesting. Consequently, long-term storage is provided as a solution during seasons of no production [START_REF] Rentizelas | Logistics issues of biomass: The storage problem and the multi-biomass supply chain[END_REF]. As per the knowledge of the authors, there are no studies which consider animal waste as the input for biopower plants.

The procurement of biomass deals with two sub-problems (a) selection of sites from where the biomass waste can be collected and, (b) determination of the procurement plan from each selected site such that the biomass equivalent to the processing capacity of the plant is collected. While the first decision is made for a few years (strategic), the latter decision is taken for shorter periods, usually weeks or months (tactical). The selection of farms would affect the subsequent procurement schedule, which may be globally sub-optimal. Thus, it would be only beneficial to solve the two sub-problems as an integrated global problem, as is done by [START_REF] Lin | Integrated strategic and tactical biomass-biofuel supply chain optimization[END_REF].

The supply chain design for biomass processing industries has not received much attention in the research community until recently. However, recently this topic has interested a lot of researchers due to the growing awareness towards efficient energy production from renewable sources. Several optimization models have been proposed for decisions at strategic, tactical and operational levels. At the strategic level, the main problems are facility location and capacity determination. These problems have been studied by [START_REF] Lin | Integrated strategic and tactical biomass-biofuel supply chain optimization[END_REF], [START_REF] Kim | Design of biomass processing network for biofuel production using an MILP model[END_REF] and [START_REF] Balaman | A network design model for biomass to energy supply chains with anaerobic digestion systems[END_REF] among others. All of these studies propose MILP models to solve these problems. However, most of these problems are focused on biomass to bioethanol supply chains. At the tactical level, the key decisions include procurement schedules, transportation fleet management and inventory monitoring. These problems have been studied mainly for bioethanol production by several types of biomass. The studies have been made on miscanthus [START_REF] Shastri | Optimization of miscanthus harvesting and handling as an energy crop: Biofeed model application[END_REF], wood residues [START_REF] Svanberg | Analysing biomass torrefaction supply chain costs[END_REF] and cotton stock [START_REF] Tatsiopoulos | Economic aspects of the cotton-stalk biomass logistics and comparison of supply chain methods[END_REF]. The studies have been missing on supply chains with animal waste which are very different at operation and tactical levels.

There have been a few studies which have tried to integrate the two levels of decision as a single problem. Several mathematical models have been proposed to find the optimal facility location and biomass collection schedules [START_REF] Eksioglu | Analyzing the design and management of biomass-to-biorefinery supply chain[END_REF], [START_REF] Zhang | An integrated optimization model for switchgrass-based bioethanol supply chain[END_REF]. In addition to the above problem, a MILP model to address the problem of farm selection, resource allocation and operating schedule has been proposed in a study by Lin et al. [16]. This is one of the most comprehensive models for biomass to bioethanol supply chains. However, this study is focused on agricultural waste and does not take into the time-windows imposed by animal waste. Figueirido et. al. [START_REF] De Figueiredo | A systemic approach for dimensioning and designing anaerobic bio-digestion/energy generation biomass supply networks[END_REF] presents a systematic approach for the design of biopower supply chains with animal waste as biomass. However, the study does not present any optimization models for the design of the supply chains. Moreover, it proposes the selection of farms through a qualitative method and does not consider the time-windows.

As far as solution methods are concerned, most biomass supply chain design problems are solved with MILP and MINLP solvers (see for example [START_REF] Eksioglu | Analyzing the design and management of biomass-to-biorefinery supply chain[END_REF]- [START_REF] Zhang | An integrated optimization model for switchgrass-based bioethanol supply chain[END_REF]) which sometimes require extremely long computation times. With the addition of time-windows the problem becomes even harder, especially for a large number of farms and finer time discretization. Thus, we need a method which can provide good results in reasonable amount of time. As far as heuristic and metaheuristic methods are concerned, no method seems to be widely used by different authors. Adaptive Large Neighbourhood Search (ALNS) has proved to be really efficient for solving optimization problems such as vehicle routing problems, arc scheduling problems, pollution routing, snow-plowing and timetabling problem (see [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF], [START_REF] Laporte | An adaptive large neighbourhood search heuristic for the capacitated arc-routing problem with stochastic demands[END_REF], [START_REF] Masson | An adaptive large neighborhood search for the pickup and delivery problem with transfers[END_REF], [START_REF] Kristiansen | Adaptive large neighborhood search for student sectioning at Danish high schools[END_REF], [START_REF] Azia | An adaptive large neighborhood search for a vehicle routing problem with multiple routes[END_REF], [START_REF] Demir | An adaptive large neighborhood search heuristic for the pollution-routing problem[END_REF], [START_REF] Sorensen | International timetabling competition 2011: An adaptive large neighborhood search algorithm[END_REF] or [START_REF] Korsvik | A large neighbourhood search heuristic for ship routing and scheduling with split loads[END_REF]). However, for solving supply chain network design (SCND) problems, the use of Large Neighbourhood Search (LNS) has been explored by a few studies (see ref. [START_REF] Eskandarpour | A large neighborhood search heuristic for supply chain network[END_REF] or [START_REF] Copado-Méndez | Large neighbourhood search applied to the efficient solution of spatially explicit strategic supply chain management problems[END_REF]), but to the best of our knowledge, LNS (or ALNS) has not been used yet to solve integrated strategic and tactical problems, especially for biomass supply chain design problems.

The literature survey clearly identifies the need for studying biomass processing supply chains for animal waste which take into account its special characteristics. In addition to that, we need a solution method which provides good results in reasonable amount of time. Thus, our study tries to fill in these gaps by providing a model for the selection of farms and the determination of procurement schedules which take into account the inherent time-windows.

III. PROBLEM DESCRIPTION

In this problem, we try to optimize the procurement plan for a single biomass processing facility, with known monthly maximal demand, at integrated strategic and tactical level. It offers several characteristics like seasonal supply, timewindows on collection trips which make the problem interesting. We only consider the upstream leg of the supply chain as the output from the plant is injected into the grid. Downstream actors are not considered as the power is injected directly into the electricity grids. We also consider the facility as a black box, and any activity within the plant is out of the scope of this model.

We consider a single biomass processing facility for a given maximal capacity equivalent to the energy output of C tons of biomass per day. The production is a continuous process and exactly C tons of biomass (animal waste, in this case) must be provided each day. The biomass can be procured from a set F of candidate farms, which can be selected to be long-term suppliers of biomass. The decision of selection or non-selection remains unchanged for the whole planning horizon T . Each farm i ∈ F is located at a distance di from the facility. The availability of biomass at each farm depends on the farm size and shows high seasonality. For each period t ∈ T , biomass availability at farm i is given as P it . The length of each period t ∈ T is n days.

Animal waste offers specific characteristics compared to the agricultural waste. Although it is not limited to small periods (e.g. harvesting), it still offers a seasonal trend. Moreover, there are conditions imposed on the collection of animal waste due to practical reasons. The production of animal waste is spread both spatially and temporally; therefore, all the animal waste is usually assembled together after a certain period of time which could be of the order of a few weeks. This fact introduces a minimum gap between consecutive collection trips, as before this period the animal waste may not be accumulated. Two consecutive collections at a farm must be separated by at least g min periods.

Moreover, animal waste usually has a higher rate of biodegradation than the agricultural waste, due to its higher moisture content. Its semi-liquid state renders it difficult to store and handle. Thus, long-term storage is avoided [START_REF] Batidzirai | Optimisation potential of biomass supply chains with torrefaction technology[END_REF], [START_REF] Dornburg | Bioenergy revisited: key factors in global potentials of bioenergy[END_REF]. This imposes a restriction on the maximum gap of g max periods between two consecutive collections at a farm.

We assume that all the manure collected beyond the gap of g max periods is rendered invaluable due to excessive moisture content and therefore, low heating value. Additionally in some cases, mutual agreement dictates a collection of at least α% of the total biomass available at any farm, in case a farm is selected as a supplier. However, to maintain generality, this constraint can be easily relaxed by setting the value of α to zero in cases where this condition does not hold true. The quantity Wit of biomass left in the farm after each collection is disposed off and is not available anymore in subsequent periods. If there is no collection at farm i during period t -1, then the quantity J it of biomass available at period t is calculated as J it = J i,t-1 + P it .

In case the total available biomass from the selected farms is less than the required capacity of nC for the period, a quantity E t of biomass is secured from non-contractual suppliers (not selected) at an additional cost or penalty of A/ton. Furthermore, storage of collected biomass at the facility is of importance as well. Biomass is degradable and loses its energy value over time due to absorption of moisture. Thus, it is really important to process the biomass before its energy content reaches a threshold below which it is deemed as unacceptable. Any consumption until K days from the day of collection is considered acceptable. Any processing beyond that would result in a reduced production of energy, and hence a loss of LEU R/ton. This cost is in addition to the standard inventory costs for the biomass within the plant, charged at HEU R/ton/period.

The optimization problem consists of selecting a set of farms and planning the procurement at each farm in the time horizon while minimizing the total costs. The various parameters involved in the problem are defined in Table I. The variables of the mathematical model are introduced in Table II.

For the given notations, the problem can be modeled by the MILP (1)- [START_REF] Sharma | Scenario optimization modeling approach for design and management of biomass-to-biorefinery supply chain system[END_REF].

Q it d i + t ∈ T (HI t + AE t + LI t ) (1)
s.t.

I t = I t-1 + i∈F Q it + E t -nC ∀t ∈ T -{0} (2)
I t ≤ KC + I t ∀t ∈ T (3)
J i,(t-1)

+ P it = Q it + J it + W it ∀i ∈ F, t ∈ T -{0} (4) W it ≤ M i X it ∀i ∈ F, t ∈ T (5)
J it ≤ M i (1 -X it ) ∀i ∈ F, t ∈ T (6) α(J i,t-1 + P it ) ≤ M i (1 -X it ) + Q i t ∀i ∈ F, t ∈ T (7) t+gmin t =t X it ≤ Y i ∀i ∈ F, t ≤ |T | -g min ( 8
) t+gmax t =t X it ≥ Y i ∀i ∈ F, t ≤ |T | -g max (9) Q it ≤ M i X it ∀i ∈ F, t ∈ T (10)
J i0 = 0 ∀i ∈ F (11) Y it , X it ∈ {0, 1} ∀i ∈ F, t ∈ T ( 12
)
I t , I t , J t , W it ≥ 0 ∀i ∈ F, t ∈ T (13)
The objective function (1) minimizes the total cost which is composed of transportation costs, additional costs of buying biomass from non-contractual suppliers, cost of holding biomass surplus to the requirements of K days and standard inventory costs. Constraints (2) are the inventory balance equations for the plant. Constraints (3) calculate the inventory surplus to the requirement of K days. Constraints (4) balance the inventory at each farm. Constraints (5) permit discarding only in case of a collection. M i is a large number th at can be set at t t =0 P it . Constraints [START_REF] De Figueiredo | A systemic approach for dimensioning and designing anaerobic bio-digestion/energy generation biomass supply networks[END_REF] indicate that whatever is left after any pickup will be lost and so it will not be considered in the successive periods. Constraints [START_REF] Pisinger | A general heuristic for vehicle routing problems[END_REF] force a minimum amount of biomass to be collected, while constraints ( 8) and ( 9) impose the time-windows on the collections. Constraints (10) allow collection only when collection is permitted and only from the selected farms. Another point that must be noted is that the model contains both binary as well as continuous variables. However, the objective function depends only on the continuous variables.

IV. EXACT APPROACH A. Case Study

There is no existing benchmark for the problem. Therefore, we apply the MILP model to a case study corresponding to a biomass processing plant located in central France. The plant can be served by a pre-screened set of farms, out of which a set of optimal farms have to be determined. In this instance, 29 farms were considered for selection. The location of the farms is showed in Figure 1. For each of these farms, the production of biomass has been forecast based on historical records, which takes into account the inherent seasonality.

The capacity of the plant is considered fixed. However, in order to understand the structure of the problem, several capacity scenarios were studied. In addition to that, the problem was solved with different time discretizations. The various parameters for the problem are detailed in Table III.

B. Computational Results

Several scenarios were run on CPLEX Optimization Studio 12.6 on a Core i5-3230M CPU 2.60GHz, 6GB RAM. The runs were made for a runtime of 3 hours for several values of capacity C and period lengths n. We analyzed 3 different capacity values, namely, 60 tons/day, 70 tons/day and 80 tons/day. In this particular case, it must be noted that the total availability of manure corresponds to around 70 tons per day. Thus, the three chosen capacities represent the cases of under-consumption, total-consumption and over-consumption, respectively. The results are presented in Table IV.

From Table IV, it can be noted that the practical difficulty of the problem depends on both the time discretization as well as the capacity of the facility. It can be easily understood that the dependence on time discretization is due to the explosion of the number of variables for the same length of the planning horizon. However, the dependence on the capacity is not that straightforward. The instance where the capacity corresponds to the total available manure is easier to solve as it is in the interest of the problem to exclude all solutions which do not include selection of all farms. While the exact method works well for lower temporal resolutions and capacities greater than the availability of the manure, it becomes highly intractable for the other cases.

V. ADAPTIVE LARGE NEIGHBORHOOD SEARCH (ALNS) METAHEURISTIC

Since the problem becomes intractable for larger instances with finer time discretization, we explore the possibility of solving the problem through a metaheuristic. The Adaptive Large Neighbourhood Search (ALNS) was proposed by Ropke and Pisinger [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF] and is a variant of Large Neighbourhood Search (LNS) framework which was first proposed by Shaw [START_REF] Shaw | Using constraint programming and local search methods to solve vehicle routing problems[END_REF] to solve Vehicle Routing Problems. The metaheuristic is very similar to the ruin and create heuristic proposed by Schrimpf et al. [START_REF] Schrimpf | Record breaking optimization results using the ruin and recreate principle[END_REF]. An extensive review of the LNS and ALNS methodologies and their applications has been provided in ref. [START_REF] Pisinger | Large neighbourhod search[END_REF]. The basic mechanism of the LNS is to destroy and repair a solution iteratively using a combination of several destroy and repair heuristics. In ALNS, the destroy and repair heuristics are chosen dynamically based on their performance in the solution process. A general scheme of the ALNS framework is depicted in Algorithm 1.

Algorithm 1 Adaptive Large Neighborhood Search At each iteration, a part of the current solution is removed through some predefined Destroy mechanism to get a partial solution s (line 8). This solution s is then reconstructed using one of the predefined Repair heuristics to reach a new feasible solution (line 9).The solution s can be then accepted or rejected based on acceptance criteria (line 13), which may be defined by the user. The Destroy and Repair mechanisms are selected based on their past performance using scores ρ -and ρ + , respectively (line 6), which are updated at the end of each 100 iterations (line 16).

A. ALNS Implementation

The ALNS metaheuristic was coded in Java. It is modular and can be applied to the problem at hand by defining the following problem specific components: initial solution, acceptance criterion, stopping criterion. The various destruction and repair operators will be detailed in the next section. 2) Getting Initial Solution: In the initial solution, we consider that all farms are selected. In this situation, the problem can be considered equivalent to a multi-dimensional knapsack problem, where the capacity of the facility has to be filled with supplies from the selected farms. However, in this case, the supplies surplus to the capacity can be stored in the inventory, with inventory cost H and penalties for supplies of more than K days of requirement. We use the greedy mechanism proposed by Dantzig [START_REF] Dantzig | Discrete-variable extremum problems[END_REF] to solve the knapsack problem for each period sequentially. The farms which high quantity of waste per unit distance are preferred. If a farm is selected for one given period, it is not selected during the following g min -1 periods to maintain feasibility. Further, in case a farm is not selected for g m ax periods, then it is selected in the next period. This ensures time constraints satisfaction.

3) Acceptance Criterion: To avoid the solution from getting trapped in a local-minimum, a classical approach is to adopt an acceptance criterion based on the principles of simulated annealing [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]. If solution s deteriorates previous solution s, it can be accepted, with probability e -f (s )-f (s) θ , where s is the current solution and θ is a positive parameter called temperature. The initial temperature is set such that a solution that is 5% worse than the initial solution has an acceptance probability of 50%. This temperature is then reduced at the end of each iteration with a fixed cooling rate of 0.03%.

4) Stopping Criterion:

The algorithm is run for 10000 iterations.

B. ALNS Operators

This section details the operators used in our implementation of the ALNS metaheuristic.

1) Destruction Operators:

The destruction operators are used to remove some parts of the solution, so that newer solutions can be explored by re-constructing partial solutions. In this particular problem, we distinguish two ways of destroying a solution: a) altering the matrix Y of farm selection, b) altering the matrix X of farm collections.

Let us first describe the destruction operators of family a). Note that after altering matrix Y , matrices X and Q must be completely rebuilt. a1) Random De-Selection: This operator randomly deselects between 1 and 20% selected farms from matrix Y i . a2) Random Selection: In addition to preserving the selected farms from the incumbent solution, this operator randomly selects between 1 and 20% of the previously nonselected farms. This allows the construction of new X and Q matrices over a wider range of farms. a3) Random Toggle This operator randomly chooses between 1 and 0.2 × |F | elements in matrix Y and toggles the values between 0 and 1. The selected farms are set to 0, while non selected farms are set to 1. a4) Worst Output/Distance Ratio Swap: This operator sorts the lists of selected (S) and non-selected (NS) farms in the ascending order of total output through the planning horizon per distance ( t P it /d t ). Then, the r th element in the selected list is set as non-selected while the p th element in the nonselected list is set as selected. The value of r and p are calculated as shown in ( 14) and ( 15):

r = S [λ p |S|] (14)
p = N S [θ p |N S|] (15)
where λ and θ are random numbers between 0 and 1, p ≥ 1 is a deterministic parameter as described by Ropke and Pisinger [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF]. a5) Worst Surplus Removal: The last farm added during each period produces a surplus I t , which is undesirable. Thus a list L containing the last added farms during each period is maintained, and sorted non-increasingly based on their value of surplus. The q th element from the list de-selected is presented in [START_REF] Lin | Integrated strategic and tactical biomass-biofuel supply chain optimization[END_REF].

q = L [κ p |L|]] (16)
where, κ is a random value between 0 and 1, and p ≥ 1 is a deterministic parameter as defined in operator a4.

Let us now detail the destruction operators in family b). For the following operators, we assume that matrix Y is not destroyed. b1) Random Erase: For this operator, we randomly erase between 40% and 60% of the values in matrix X. We choose a high range of destroyed values in order to allow the heuristic to explore different solution regions. b2) Random Erase and Swap: This operator randomly erases x e % of the solution matrix. Out of the remaining (100 -x e )% of the matrix, we toggle 15% of the values. Only feasible toggles are performed, which do not violate the g min rule. The g max rule can be compromised in this phase of the algorithm and would be taken care of in the repair phase. While rule b1 allows preserving a part of the solution and building around it, there are still chances of patterns being repeated. The rule b2 forces the solution to have some characteristics opposite to the incumbent solution.

b3) Removal of collections with surplus: This operator forces the value of collections generating surplus to 0. The removed collections, along with the collections from the nondestroyed solutions, are kept intact for the repair phase.

2) Repair Operators: While the destruction operators allow the algorithm to navigate through a wider section of the solution space, the repair operators follow the path of intensification. They are more focused on reaching the best solution for a given partial solution. The following repair operators were employed during the study: c1) Greedy Repair This operator is almost the same as the construction heuristic used for the initial solution. However, in case of solutions for which the collections were partially destroyed, we consider the capacity of the period as a partiallyfilled knapsack. Then, we try to fill this capacity by the best feasible collection. The best collection is the collection with the highest value of biomass available per unit distance. c2) Noised Greedy Repair A noise of η i , very close to 1, is added to the value ν i of biomass per unit distance for each farm i ∈ F , such that the updated value for farm each farm is η i × ν i . Now, the knapsack is again filled with the feasible farms with the best noised values.

3) Selecting Destroy and Repair Operators: The destroy and repair operators are selected through an adaptive roulette wheel procedure. Each destroy (or repair) operator j is given a score ρ - j (or ρ + j ) based on its performance during the previous iterations. Higher value of ρ - j (or ρ + j ) suggests greater success of the operator for the given instance. The probability of selecting an operator from m different operators is given as The value of ρ j is updated at the end of each segment. We define a segment as 100 iterations. During each segment, we sum up the scores received by each operator for its contribution to a solution. An operator receives a score σ 1 if it provides a new global best solution, σ 2 if it provides a solution which has never been accepted before with objective value better than the current solution, and σ 3 if it provides a solution which has never been accepted before but with objective value worse than the current solution. Thus, at the end of each segment, the value of ρ j is updated as in equation [START_REF] Zhang | An integrated optimization model for switchgrass-based bioethanol supply chain[END_REF].

ρ j ← ρ j (1 -η) + η π j Ω j (17)
where, π j is the sum of scores received by operator j during the last segment, and Ω j is the number of times it was selected. The reaction factor η ∈ [0, 1] controls the rate of adjustment of ρ.

C. Assessment of the ALNS metaheuristic

Due to the stochastic nature of the heuristic, 10 runs for each scenario were made and the average computational results are provided in Table V. Columns 2 and 3 report the upper bounds obtained with Cplex (after 3 hours of computation) and the ALNS metaheuristic, respectively. Column 4 is the relative difference between the results on Columns 2 and 3. Column 5 indicates the computation time for the 10000 iterations of the ALNS.

It can be seen that the ALNS method shows significant improvement over the solution provided by CPLEX as we move towards the more constrained instances. The improvement is as high as 16% with capacity 60 and time discretization of 7 days. Thus, it can be seen that ALNS, in general, provides better and quicker solutions.

VI. CONCLUSION AND FUTURE WORK

The optimization of supply chain for biopower production from animal waste is a relatively unexplored field from an optimization research point of view. In this study we address the optimization of biopower supply chains with animal waste as the source of energy. We propose a mathematical model which takes into account the practical constraint of collection time-windows imposed by the inherent nature of the animal waste. We also consider other practical constraints which often tie facilities into collecting a minimum proportion of the total production. Then, we propose an ALNS metaheuristic that can solve the problem very efficiently.

The proposed model can be further improved by inclusion of transportation fleet management and a finer modeling of methane loss. Possible extensions include mixed sources of biomass (animal, energy crops, wood, industrial waste, etc.) or bi-objective optimization of both the logistic costs and the environmental impact.

Fig. 1 .

 1 Fig. 1. Location of the 29 pre-screened farms.

1 )

 1 Solution Encoding: A solution can be encoded as a set of three interdependent matrices Y , X and Q. Y is a binary matrix of dimension |F |×1 representing the selection of farms. X is a binary matrix of dimension |F | × |T | representing the collection periods at each farm. Q is a real positive matrix of dimension |F | × |T | representing quantities collected at each period at each farm.

TABLE I .

 I PARAMETERS USED IN THE MATHEMATICAL MODEL. Length of each time period (in weeks) C Processing capacity of the facility per day (in tons) K Number of periods before biomass start losing its power (in weeks) Pit The biomass availability at farm i ∈ F during period t ∈ T di Distance of farm i ∈ F from biopower plant (in km)

		Parameters
	F	Set of farms
	T	Set of time periods in the planning horizon
	n
	gmin	Maximum time gap between two consecutive collections for each farm
		(in weeks)
	gmax	Minimum time gap between two consecutive collections for each farm
		(in weeks)
	H	Inventory cost in the biopower plant (in /ton/week)
	Tc	Cost of transport per ton-distance biomass collected (/km-ton)
	A	Additional cost of buying biomass from non-contractual suppliers (/ton)
	L	Cost (loss of value + handling + storage) incurred to store biomass for
		more than K days (in /ton/week)

TABLE II .

 II PARAMETERS USED IN THE MATHEMATICAL MODEL. Continuous material flow variables (in tons) It Total inventory of biomass at the biopower plant during period t ∈ T Jit Total availability of biomass at farm i ∈ F during period t ∈ T I t Total amount of biomass that stays for more than K days at the facility during period t ∈ T Et Quantity of additional biomass bought during period t ∈ T Wit Amount of biomass discarded after each collection Binary and Integer Decision Variables

Yi 1, if a farm i ∈ F is selected; 0, otherwise Xit 1, if a collection is made from a farm i ∈ F during period t ∈ T Qit quantity of biomass collected from a farm i ∈ F during period t ∈ T (in tons) max T c t∈T i∈F

TABLE III .

 III PARAMETERS

	Parameter	Value
	|F |	29
	gmin	4 weeks
	gmax	13 weeks
	K	2 weeks
	α	85%
	n	7 or 15
	C	60, 70, 80

TABLE IV

 IV

	.	UPPER BOUND (UB), LOWER BOUND (LB) AND GAP
	AFTER 3 HOURS FOR SEVERAL VALUES OF C AND N.
		C	LB	UB	% Gap
			n=15 (time period of two weeks)
		80	51322.25	51568.60	0.48%
		70	33670.21	34010.28	1.01%
		60	19423.02	19720.19	1.53%
			n=7 (time period of 1 week)
		80	51958.30	53285.70	2.49%
		70	30840.80	35187.60	12.35%
		60	16964.1	20838.7	18.59%

1 :

 1 InitialSolution ← getInitialSolution() 2: BestSolution ← InitialSolution 3: CurrentSolution ← InitialSolution 4: ρ -(1, ..., 1) and ρ + = (1, ..., 1) 5: count ← 0 6: while termination criteria not met do Select Destroy and Repair Operator using ρ -& ρ +

	7:
	8:	s ← CurrentSolution
	9:	s ← Destroy(s)
	10:	s ← Repair(s)
	11:	if s ≤ BestSolution then
	12:	BestSolution ← s
	13:	CurrentSolution ← s
	14:	else if acceptance(s, CurrentSolution) then
	15:	CurrentSolution ← s
	16:	end if
	17:	count ← count + 1
	18:	if count mod 100 = 0 then
	19:	update ρ -and ρ +
	20:	end if
	21: end while
	22: return BestSolution

TABLE V .

 V COMPARISON OF UPPER BOUNDS ACHIEVED BY CPLEX AND THE ALNS ALGORITHM.

	C	CPLEX	LNS	% Gap	CPU LNS (sec)
		n=15 (time period of two weeks)
	80	51568.60	51523.24	0.09%	141.4
	70	34010.28	33896.40	0.33%	139.8
	60	19720.19	19444.14	1.40%	140.2
		n=7 (time period of 1 week)
	80	53285.70	52489.30	1.49%	142.6
	70	35187.60	30897.36	12.19%	143.6
	60	20838.70	17411.8	16.44%	141.7