A Ghizzo

P Bertrand

T W Johnston

N Besse

Lecture note ISSS9 The 9 th International School of Space Simulations: Eulerian Vlasov models for plasma simulation: review and advances

We present some problems for which Vlasov simulations could bring real break-thought to study short pulse high intensity laser-plasma interaction for plasmas relevant to conditions met in the National Ignition Facility or the Laser MégaJoule under construction in France. Such a Vlasov code allows a fine description of parametric -like instabilities in underdense and overdense plasmas and their saturations, particle trapping, particle acceleration and nonlinear resonant wave-particle interaction till plasma conditions relevant to the fast Ignitor concept.

Introduction

Halfway between the N body problem and the well-known hydrodynamical (fluid) model, the Vlasov equation (supplemented by the Poisson or Maxwell equations) describes self-consistent particle dynamics in a variety of physical systems, including space and laboratory plasmas, stellar dynamics or thermonuclear fusion. The numerical integration of the Vlasov equation is one of the key challenges of computational plasma physics. Since the early days of this discipline, an intensive work on this subject has produced many different numerical schemes, which however, can be bunched together in two main groups. On the one hand, Particle-In-Cell (PIC) codes have proven to be useful in studying plasma dynamics even in 2D or 3D problems and complex geometries and on the other hand the famous Vlasov models.

As far as collective phenomena are our concern, it is worth remembering that in a Vlasov plasma, Coulomb interaction between charged particles are actually replaced by a mean field calculated from Poisson or Maxwell equations using charge density (or current density) where the microscopic fluctuations (due to the fact that a plasma is not a continuum) are averaged over the Debye length. This mean field concept is thus the basic idea which points to the concept of PIC model. In PIC codes, particle trajectories are computed from a field prescribed on a fixed grid with a mesh size of order the Debye length, which is indeed equivalent to compute the characteristics curves as given by the Vlasov equation. At the end of the time step, the charge of each particle is redistributed among the neighboring mesh points, allowing to solve Maxwell's equations. This method yields satisfying results with a relatively small number of particles. However only a few particles per cell have been used in PIC codes leading to a high level of numerical noise (due to the invidual effects), especially in regions of phase space where the density is low. These PIC codes, for situations relevant to laser-plasma interaction or gyrokinetic model may largely overestimate the plasma heating because of their intrinsic numerical heating.

We present here the basis of the numerical schemes used for solving Vlasov_Maxwell system. Some problems for which Vlasov simulations could bring real break-thought to study high intensity laser plasma interaction. It is beyond the scope of this paper to present an exhaustive catalogue of all the results obtained in this research field. Moreover, we would try to help the reader to select the problems in which Vlasov simulations could be pertinent both from the numerical and physical points of view.

The structure of the paper is organized as follows. Part 2 is devoted to phase space properties and the difference between the conservative and advective forms of the Vlasov equation. In part 3 we quickly survey th numerical schemes for the interation of the Vlasov equation in an eulerian approach and in particular the notion of time spliiting.. Sect. 4 discusses how to choose between PIC (particle-in-cell model) and Vlasov codes for a given phnysical problem. Part 5 we present the important role of the minority particle populations with particular emphasis on trapped particles in nonlinearLandau damping and finally an application of wave-partcicle interaction in part 6 concerning laser-plasma interaction at low laser intensities. The complex problem of eulerian Vlasov solvers in the relativistic regime is studied in part 7 and we conclude in section 8.

2 Phase space properties for a hamiltonian system: the Vlasov model.

We consider the motion of N particles in phase space (x, p) interaction through selfconsistent Coulomb forces and external fields. To simplify the notation and without any loss of generality we supposed an one-dimensional space. Thde traditional approach to plasma kinteic theory, in the case of electrostatic systems, starts from the Vlasov-Poisson system coupling the particle distribution function f (x, p, t) with the electric field E (x, t). Thus f (x, p, t) dxdp is the number of particles at time t having positions in the range between x and x + dx and momenta in the range p and p + dp.

The temporal evolution of f leads to a description of the system more detailed than a fluid description but indeed less detailed than following the trajectory of each individual particle. (i) The first phase space property concerns the conservation of the phase space volume dS 0 = C 0 dxdp of a close contour C 0 defined at time t 0 . At a later time t the particles on the initial contour have moved and define a new contour C, the volume of which is C dxdp = dS. If the motion is Hamiltonian, i.e. if ẋ = ∂H ∂p and ṗ = -∂H ∂x , then we have:

dS = dS 0 (1)
The second fundamental property is that these particles which at time t 0 are inside the contour C 0 are exactly the same, which at time t, are located inside C. Now letting dS 0 and dS tend to zero, we can define a phase space density f (x 0 , p 0 , t 0 ) and f (x, p, t). From f (x, p, t) dS = f (x 0 , p 0 , t 0 )dS 0 and dS 0 = dS, we deduce then f (x, p, t) = f (x 0 , p 0 , t 0 ) or

Df Dt = 0 (2)
The kinetic limit reads as the conservation law for the distribution function f , leading to the Vlasov equation:

∂f ∂t + dx dt ∂f ∂x + dp dt ∂f ∂p = ∂f ∂t + v ∂f ∂x + eE(x, t) ∂f ∂p = 0 (3)
Since the plasma is a Hamiltonian system, it is then possible to write the Vlasov equation in aconservative form:

∂f ∂t + ∂ ∂x (vf ) + ∂ ∂p (eE(x, t)f ) = ∂f ∂t + div X (U f ) = 0 (4)
where X = (x, p) stands for the phase space coordinates and U = (v, eE) is a divergence-free advection field. Now when and under what conditions is the Vlasov model valid? This is indeed a complex question. In a Vlasov equation the particles' interaction is entirely described by a field computed via the Poisson equation:

∇.E = e ε 0 f dp - en 0 ε 0 (5)
or more generally by the Maxwell's equations in the electromagnetic case. Discrete particle (charge e, mass m e and of density n 0 ) imply the existence of collisions inside the plasma, which may become strong when the impact parameter is shorter than e 2 /ε 0 k B T and where k B T is the thermal energy (k B being the usual Boltzmann's constant). This gives a mean free path of

= n -1 0 k B T ε 0 e 2 2 ( 6
)
and a collision frequency ν e v th (7) with v th = k B T /m e which is the thermal velocity. Using the well-known plasma frequency expression ω p = n 0 e 2 /m e ε 0 , we see that a condition of validity will be given by:

ν e ω p 1 (8)
or equivalently

ν e ω p v th = 1 n 0 λ 3 D (9)
where λ B = k B T ε 0 /n 0 e 2 is the (electron) Debye length. Consequently from (8) and (9) the Vlasov equation (used in the collisionless limit) is valid for the graininess parameter g = 1/n 0 λ 3 D which tends to zero (i.e. for n 0 λ 3 D 1). In particular in a one-dimensional case, rather than collisions it is more convenient to speak about fluctuations connected to the discrete nature of phase space fluid while both concepts -collisions and field fluctuations-agree in real 3D space (note that in the exact onedimensional N body problem, in which particles are described by charged plans, the concept of collision has no meaning). An usual way to characterize the plasma is to define a dimensionless parameter g in such way, that if we perform the dichotomy experiment (first introduced by Rostoker and Rosenbluth (1960) in Ref. [1], which consists in cutting in two parts each particle, this parameter g will be also divided by two. A systematic analysis was developped by Lotte and Feix (1984) [2] for different systems. They look for the length and time sacle invariant under this dichotomy. It can then be shown that λ D and ω -1 p are the length and time scale invariant under this dichotomy and that the Vlasov equation remains invariant and we get g = 1 n 0 λ d D (10) where d is the dimentionality of the plasma (d = 1, 2 or 3). As a matter of fact the purely collective approach-as Vlasov model-are approximation of the N -body problem; A simple calculation delineates this problem. A plasma of length L allows (L/λ D ) 3 collective modes while the importance of the individual effects is given by n 0 λ 3 D . Fusion plasmas exhibit a graininess factor g of order 10 -6 to 10 -8 and space plasmas have still smaller values of g, while L/λ D ∼ 10 3 -10 4 . Obviously we cannot treat 10 18 particles (since N = n 0 λ 3 D (L/λ D ) 3 . This number corresponds indeed to L/λ D = 10 3 and n 0 λ 3 D = 10 9 ). It was consequently realized very quickly that a numerical solution of the Vlasov equation -which is based on the limit n 0 λ 3 D → +∞ -was a more appropriate approach.

At this stage, two remarks must be pointed out:

• First, particle codes, which correspond also to the Vlasov limit, reintroduced an individuality effect. particles are considered as volume elements of phase space gathered at the same points (x i , p i ). An efficient smoothing is then introduced either by considering extended in space overlapping particles or by filtering the short wavelengths of the particle among their neighboring points of the configuration mesh (particle-in-cell). Although these treatments decrease efficiently the unwanted -because exagerated-individual effects they do not change the scaling of these effects and dividing these effects by two still imposes a doubling of the total number of particles. The reason is that, to kill these effects, we must smooth the fields not only the interparticle distance but also on the Debye length, but in the last case we begin to modify the possible collective behavior. Thus, even in numerical treatment n 0 λ d D of order 10 3 must be used.

• Second, The Vlasov approach, if it deals with the correct treatment of the purely collective plasma, does it by introducing a phase space i.e. a product of the configuration space by the momentum (or velocity) space. A sampling of this phase space for a plasma of dimension d (and length L) requires

(L/λ D ) d N d p ,
where N p is the number of samples needed for each momentum component for a precise description of the distribution f (x, p, t). Typically N p = 10 2 or 10 3 is needed. Thus a 3D case (i.e. 6D in phase space) remains a difficult problem and a natural challenge.

Eulerian Vlasov codes

An explanation of the fact that little effort was spent on Vlasov codes during the '80s was connected to the limited size of the available computer and up to very recently these codes could only treat 2D or 3D phase space. But another explanation must be also discussed. Previous attempts on the numerical solution of the Vlasov equation were not able to solve, in a practical way, the complex problem of filamentation in velocity space. The source of the problem lies in the treatment of the free streaming term of the Vlasov equation:

∂f ∂t + v ∂f ∂x = 0
which has a general solution in the Fourier k-space in the form:

f (k, v, t) = f (k, v, 0) exp(ikvt) (11)
Thus from Eq. ( 11), it is clear that f (k, v, t) oscillates with respect to v at frequency kt. When the time t is such that this frequency reaches the inverse of the velocity mesh size v, i.e. t max (k v) -1 , we cannot anymore followthe exact solution of f . This is the unavoidable filamentation phenomena with its related propagation to large Fourier modes in velocity. Thus in the treatment of f bu eulerian codes we may have to interpolate and this operation may lead to numerical instability. The code must avoid these unstabilities and provide a smoothing of the details of the velocity space, hoping that the details will not influence the subsequent evolution of the plasma. Knorr (1963) [3] gave the first solution through a spectral method (a double Fourier transform both in space and velocity x, v) but considered only weakly nonlinear problems close to the linear theory. The choice of the Fourier transform was rather obvious for the configuration space: it has the advantage that in a linear theory all k components are independant and therefore slightly nonlinear problems could be expected to be described by a small number a Fourier modes. The same choice for the velocity space was less obvious and another treatment was introduced by Feix and Grant (1967) [4] and Armstrong (1967) [5] Fourier transforms were kept naturally in x but Hermite functions were introduced in the v space:

f (x, v, t) = 1 √ 2π +∞ s=0 a ns e ink 0 x H s (v) e -v 2 /2 (12)
where we have

H s (v) = (-1) s exp v 2 2 d dv s exp - v 2 2 ( 13
)
H m (v) H n (v) exp - v 2 2 dv = √ 2πn!δ nm (14)
The Hermite series has advantages (it does not introduce spurious periodicity in v in particular it is possible to recover landau damping with a reasonable number of Hermite polynomials) but has the great disadvantage of a slow convergence due to the fact that, for large n, H n (v) cos ( √ nv) or sin ( √ nv). Consequently large number of Hermite polynomilas -up to 10 4 in some calculations was needed, nevertheless the treatment was efficient enough to take into account small nonlinearities, initially located on long wavelengths (kλ D < 1). More recently Schumer-Holloway (1998) [6] have used rescaled orthogonal basis-like the so-called scaled Hermite basis to provide long time numerical stability while using a small set of basis functions. Note also, that an interesting way to remedy the filamentation problem was proposed by Klimas (1987) [7], who introduced a filtering term into the Vlasov equation to eliminate the high frequencies (introduced by filamentation effects) at the price of a smoother approximation. But in discussion given in Figua et al (2000) [8] it was shown that the efficiency of the method is indeed very limited.

The solution for the correct treatment of the Vlasov-Poisson equations was finally given by Cheng and Knorr (1976) in ref. [9]. The method is the basis for modern Vlasov codes. It introduces two ideas:

• The time splitting, into two steps. In the first we treat the free streaming part:

∂f ∂t + v ∂f ∂x = 0 (15)
and then the acceleration part:

∂f ∂t + eE(x, t) m e ∂f ∂v = 0 (16)
• The reconstruction of the distribution function on a mesh in phase space. Both Eqs. (15 and (16) have analytical solutions and in section 3.1 we will give the detailed algorithm.

Although the solution was given in 1976 -the method was not systelatically used before 1982 with the study of different problems in the electrostatic regime (nonlinear Landau damping [START_REF] Feix | Eulerian codes for the Vlasov equation[END_REF], BGK holes stability or two stream instability [10], Wigner equation describing quantum [START_REF] Manfredi | Transport equations for quantum plasmas[END_REF] and also in electromagnetic regimes (especially laser plasma interaction models). Since that time, (semi-lagrangian) Vlasov codes are slowly taking their place among the plasma simulation tools.

Resolution via the Splitting Scheme

We outline here some physical meaning of the splitting scheme method in the case of the Vlasov-Poisson system. We want to solve the Vlasov equation for the onedimensional electron plasma distribution function f (x, v, t) using usual normalized quantities

∂f ∂t + v ∂f ∂x + E(x, t) ∂f ∂v = 0 (17)
where the electric field is given self-consistently by the Poisson equation:

∂E ∂x = +∞ -∞ f (x, v, t)dv -1 (18)
starting from a given initial condition f (x, v, t = 0) and assuming, for instance, the existence of periodic boundary conditions. The Vlasov equation ( 17) is then integrated in the original phase space by applying a splitting scheme (indeed a method of fractional steps), which consists of treating the convective term and the acceleration term separately, as previously mentionned. The crucial point is the representation of the electric field E (x, t) by a succession of Dirac pulses. Thus we have to solve the Vlasov equation (17) in which the electric field E(x, t) has been replaced by the term E * (x, t)given by:

E * (x, t) = E (x, t) t n δ t -t n+ 1 2 ( 19
)
where t n = n t. Denoting t -

n+ 1 2 and t + n+ 1 2
as the time before and after the Dirac pulse, the integration during the time interval t from t n = n t to t n+1 = (n + 1) t is straighforward and can be divided into four steps:

1. We solve the free particle motion for t n ≤ t ≤ t -

n+ 1 2
and obtain the solution in the usual form:

f * (x, v) = f n x - v t 2 , v (20)
2. Compute the electric field at time t -

n+ 1 2
by substituting f * in the Poisson equation (18).

3. Compute, for the acceleration term, the quantity:

f * * (x, v) = f * (x, v -E * t) (21)
4. Repeat again the step 1 to obtain finally:

f n+1 (x, v) = f * * x - v t 2 , v (22)
It must be pointed out that applying successively this sequence, we obtain now

f n+1 (x, v) = f n x -t v - E * t 2 , v -E * t where E * = E x -v t 2
, which is equivalent to the following integrated equations of the characteristics of the Vlasov equation:

x n = x n+1 -t v n+1 -E x n+ 1 2 , v n+ 1 2 t 2 (23)
v n = v n+1 -E x n+ 1 2 , t n+ 1 2 t ( 24
)
The continuum from of the characteristic equations, is of course, given by

ẋ = v and v = E (x, t) (25)
Thus, from a numerical point of view, the splitting scheme is equivalent to an integration of the Vlasov equation along their characteristics and is correct to the second order in time step. From a physical point of view, the fact that Eqs. ( 20) and (21) give an exact solution of the Vlasov equation (17) with the replacement of E by the field E * provide some interesting information on the size of t, as already pointed out by Ghizzo et al (1988) (see [10]) . For instance, linearizing the Vlasov equation with E * around a homogeneous equilibrium (i.e. f (x, v, t)

= n 0 F 0 (v) + f 1 (x, v, t))
and performing the usual Laplace (t → s) in time and Fourier in space (x → k) transforms yields:

-ik

   E (k, s) + i k n E k, s - 2iπn t +∞ -∞ dF 0 dv dv s -ikv    = +∞ -∞ f 1 (k, s, t = 0)
s -ikv dv (26) Due to the appearance of the sum n on the left-hand side, the plasma dispersion relation cannot be revovered. But, since we knom that the frequency spectrum of the electric field does not exceed a few ω p (choosing the time step t such that we have tω p 1, results in a zero contribution of all terms of the discrete sum, except for n = 0, which allos us to recover the usual landau dispersion relation.

Interpolation methods

Now, we have to consider the problem of performing the shifts in x (from (20) and in v (see Eq. (21). Different methods can be used. It can be pointed out that during each time step, the Vlasov equation takes on the form of an advective equation (one over x and one over v). Cubic spline interpolation in the velocity space and Fourier interpolation in the configuration space (without using Fast Fourier transforms) have been used by Cheng and Knorr (1976) (see ref. [9]), which requires, in the x direction, an execution time proportional to N 2

x (where N x is the number of points in the spatial direction). Because this fact increases the computational effort when a large number of points is used, a cubic spline interpolation can then be used also in this direction. On the other hand, the fact taht a shift can be expressed in the corresponding Fourier space (by just changing the corresponding phase) suggests the use of fast Fourier Transforms (FFT).

Fourier in x /Fourier in v interpolation method

Let us recall briefly the main features of the Fourier in x / Fourier in v interpolation scheme. By performing a Fourier transform in the xspace (k being the conjugated position variable), Eq. ( 20) leads to an exact solution in the form:

f * (k, v) = f n (k, v) exp -ikv t 2 (27)
while a Fourier transform in the v space (λ being now the conjugate of the velocity space variable) applied to Eq. ( 21) gives:

f * * (x, λ) = f * (x, λ) exp (-iλE * t) (28)
Thus this sequence consists to apply a FFT, alternatively in each direction (x and v ), to multiply by the corresponding phase, to finish to make an inverse FFT. This last code has been used by [START_REF] Ghizzo | [END_REF] [START_REF] Ghizzo | [END_REF]10] to study the nonlinear behavior of BGK waves: the code has been found to be extremely stable over time of order 10 3 ω -1 p and this scheme is well adapted to treat the quantum Liouville equation (see Suh et al (1991 in [14]), which gives the time evolution of the Wigner function. Such a representation plugs a quantum mechanical problem in the familiar phase space environment. Many good review articles on the Wigner formalism already exist in the literature (see [START_REF] Feix | An universal model, The Vlasov equation[END_REF] [START_REF] Feix | An universal model, The Vlasov equation[END_REF]), to which the interested reader may refer for further details.

First the so-called Wigner equation is given by

∂f w ∂t + p m e ∂f w ∂x = ie 2π 2 V x + λ 2 , t -V x - λ 2 , t (29)
f w x, p , t exp i (p -p ) λ dλdp

In this equation, f w usually represents the Wigner function, which is defined as follows:

f w (x, p, t) = 1 2π dλρ x - λ 2 , x + λ 2 e -ipλ/ (30)
where ρ (x, y) is the density matrix. We recall that, for a pure quantum state, the density matrix can be written ρ (x, y) = ψ * (x) ψ (y); ψ being the usual wave function satisfying a Schröedinger equation. The splitting scheme can be usefully applied also to the numerical solution of the Wigner equation. We again split the equation into two parts and first solve the kinetic part.

Step 1 and 4 are straighforward. Taking the Fourier transform of step 1 yields to the following solution:

f * w (k, p) = f n w (k, p) exp -ik p m e t 2 (31)
Taking the inverse transform we get f * w (x, p). To perform step 3 (the potential step) we first notice that the integral (with respect to p ) in the Wigner equation (29) has the form of a convolution product. Therefore we perform a Fourier Transform in pspace (λbeing now the Fourier conjugate of p) and get the solution:

f * * (x, λ) = f * (x, λ) exp - ie t V x + λ 2 , t -V x - λ 2 , t (32)
Notice, at this point, that the Fourier transform on p turns the treatment of the potential term into a problem nearly as simple as the classical one that means A Vlasov model. It turns out that the nonlocal character of the quantum interaction appears very clearly in replacing in Eq. (32) a derivative (∂V /∂x = -E) by its centered finite difference counterpart.

Disadvantage of the method

While cubic spline interpolations seems to be well adapted for simulations with periodic and causal boundary conditions, the Fourier method can lead, in the treatment of the velocity space, to the occuring of aliasing in the Fourier transform (the Fourier transforms imply that the distribution function is periodic in velocity, but regardless of its initial gaussian shape, it will be generally evolue to a function that is non periodic. Consequently to avoid aliasing a large nmber of values of the distribution function above a given cut-off in velocity is set equal to zero. But this usual method for controlling this edge noise by filling the grid with zero's near its velocity boundaries, leads to a non-negligible increase of the memory, togetherwith the computational CPU time.

Flux conservative method

The starting point of this method is the Flux Balance method (FBM) (see for instance Fijalkow (1999) in ref. [16]), discretizing the Vlasov equation in the conservative form. This is a technique based on application of the ENO scheme method (for more details see ref. [17]). The basic idea of this method is to compute the average of the Vlasov equation solution in each cell of the phase space grid by a flux conservative law. We first observe that, by a standard time Splitting scheme, we can restrict ourselves, without loss of generality,to a one-dimensional case:

∂f ∂t + ∂ ∂x (u (x, t) f ) = 0 (33)
Again we can define the characteristic curves solution as

dX dt (s) = u (X (s) , s) (34)
and X (t) = x. Let us denote by X (x, s, t) the solution of (33) and define the Jacobian J (x, s, t) = ∂ x X (x, s, t) and the solution of (33) can be written in the standard form:

f (x, t) = f (X (x, s, t) , s) J (x, s, t) = f (X (x, s, t) , s) ∂X ∂x (x, s, t) (35)
which describes the conservation of the particles along the characteristics. Let us introduce a finite set of mesh points x i+ 1 2

and assuming the values of f at time t n+1 by integrating f on each subinterval (cell). The flux conservation becomes then

f n+1 i = f n i + φ i-1 2 (t n ) -φ i+ 1 2 (t n ) x ( 36
)
were

f n i = 1 x x i+1/2 x i-1/2 f (x, t n ) dx (37)
and

φ i+ 1 2 (t n ) = x i+1/2 X(x i+1/2 ,tn,t n+1) f (x, t n ) dx (38)
The main step is then to choose an efficient method to reconstruct f from the values in each cell. Fijalkow (1999) in [16] only used a linear interpolation but this method does not give a positive scheme and does not control spurious oscillations. It must be pointed out, that the initial scheme proposed by Boris and Book (1973) [18] or the use od classical slope limiters is too dissipative to give an accurate description of the distribution function. An improved version is the positive and flux conservative (PFC) method (Filbet, Sonnendrücker and Bertrand, 2001 in [19]), which is not only conservative of course, but also preserves the positivity and the maximum value of f . Different test problems like free particle advection, plasma wave echoes and nonlinear Landau damping, have been studied. Among the different schemes presented in Filbet et al (2003) in [20], there is no clear winner, each scheme having its pros and cons. An interesting study for the case of Vlasov-Maxwell system can be found in Mangeney, Califano et al (2002) in [21]. Furthermore new reconstruction methods, namely weighted essentially nonoscillatory (WENO) have been recently introduced (Shu 1998 in ref. [START_REF] Shu | In advanced numerical approximation of nonlinear hyperbolic equations[END_REF]) . They can be applied to any hyperbolic system of conservation laws, and clearly the Vlasov equation belongs to this framework. They have the double advantage of being able to deal with shocks or steep gradients while achieving high order accuracy in the smooth regions. This makes them the most precise methods available for conservation laws. Preliminary results for the Vlasov-Poisson system are very encouraging (see for instance Labrunie S., Carillo, Bertrand , 2004 in [23]).

P.I.C. versus Vlasov codes

The questionarises now of the choice between PIC and Vlasov codes, even though PIC simulations (like any Monte Carlo) become numerically more interesting as the dimension increases. First, it must be pointed out that solving Poisson's equation (or Maxwell's equations) requires the same computational effort for both PIC and Vlasov models, and needs the same spatial grid, with a mesh size x of the order of the Debye length λ D . On the other hand solving the Vlasov equation is nothing else but solving the characteristic equations for one particle (since f is the one-particle distribution function). This is completely dequivalent, in principle, to solving the motion equations for one particle in a PIC code. Therefore it is lear that pushing one particle in a PIC solver needs the same numerical effort as reconstructing a phase space mesh point in an Eulerian (Vlasov) solver.

Thus the ratio between the numerical cost for a PIC code and a Vlasov code (CPU time as well as memory requirement) will scale as the ratio

N vlas N part
where N vlas is the total number of mesh points in phase space (for the vlasov equation) and N part is the total number of "super-particle" in the corresponding PIC methods. Defining N V as the number of grid points required for the sampling of the distribution function in the velocity space, N vlas can be written simply as the product of the spatial grid by the velocity grid

N vlas = L x dx N d V V ( 39
)
where d V (d x ) is respectively the dimension of the velocity space (d V = 1, 2 or 3) (and of the configuration space with d x = 1, 2 or 3). As we have previously pointed out, a plasma of length L allows (L/λ D ) dx collective modes and the importance of the individual effects is usually described by the factor g or its inverse n 0 λ dx D . Therefore, for a PIC code, N part can be written in the following form

N part = n 0 L dx = n 0 ( x) dx L x dx ( 40
)
and we see the occuring of the factor n 0 ( x) dx which is nothing else but the inverse of the the graininess parameter due to particle discreteness of the PIC model

g P IC = 1 n 0 λ dx D ( 41
)
This relation deserves a few comments. A PIC code exhibits a paradoxal situation where the introduction of a spatial grid x ∼ λ D allows one to describe the collective Table 1: PIC versus Vlasov code: ratio

N vlas /N part for N V = 100 d V = 1 d V = 2 d V = 3 g P IC = 10 -2 1 10 2 10 4 g P IC = 10 -4 10 -2 1 10 2 g P IC = 10 -6
10 -4 10 -2 1 effects through the mean field approximation but at the price of the reintroduction of individual effects due now to the finite number of superparticles. These individual effects are often considered as a numerical noise whose measure is given by the factor g P IC as given by Eq. ( 41). Thus N part may be now written as

N part = g -1 P IC L x dx (42)
Finally from Eqs. (39 and (42) we get the ratio

N vlas N P IC = g P IC N d V V (43)
To illustrate this relation, consider a typical value N V = 100, which is quite reasonable to have a fine sampling of the velocity space (we will come back to this point in section 7). From (43) we can draw the table I. Although these numbers must be considered with caution, they ive the general ideas for deciding the use of PIC or Vlasov Eulerian code. Below the principal diagonal Vlasov codes offer a better resolution at a smaller price, while above the diagonal, PIC codes must be preferred.

• The crucial point is the dimension of velocity space and the value of the factor g required for describing the physical problem.

(i) For 1D velocity space, which is indeed the case for gyrokinetic modeling since in that case the only kinetic variable is the longitudinal velocity v (parallel to the magnetic foeld), Vlasov codes must be preferred.

(ii) For 2D velocity space the question is still open and depends on the level of noise which is bearable in the PIC code: if a small level is needed to study the onset of instabilities in laser-plasma interaction (Ghizzo, Huot and Bertrand, 2003 in ref. [24]) (with g P IC smaller than 10 -4 , for instance) a Vlasov code has to be used. On the contrary, for gross phenomena (where a small g P IC is not needed) a PIC code will do the job at a lower price.

(iii) For 3D velocity space unless there is a need of ultra low noise with very small g P IC , PIC codes remain the only possibility for a long time.

• The use of a low g P IC value is needed when kinetic effects are dominant, and in particular for resonant wave-particle interaction, as considered in section 6. This points to the use of Vlasov models. Moreover, in the PIC models, particles are usually loaded so as to reproduce the velocity distribution function f . This means that the velocity space resolution on a tail region of f , where usually Landau resonance or wave-particle interactions take place, would be coarser than in other regions where f peaks. On the contrary, in a Vlasov code, phase space resolution (provided the sampling is uniform) is always guaranteed, whatever the velocity profile of f .

• the appearance of Vlasov codes have allowed a precise observation of the structure of the phase space especially in these regions where f was very small. It is interesting to notice that while the importance of phase space holes (which may characterize Bernstein-Green-Kruskal-BGK-like waves [25]) was well reconized in phenomena like the two-stream instability, their importance in other phenomena had to wait the undertaking of the Vlasov codes (see for instance [START_REF] Ghizzo | [END_REF] in ref. [START_REF] Ghizzo | [END_REF]).

Role of minority populations: from nonlinear Landau damping till KEEN waves

Landau damping is the oldest and most famous kinetic effect exhibited by the Vlasov equation. It may be worthwhile to mention, that Vlasov in the linearized treatment of this equation missed an aspect connected to analytical continuation in the complex plane. Landau, in 1946 (see ref. [26]), gave the correct but obscure treatment of the problem. Consequently a long discussion took place until 1962, some authors cleaming that it was a mathematical artefact with no physical aspect until Dawson (1965) [27] and independly Eldridge and Feix (1962 in ref. [28]) showed its existence through computer experiments. A more direct check was finally given by the eulerian codes of Knorr (1963) (Fourier-Fourier) and the Hermite-Fourier algorithms of Feix, Grant et al in [29] (1967) and then spline interpolation by Shoucri and Gagné (1978) in ref. [30]. These authors notice also that an increase of the nonlinearity (i.e. the size of the perturbation) brings a stopping of the damping after some oscillations with possible modulation of the amplitude of the oscillations, and finally, the establishment of a steady undamped oscillation. A clear physical picture of this phenomena has waited the appearance of graphical logicial allowing a representation of the 2D phase space f (x, v) and in 1989 Bertrand, Ghizzo et al in [START_REF] Bertrand | Nonlinear phenomena in Vlasov Plasma Doveil Eds[END_REF] showed phase space pictures obviously, could not have been obtained particles codes.

Such a Landau damping problem in the nonlinear regime, is one of the basic test used today to check the validity and the accuracy of new algorithms for Vlasov solvers.

The problem is the following: the plasma is unidimensional, electrostatic with only electrons moving (ions forming a homogeneous neutralizing background). The problem is given by the following initial condition

f (x, v, t = 0) = F 0 (v) (1 + αcosk 0 x) (44)
where

F 0 (v) = 1 √
2πv th e -v 2 /2v 2 th is the Maxwellian distribution function. Both theoretical and numerical studies indicate that after an initial Landau damping, the amplitude of the wave changes in time in an oscillatory manner, due to a periodic exchange of energy between the electrons trapped in the potential well of the plasma wave, and the wave itself. As time goes on (tω p 300) the amplitude oscillations tend to constant level. A rough explanation was given through the use of quasilinear theory arguments. The help of the high resolution phase space diagnostics of the Vlasov code allows to give a much more precise explanation. A first series of experiments was performed with α = 0.1, k 0 λ D = 0.3 in a periodic box of length L = 2π/k 0 . (The Landau damping is then γ L /ω p 0.018 with a phase velocity of v ϕ 3.9v th ).

The upper part (v 0) of the phase space is represented in Fig. 1, using grey shading for distribution function values less than 10 -3 .

1

We see learly on the first pictures of Fig. 1 the overtaking and the subsequent breaking of the plasma wave in the phase space, followed by the formation of a hole (a BGK-type structure) with a filamentary structure winding around this hole. The rotation of these filaments inside the hole structures have indeed the bounce frequency and are responsible of the oscillations of the amplitude of the electric energy. After a long time, the filaments have disappeared due to the phase miwing, and the persistence of the oscillation is simply related to the "balistic" motion of the holes ( each hole of positive v hole velocity is associated with another hole of negative velocity). This hole's motion at constant velocity associated to the spatial periodicity with wavenumber k 0 , ives the corresponding frequency ω = k 0 v hole , obviously close to the Bohm-Gross frequency since v hole is close to the phase velocity of the wave v ϕ . Note that the formation of the holes corresponds to the presence of a small population of trapped particles but which plays the fundamental role in the stopping of the Landau damping. Reducing α with the same k 0 λ D gives rise to smaller holes and consequently a smaller asymptotic energy level. Note that the presence of the holes was obtained with a modest numerical effort (the treatment of N x N V = 64 × 256, i.e. 16384 mesh points). With this number of particles, particle codes hardly suggest the holes.

Many more interesting results have been obtained with the Vlasov codes especially concerning the fusion of holes. Figs 2 and 3 describe the coalescence of two holes obtained initially by nonlinear Landau damping leading to a smaller structure. This simulation corresponds to the following initial conditions:

f (x, v, t = 0) = 1 √ 2πv th exp - 1 2 v v th -βsin2k 0 x -αsink 0 x 2 (45)
with k 0 = 025, β = 2 and α = 0.05. This fusion leads to a wave number divided by two, since we have usually an acceleration term of the new hole, the fequency of the driven oscillation kv k is not quite divided by 2 but jumps from ω = 1.15ω p to ω = 0.82ω p ,i.e. well below the plasma frequency -a characteristic of low-frequency Kinetic Electrostratic Electron Nonlinear wave (KEEN) (for more details see ref. [32]). The concept of these waves are a novel version of a non-stationary Bernstein-Greene-Kruskal (BGK) nonlinear waves (1957), (see ref. [25]) with electrons trapped in the wave troughs and characterized by self-sustained electron holes in phase space.

This brings us to an important question of the physical origin of these lowfrequency kinetic waves.

Kinetic wave concept is however becoming of wider and wider importance, and one of these reasons is the remarkable set of results obtained by Montgomery et al (2001) , in a series of laser-plasma interaction experiments [33]. These experiments were aimed at improving our knowledge of Stimulated Raman Scattering (SRS) from Electron Plasma waves (EPWs) from a single speckle of laser light and employed laser scattering as a key diagnostic. In addition to the extected EPW scattering signal, there was also an unexpected signal which was much smaller than that associated with SRS-related EPW (∼ 10 -3 lower intensity) but nonetheless well out the thermal noise. The extremely surprising aspect of this scattering result was that the frequency of the relatively weak electrostatic density wave associated with this signal was only a modest fraction (∼ 0.37) of the plasma frequency, which put it right in the so-called frequency gap where no electrostatic wave ought to be, at least according to linear wave theory for a plasma with Maxwellian electrons. Note that the observed wavevector was at a modest fraction (0.26) of the reciprocal Debye length, giving the electrostatic wave a phase velocity not far above the electron thermal velocity.

To explore the physics of such generation and to make contact with a possible laboratory experiment using the ponderomotive force (PF) between two opposing laser beams appropriately separated in frequency, it was decided (beginning in 2001) and published by T.W. Johnston et al (2009) in [32] to undertake a study of the effects of applying a PF drive to a collisionless uniform Maxwellian plasma in an appropriate computer model for an 1D electrostatic Vlasov electron plasma. By varying the carrier frequency of the drive, it has been found that not only can KEEN waves be sustained nonlinearly and self-consistently, but that a wide range of frequencies, hitherto thought to be inaccessible for coherent collective wave excitation (in a band gap between linear Ion Acoustic Wave (IAW) and Langmuir wave) were all legitimately sustainable as well. The simulations demonstrated successful generation of indefinitly self-sustained low-frequency waves, produced by using a PF drive which (as it turned out) needed to be applied for a limited time only. In effect the successful drive results in a trapped electron distribution well distributed in the trapping well so as not to remove energy from the wave. The probe results for the case when the drive is sufficiently strong to produce a sustained KEEN wave are shown in Fig. 4, obtained with a drive of a 0 = 0.05 (where a 0 is the normalized ponderomotive drive amplitude). The available simulation results suggest that a criterion of necessary drive strength and duration for producing a surviving KEEN wave can be established in terms of the characteric electron bounce (or trapping) time τ B relevant to the formation of an electron trapping vortex in phase space. It seems that the electron plasma response for untrapped electrons proves to be well approximated by linear theory over a sufficiently long (up to τ B /2) time during the drive. The resulting normalized electron bounce (or trapping) time

τ B ω p = 2π [(1 + χ e ) /kλ D a 0 ] 1/2
where χ e (ω/ω p , kλ D ) is the electron susceptibility. Thus the factor 1+χ e is a large correction at the low frequencies for the case at hand. The plasma response is essentialy that of an imperfect but fairly effective shield. For our usual parameters ω/ω p = 0.37, kλ D = 0.26 we have χ e = -1.21 + 9.59i and the normalized trapping time is then τ B ω p = 2π [9.6/(O.26a 0 )] 1/2 = 38.2/ √ a (this is the formula used to indicate τ B in Fig. 4).

The evolution of the phase space snapshots is shown in Fig. 5, at the same times t/τ B = 1 2 , 1, 2, 10 with the total-force potentials superposed. For the evolution of the distribution function the t/τ B seems to work quite well, so one can describe the evolution of a KEEN wave fairly generally as a function of t/τ B . The development of the distribution function up to t/τ B = 1 2 reflects a very normal early evolution of the distribution function. Between t/τ B = 1 2 and t/τ B = 1 one sees the development of a well-organized vortex in phase space which would seem like to survive, just as the drive is being terminated. With the drive already turned off for some time (about τ B ) at t/τ B = 2, the phase space vortex is showing some complicated evolution with fine ripples and structures. By t/τ B = 10 the distribution function in Fig. 5 appear to have evolved to relatively smooth states which resemble BGK-like states, at least superficially. Thez question of whether there is a continuing evolution to an ideal BGK state of whether fluctuations about this ideal state do not die out with time must be left to investigations with more sophiticated diagnostics than those employed here.

Stimulated Raman scattering (SRS), wave action and manley-Rowe partition

The acceleration of charged particles by fast large amplitude longitudinal electron plasma wave trains has received considerable attention in recent years. One of the potentially most effective ways to produce such a wave uses the forward Stimulated Raman Scattering (SRS) process, long of interest for laser fusion. For these problems, an eulerian relativistic Vlasov code renders possible a detailed examination of the low density regions of phase space (see ref. [34,35]). For these problems particle simulations codes suffer from poor statistics. This is because the PIC codes lack enough simulations particles to display the detailed phase space structure of the distribution function which is often obtained in those regions of phase space where particle and phase velocities are comparable and where the trapping occurs.

Here we present the Vlasov simulation results on action evolution for a periodic plasma and comparisons with the canonical three-oscillator model. In order to handle the plasma wave velocities relevant to future particle accelerators, the model equation must at least be relativistic in the acceleration direction (x). For the linear (y direction) laser polarization intensities we consider here, the transverse dynamics can be economically included in the nonrelativistic transverse cold-fluid approximation. The Vlasov equation in the acceleration (x) direction is thus:

∂f ∂t + p x m e γ ∂f ∂x + e (E x + u y (x, t)B z ) ∂f ∂p x = 0 (46)
where γ = 1 + p 2 x /m 2 e c 2 is the Lorentz factor. In the y direction, for reasons of computer economy, the cold nonrelativistic approximation is used:

∂u y ∂t = eE y m e (47)
To assist in the understanding of kinetic effects met in laser-plasma interaction, Manley-Rowe partition is very useful. Resonant wave-particle interaction plays an important role in plasma physics and is responsible for a number of phenomena a simple classic fluid model cannot indeed describe. Among laser-plasma phenomena, SRS is a very good candidate to check the pertinence of Vlasov kinetic simulation rather than fluid modeling and further the use of a Vlasov code rather a PIC code. Remember first that SRS is a three-wave parametric instability which involves the decay of a pump wave (ω 0 , k 0 ) into another electromagnetic wave (ω s , k s ) and an EPW (ω e , k e ) satisfying the phase matching relations:

ω 0 (k 0 ) = ω s (k s ) + ω e (k e ) and k 0 = k s + k e (48)
We focus here on the case of a forward propagation of the (Stokes) scattered electromagnetic wave (i.e. with a wave number k s > 0) the so-called forward stimulated Raman scattering SRS-F. Before introducing Vlasov simulations it is worth recalling some predictions obtained from a simpler fluid model and in particular concerning the Manley-Rowe partition. We obtained for the longitudinal electric field component E x and the potential vector A ⊥ = A ⊥ (x, t)e y (i.e. considering a linearly polarized electromagnetic wave), taking the Lorentz factor close to one:

∂ 2 E x ∂t 2 + ω 2 p E x -3v 2 th ∂ 2 E x ∂x 2 = - 1 2 eω 2 p m e ∂A 2 ⊥ ∂x (49) ∂ 2 A ⊥ ∂t 2 + ω 2 p A ⊥ -c 2 ∂ 2 A ⊥ ∂x 2 = eA ⊥ m e ∂E x ∂x (50)
Now making use of the usual assumption of a slowly varying envelope for each wave the last part of the calculation involves a multiple Krylov-Bogolioubov-Mitropolski expansion of Eqs. ( 49) and (50) in the conventional form:

A ⊥ (x, t) = 1 2 εA 0 (x 1 , t 1 )e i(k 0 x 0 -ω 0 t 0 ) + 1 2 εA s (x 1 , t 1 )e i(ksx 0 -ωst 0 ) + c.c. (51)
E x (x, t) = 1 2 εE e (x 1 , t 1 )e i(kex 0 -ωet 0 ) + c.c. ( 52
)
together with the following expansion where ε → 0:

∂ t = ∂ t 0 + ε∂ t 1 and ∂ x = ∂ x 0 + ε∂ x 1 (53)
To the first order in ε we get the three-wave envelope model. Defining the complex action amplitude a 0,s,e such that the action density S = aa * is given from the energy density W by S = W/ω, we obtain then:

∂ ∂t + v g0 ∂ ∂x a 0 = -Λa s a e (54)
∂ ∂t + v gs ∂ ∂x a s = Λa 0 a * e (55)
∂ ∂t + v ge ∂ ∂x + η a e = Λa 0 a * s ( 56
)
where a i = A i 0 ω i /2 for i = 0, s and a e = ω -1 p E e 0 ω e /2. Here the coupling coefficient is Λ = (e/2m e ) (2 0 ω 0 ω s ω e ) -1/2 k e ω p and η a phenomenological damping term. Integrating Eqs. (54) to (56) over a box of length L yields to the usual Manley-Rowe relations:

D Dt (N 0 + N s ) = F 0 + F s (57) D Dt (N 0 + N e ) + 2ηN e = F 0 (58)
where we have defined the quantity +∞ -∞ a i a * i dx = N i with i = 0, s, e. N 0 and N s represent the mean number of pump and scattered photons in the slab per unit area while N e is the mean number of plasmons per unit area. The quantities F i = v gi (S i (x = 0, t) -S i (x = L, t)) / denotes the photon flux of type i entering the slab at x = 0 minus that exiting at x = L. Moreover, the periodicity condition in space leads to the well-known Manley-Rowe equations written in term of action densities S i rather than quasi-particles. Eqs. (57) and (58) become when η = 0: We focus now on the numerical results obtained through Vlasov semi-lagrangian simulations in a periodic case. The time behavior of the action densities from a Vlasov simulation for SRS-F is shown in Fig. 6 Kinetic effects can be observed by introducing a small hot electron population in the initial distribution function which allows to implement particle trapping. The plasma was indeed chosen with two electron temperature components, the majority (95% in electron density) component with a 15keV temperature (high enough for electron Landau damping to subdue the usually rapidly growing but here unwanted SRS-B) and a minority (5%) at higher temperature (of 100keV ) to enhance wave-particle interaction around the phase velocity in phase space. The action is transferred back and forth between the pump and the daughter waves in a conventional fashion with an accumulating loss due to the plasma wave and the action transfer to trapped particles. As expected by Eq. ( 59), the first exact invariant C s = S 0 + S s is well conserved, while the pump-plus-plasma wave sum C e strongly decreases (see Fig. 7). This is clearly the signature of kinetic effects due to the wave-particle interaction involving a plasma wave that can trap electrons and thus lose action. By taking into account this "kinetic" action loss S kin due to the occurring of the electron trapping we might expect to observe a good conservation of the second (exact) Manley-Rowe invariant. This kinetic contribution in density action may be determined by integrating the kinetic energy of the trapped particles which are located above the lower separatrix in phase space. We have shown in ref. [35] that this integral is equal to the action density S kin may be estimated by:

C s = S 0 + S s = C s (t = 0) = const
S kin = W kin ω e = m e c 2 ω e L L 0 dx +∞ p low (γ -1) f (x, p x , t) dp x (61)
where W kin is the relativistic energy density of particles with a momentum greater than the lower separatrix trajectory directly computed from the Vlasov code. It is clear that the total action sum C e = S 0 + S e + S kin gives now a very good action conservation, as can be seen in in Fig. 8). Fig. 9 shows the x -p x phase space representation of f afforded by our Vlasov code, at six different times during the evolution when the hot electron population is introduced. Our idea was to follow test-particle trajectories in the potential created by the plasma wave. These calculations can be very simple if we introduce two hypothesis:

(i) We have neglected the ponderomotive force in the Vlasov equation (46), which is valid in the case at small laser flux i.e. when a 2 osc 1.

(ii) The plasma wave is stationary in a frame moving at the phase velocity of the plasma wave v ϕ . From the expression of the Hamiltonian in the laboratory frame:

H lab = eΦ lab (x -v ϕ t) + m e c 2 (γ -1) (62)
where Φ lab data are directly obtained by our Vlasov code, the new hamiltonian in the frame of the plasma wave may be determined from the Lorentz transform. A little algebria leads to the following form:

Figure 7: For SRS-F in periodic plasma and in a weak kinetic regime, the introduction of a small hot electron population (with 5% in density at 100keV ) modifies the particle dynamics. We do not obtain here the conservation of the plasma plus pump action due to kinetic effects. Physical parameters are identical to those used in Fig. 6.

H w = γ ϕ eΦ lab (x -v ϕ t) + m e c 2 (γ -1) -v ϕ p x (63)
Here γ ϕ is the Lorentz factor corresponding to the phase velocity. Using the property that this new hamiltonian is now a constant of the motion, it is now possible to determine the separatrices. At the X-points we have indeed γ = γ ϕ and the hamiltonian takes the form:

H sep w = γ ϕ eΦ 0,lab + m e c 2 γ -1 ϕ -1 (64)
where Φ 0,lab is the maximum potential amplitude. If H w < H sep w the particle is trapped; with the opposite inequality it is of course detrapped. From the separatrix orbit condition, one can readily solve the resulting quadratic condition for normalized separatrix momentum value:

p ± (x) m e c = u ϕ 1 + eγ ϕ m e c 2 ∆Φ ± 1 + eγ ϕ m e c 2 ∆Φ 2 -1 (65)
where u ϕ = γ ϕ β ϕ and β ϕ = v ϕ /c and ∆Φ = Φ 0,lab -Φ(x).

Figure 8: For a periodic SRS-F case, ime evolution of wave action, as in Fig. 6, together with C e and S part where S part accounts for particles accelerated above the lower separatrix (see Fig. 9) and C e is now the total action sum (pump plus plasma plus particles).

In Fig. 9 we present the separatrices calculated in this way superimposed on the phase space representation of f . Indeed the plasma wave traps a significant fraction of particles. This resulting trapped particle population behaves as a "macro-particle" moving at the phase velocity v ϕ = ω e /k e = 0.884c. A striking agreement can be seen between the separatrices ans the phase space trapping limits for growing wave: the spiral structure inside the trapping limits reflects the monotonic increase in wave energy.

At this stage, two remarks can be pointed out:

• First, particles interacting resonantly with large-amplitude coherent plasma wave can be trapped and subsequently detrapped. The long-time dynamics of such accelerated particles (located around the phase velocity in the resonance region) may become chaotic throughout a large region of the phase space when repeated wave-particle interactions occur. This chaotic behavior is located near the separatrices in phase space. Hamiltonian stochasticity may be obtained even for a simple plasma wave but with slowly time-periodically oscillating amplitude. Several Vlasov simulations, carried out by A. Ghizzo et al (2009) in [36], show that particle dynamics may become diffusive and that the corresponding spreading of particle trajectories (leading to the chaotic diffusion) is here due to the separatrix crossing, i.e. successive trapping and detrapping by the slowly oscillating plasma wave. If the trapped particles located near the separatrices have just adiabatic motion in the well, they would have to "cross" the separatrix when the amplitude of the field decreases. It must be pointed out that this separatrix crossing is not only located near the X-points but along the upper separatrix in phase space.

• The second point concerns the so-called kinetic regime of Stimulated Raman Back-scattering. Previous Vlasov simulation have shown the importance of trapping effects on plasma evolution. Such kinetic effects corespond to regimes which are relevant to plasmas conditions expected in ignition designs to be elded at the National Ignition Facility (NIF), as well as the laser mégaJoule project (LMJ) in France.

We focus now on this last point. In the weak Landau regime (kλ D ≤ 0.3), it is well known that the main mechanism responsible of the saturation of SRS-B is the Langmuir Decay Instability (LDI), but the evidence that LDI limits the amplitude of SRS-B was only gathered in that regime. In the regime where kλ D ≥ 0.3 (i.e., at high plasma temperature or low density), referred to here as the kinetic regime of the instability, the situation proves to be more complex. SRS-B reflectivity presents a bursting behavior, building up and breaking down in cycles. Here trapped electrons play a major role, so the nonlinear behavior of the plasma is dominated by kinetic effects and not by fluid effects, through an initial fluid phase and its transition to the kinetic one was recently identified by Kline et al. and presented in ref. [37] . At low kλ D values, LDI was observed via the existence of multiple EPWs while, as kλ D increases, the amplitudes of the daughter waves from the LDI cascade drop below detection threshold a single frequency-broadened EPW spectrum was observed.

It is clear that in that kinetic regime, trapping effects and the associated nonlinear shift in frequency (due to the nonlinear mode coupling or particle trapping ) are expected to dominate since the LDI process is too heavily Landau damped to compete. Kinetic effects may however take different aspects. In the Vlasov simulations discussed in Albrecht-Marc and Ghizzo (2007) in [38] , the main process that limits the Langmuir wave three-wave resonant frequency-wavevector chirping of SRS-B in its "kinetic" form (i.e. when LDI does not dominate) is provided by the process which begins with pair-wise merging of phase space holes and was here labeled "inverse cascade" (in which there is a continuous flow toward lower wavenumbers and lower frequencies, although, perhaps, "downward chirp" might be more appropriate). This is the process which is seen to break up the chirped resonance tracking in the kinetic regime, after which dramatic electrostatic wavenumber lowering and broadening is apparent.

The structures resulting from SRS-B instability growth are characterized by the formation of strong trapping vortices have in effect modified the plasma enough so that the linear Landau damping has been destroyed. After these vortices are well formed the waves are found to be very stable and and self-sustained. (In other simulations they indeed persist over a long time after the pump was turned off.) These waves were first generated as classic Langmuir waves plasma waves associated with SRS-B, with enough energy invested to overcome the considerable linear Landau damping and with normalized frequency and wavenumber here nonlinearly lowered by Morales and O'Neil retuning (see ref. [39]) to become Langmuir-Morales-O'Neil (LMO) waves with normalized frequency and wave vector of ω B = ω B -δω 0.25ω 0 and k B = k B + δk 1.65ω 0 /c. (It may well be the case that enough drive to obtain well-organized phase-space vortices in this kinetic regime (k B λ De ≥ 0.3) will automatically require significant LMO re-tuning, i.e., the two will go hand-in-hand.)

To aid in understanding what is happening, we present in Fig. 10 snapshots of the electron distribution function (p x versus x ) for the time tω 0 = 5600. The x-space has been divided in 32 cells with each all having a length of 125 cω -1 0 . This simulation of SRS-B corresponds to an inhomogeneous profile with injected counter-propagating pump (ω 0 , k 0 ) and seed (ω s , k s ) light with intensities of I seed = 10 -4 I pump . The physical system consists of a plasma parabolic profile of 3900cω -1 0 length, surrounded by a vacuum of 50cω -1 0 on both sides. The total length of the system is then L = 4000cω -1 0 corresponding to 609.7λ 0 (or 214 µm with an incident pump wavelength of λ 0 = 0.351µm). The quiver momentum of the pump wave (assumed to be here linearly polarized) is a osc = 0.025 which corresponds to an intensity of the pump laser beam of I pump 7 × 10 15 W cm -2 . Here we choose T e = 2keV for the electron temperature. Ions were kept fixed to simplify the presentation but mobile ions do not change the physical processes considered in this paper since LDI does not play a fundamental role for high values of the parameter k e λ D . The top panel of Fig. 9 clearly exhibits the trapping and formation of holes. As a consequence of the wave number shift ke , there is a change in the nature of the initially EPW, which is no longer a simple Langmuir wave at later times, but as we hypothesize, a nonlinear wave characterized by trapping effects in the form of (non self-sustained) phase space holes. The wave vector was found to be close to k e c/ω 0 1.58 in the top panel of Fig. 10. These new waves were first generated as classic Langmuir waves associated with SRS-B and propagating away from the laser, but with normalized frequency and wavevector here nonlinearly lowered by Morales and O'Neil retuning. The second cell shown in the middle panel in Fig. 10 is located now in the right region after the top of the parabolic profile. We see the emergence of a more turbulent state. The vortices begin to coalescence by pairwise vortex merging and lose their coherence. The bottom panel shows a third cell now located more deeper in the right region of the plasma: one see 23 trapping structures (resulting from such a previously pairwise vortex merging) corresponding to an average value of the wave vector of k V M c/ω 0 1.15.

The main process that limits the Langmuir wave three-wave resonant frequencywavevector chirping of SRS-B in its kinetic regime is provided by the process which begins with pairwise merging of holes. This is this process which is seen to break up the chirped resonance tracking in the kinetic regime (for more details see ref. [47]). The vortex merging of EPWs of shifted frequencies leads to the formation of lower wavenumber BGK-type structures (KEEN waves) whose the self-sustained character allows to persist over a long time. It must be point out that such a self-sustaining character of the resulting holes from the initial merging was numerically observed using Vlasov simulations after the pump wave was stopped. Because the structures resulting from the SRS-B instability growth were first generated as classic Langmuir waves, with enough energy invested to overcome the considerable linear Landau damping and with normalized frequency here nonlinearly lowered by Morales and O'Neil retuning, we have called these new waves Langmuir-Morales-O'Neil (LMO) waves. Rather than continuing this retuning three-wave resonance progress indefinitely, these coherent LMO waves are unstable and began to merge, a process that broke up the three-wave resonance. The resulting waves (at k V M c/ω 0 ∼ 1) are found to be very stable and self-sustained. 7 The relativistic 1D Vlasov equation

Using canonical Invariants to reduce the dimension of the velocity space

To describe the behavior of an high -intensity electromagnetic wave propagating in a relativistic electron gas in a fixed neutralizing ion background, we need now to solve the relativistic Vlasov equation. The details of the model were extensively described by Huot, Ghizzo and Bertrand (2003) in [40], for an one-dimensional plasma, but nonetheless we find it useful to recapitulate some of the main steps in our model and in particular the possibility to use canonical invariants to reduce the dimension of the momentum space. Even in the 1D spatial case, the dimension of the momentum space is high (here d v = 3). But when the characteristics of the Vlasov equation possess an exact invariant, a drastic reduction in phase space can be done. Even for the plasma wave propagating let us say along the x-direction, we have to solve a Vlasov equation for a 4D distribution function F (x, p, t):

∂F ∂t + p x m e Γ ∂F ∂x + e E+ p × B m e Γ ∂F ∂p = 0 (66)
where the Lorentz factor has the conventional form:

Γ 2 = 1 + p 2 m 2 e c 2 (67)
Let us consider the Hamiltonian of a particle in the electromagnetic field (E, B), in the relativistic regime:

H = m e c 2 Γ + eφ (x, t) (68)
with

Γ = 1 + (P c -eA) 2 m 2 e c 2 1 2 (69)
and where φ is the electrostatic potential, A the vector potential and P c the canonical momentum connected to the particle momentum p by

P c = p + eA (70)
Choosing the Coulomb gauge (divA = 0) we have then A = A ⊥ (x, t) i.e. which is purely transverse. The reduction occurs here since one has a class of initial conditions that are invariant under the dynamics. The Hamilton equation writes, along the longitudinal direction (which is also the direction of propagation of the electromagnetic wave):

DP cx Dt = - ∂H ∂x (71)
and in the transverse direction

DP c⊥ Dt = -∇ ⊥ H = 0 (72)
This last equation (72) leads to the property of the conservation of the transverse canonical momentum P c⊥ = const = C. If one assumes that the initial distribution function is located on (Dirac) delta-function spikes, at some locations, then the spikes are maintained under the hamiltonian dynamics. The full distribution function can thus be written as a sum of Dirac δ-distribution:

F (x, p x , p ⊥ , t) = Q j=1 f j (x, p x , t) δ [p ⊥ -(C j -eA ⊥ (x, t))] (73)
Therefore we can, without loss of generality, consider a plasma initially prepared so that particles are divided into Q bunches of particles, each bunch j having the same initial perpendicular canonical momentum P c⊥ = C j . For one bunch, it is possible to take this constant to zero. In other words, it means that all electrons at given (x, t) have the same p ⊥ = -eA ⊥ of the canonical momentum and the 4D distribution function can be reduced to a 2D distribution function f (x, p x , t). Upon solving the reduced Vlasov equation, one obtains an exact solution of the Vlasov model, albeit one with a singular initial condition, corresponding to cold plasma distribution in the perpendicular direction. Because of their exactness, both types of reduction inherit a hamiltonian structure from their parent models. Using the corresponding hamiltonian form, which writes here as:

H = m e c 2 1 + p 2 x m 2 e c 2 + e 2 A 2 ⊥ (x, t) m 2 e c 2 + eφ (x, t) (74)
we obtain for the reduced Vlasov equation, using the Liouville's theorem:

Df Dt = ∂f ∂t + [f, H] = 0 (75)
where the brackets are those of the Poisson bracket. Finally Eq. ( 75) reads

∂f ∂t + p x m e γ ∂f ∂x + eE x - e  m e γ ∂A 2 ⊥ ∂x ∂f ∂p x = 0 (76)
where the new expression of the square of the Lorentz factor is then given by

γ 2 = 1 + p 2 x m 2 e c 2 + e 2 A 2 ⊥ (x, t) m 2 e c 2 (77)
Eq. ( 76) is completed by the (longitudinal or "plasma") electric field equation:

∂E x ∂x = e 0 (n(x, t) -n 0 ) with n(x, t) = +∞ -∞ f (x, p x , t)dp x (78)
and by the vector potential equation

∂ 2 A ⊥ ∂t 2 -c 2 ∂ 2 A ⊥ ∂x 2 = J ⊥ 0 (79)
with

J ⊥ (x, t) = - e 2 A ⊥ m e +∞ -∞ f γ dp x (80)
Both from an analytical and numerical points of view the use of exact invariants allows obviously a reduction of the complexity of having to deal with the full Vlasov equation in a 4D phase space. Extension of the method to several particle bunches was made by Bertrand et al (2005) in [41].

The problem of the time spitting and semi-lagrangian scheme

Now, hoping the reader is convinced of the virtues of Vlasov codes, we would like to sketch the main features of the state of the art. Let us start with the well-known Cheng and Knorr scheme described above in section 3.

In the simple case considered above, the advection term v (respectively eE(x, t) in the electrostatic case) does not depend on the variable along which the shifts are performed x (respectively v), and the computation of

f * (x, v) = f n (x -v t/2, v) or f * * (x, v) = f * (x, v -eE * t/m)
are straightforward. This method proved to work very well in this electrostatic case. The advantage is that the scheme can use larger time steps than explicit Eulerian ones and is noiseless in comparison to particles-in-cell codes, but the price to pay is to reconstruct a regular grid using interpolations. However, this scheme could not be applied as easily to more complicated Vlasov problem, for instance to gyro-kinetic Vlasov models or the relativistic Vlasov equation (76), in which the Lorents factor (77) couples the different components of momentum p.

Let us consider the Vlasov equation written in the advective form:

∂f ∂t + U (X, t) .∇ X f = 0 (81)
where X stands for the phase space coordinates and U is a divergence-free advection field. Note that since U is divergence free, Eq. ( 81) can also be written in an advective (transport) form:

∂f ∂t + div X (U(X,t)f ) = 0 (82)
Splitting the components of X into two sets X 1 and X 2 , Eq. ( 82) can be rewritten in the following form:

∂f ∂t + div X 1 (U 1 (X 1 , X 2 ,t)f ) + div X 2 (U 2 (X 1 , X 2 ,t)f ) = 0 (83)
Moreover, it is well-known (see for instance ) that solving separately

∂f ∂t + div X 1 (U 1 (X 1 , X 2 ,t)f ) = 0 (84) ∂f ∂t + div X 2 (U 2 (X 1 , X 2 ,t)f ) = 0 (85)
keeps the second order accuracy for the whole Eqs. (81) and equivalently (83) by alternating the solvers. This is the basis of the time-splitting scheme introduced by Cheng and Knorr. Indeed it is important to point out that the semi-lagrangian scheme does not solve Vlasov's equation in its conservative form, but in its advective form to make full use of the backward characteristic method, provided that we have:

div X 1 (U 1 (X 1 , X 2 ,t)) = 0 and div X 2 (U 2 (X 1 , X 2 ,t)) = 0 (86)
Note that in PIC model, the Vlasov solution is obtained using the trajectory of each particle using the forward characteristic method. We have then in the "eulerian" description:

∂f ∂t + U 1 (X, t) .∇ X 1 f = 0 and ∂f ∂t + U 2 (X, t) .∇ X 2 f = 0 (87)
We can then introduce the characteristics of Eq. (82), which are the solutions of the dynamical system:

DX Dt = U (X(t), t) (88)
Let us denote by X(t; x, s) the solution at time t whose value is x at time s. Taking X(t) to be a solution of Eq. ( 88) we have then:

D Dt (f (X(t), t) = ∂f ∂t + Ẋ (X, t) .∇ X f = ∂f ∂t + U (X(t), t) .∇ X f = 0 (89)
which means that f is constant along their characteristics. This may also be written as

f (X(t; x, s), t) = f (X(s; x, s), s) = f (x, s) (90)
for any time t, s and phase space coordinate x. Note that if the conditions (86) are not fullfilled, the resolution of Eq. ( 83) by splitting is now equivalent to solve advective equations with a source term:

∂f ∂t + U 1 .∇ X 1 f = -f div X 1 (U 1 ) (91)
and

∂f ∂t + U 2 .∇ X 2 f = -f div X 2 (U 2 ) (92)
We may in that case introduce a cumulative error at each time step, resulting in poor density conservation irrespective of the numerical procedure for solving the time splitted equations. Thus in ref. [40] we have demonstrated that a necessary condition for the use of the famous time splitting scheme to preserve the conservative character of the Vlasov equation is that the advective fields U 1 and U 2 are both divergence free. The semi-lagrangian scheme was then extended in higher phase space dimensions by Ghizzo A., Huot F., Bertrand P. in [24].

To remedy this problem we have two solutions:

• The first is to make 2D (or even 3D) semi-lagrangian advections with the need of 2D (or 3D) interpolations which are now well in hand from the computational point of view (see Sonnendrücker et al 1998 in [42]).

• The other idea is to keep the Vlasov at the conservative level (84) and (86) allowing us to keep time splitting in 1D advections, and perform these advections in conservative form. This is the basis of another family of finite volume numerical schemes for the Vlasov equation (see, for instance, Filbet, Sonnendrücker and Bertrand (2001) in [19]).

Parametric instabilities in relativistic regime

A well-known solution of Eqs (76) to (80) is the Akhierzer-Polovin solution in [43], corresponding to the propagation of a circularly polarized transverse electromagnetic wave. In that case it is obvious that:

A 2 ⊥ = const (93)
and we define the constant normalized potential vector amplitude

a 2 0 = e 2 A 2 ⊥ (x, t) m 2 e c 2 (94)
Furthermore, it is assumed in Akhiezer-Polovin theory that the exact solution corresponds to wave propagation in a homogeneous plasma, i.e.

f (x, p x t) = n 0 F 0 (p x ) (95)
where F 0 (p x ) is a one variable function normalized to unity. It is obvious that Eq. ( 76) is satisfied with a zero longitudinal electric (plasma) field E x = 0. Eq. ( 79) may be written in the standard form

∂ 2 A ⊥ ∂t 2 -c 2 ∂ 2 A ⊥ ∂x 2 + ω 2 (x, t)A ⊥ = 0 (96)
In Eq. ( 96) ω p is now a constant. Finally A ⊥ obeys a linear homogeneous wave equation from which we can easily get the dispersion relation

ω 2 = ω 2 p γ 0 + k 2 c 2 (97)
where the Lorentz factor γ 0 was given by:

γ 2 = 1 + a 2 0 ( 98
)
From Eq. (94) a 2 0 can be related to the power momentum p osc in the standard expression:

a 2 0 = p 2 osc m 2 e c 2 = I λ 2 0 1 .368 × 10 18 ( 99
)
where I is the laser intensity given in W cm -2 units and λ 0 in µm. Thus increasing laser intensities lead to relativistic increase of the electron density and the relativistic plasma frequency ω p γ -1 . Therefore large intensity wave propagation in homogeneous plasma is thus allowed with electron density greater than the electron critical density n crit : this is the well-known induced transparency scenario (n crit = 1.1 × 10 21 λ 2 0 (µm), n crit being expressed in cm -3 ).

Now the question arises of the stability of these modes. large amplitude light (ω 0 , k 0 ) couples any existing plasma fluctuation (ω, k) to a hierachy of electromagnetic sidebands (ω 0 + lω, k 0 + lk). As seen in previous section 6, SRS (with l = -1 for Stokes and possibly l = +1 for anti-Stokes electromagnetic sidebands) is the dominant mechanism. At large intensity, the different l becomes strongly coupled with the growth of several plasma modes (ω, k), known as parametric instabilities (Guérin et al 1995 in [44] ) The simulation that we will describe now has been performed using a ratio of the electron mean density to the critical density of n 0 /n crit = 1 (n crit being usually defined as m e ω 2 0 ε 0 /e 2 ) and a quiver momentum value of a 0 = eE 0 /m e ω 0 c = √ 3 which corresponds to a laser intensity of Iλ 2 0 = 8.2 × 10 18 W cm -2 µm 2 (where λ 0 is the pump wavelength). We start the simulation with an initial electron temperature of T e = 3keV . The initial large amplitude pump wave (ω 0 , k 0 ) was taken in the form:

E = (0, E 0 cosk 0 x, E 0 sink 0 x) (100)
(ω 0 , k 0 ) satisfy the relativistic dispersion relation for circularly polarized waves (98).

Choosing k 0 c/ω p = 1/ √ 2, we obtain ω 0 = ω p , i.e. a ratio of n 0 /γ 0 n crit = 0.5., i.e. well below the relative plasma density. The choice of the wavenumber k 0 = 4π/L = 2∆k in effect determines the box length L in term of c/ω p . The corresponding plasma wave predicted by theory is then very close to 2 ∆k. Here we have solved the dispersion relation, i.e.

D + D -= ω 2 p a 2 0 2γ 3 0 k 2 c 2 D p -1 (D + + D -) (101)
where D p , D ± correspond respectively to the dispersion relation of the EPW and of the electromagnetic waves in the relativistic plasma. We have here the usual relations:

D p = ω 2 - ω 2 p γ 0 (102)
D ± = (ω ± ω 0 ) 2 -(k ± k 0 ) 2 c 2 - ω 2 p γ 0 (103)
By solving the dispersion relation in the case of a cold plasma, the growth rate of the relativistic parametric instability, induced by a high-intensity pump wave, can be evaluated as a function of the wave number; the result is shown in Fig. 12.

We see clearly that that the most unstable mode is located at kc/ω p = 1.40 and corresponds to a growth rate of Γ/ω p 0.409. In Fig. 12, we have plotted the time evolution of the most unstable plasma mode (here mode 2) on a logarithmic scale: the curve indicates a growth rate around of Γ num /ω p 0.406 in good agreement with the expected value predicted by the linear theory (for a cold plasma).

By comparison the theoretical growth rate obtained from the numerical solution of the dispersion relation (indeed here the maximum value) is Γ/ω p 0.409, in very good agreement with the numerical value obtained by the Vlasov code. But the more striking point is that we did not introduce any initial density perturbation in our simulation. The instability starts just from the round-off-errors, which for a 64-bit machine is of order 10 -15 . Indeed we can follow the growth of the unstable plasma mode over 15 decades! This result clearly demonstrates that the low noise character of the Vlasov code allows a very powerful and extremely precise study of the growth of instabilities over a large number of decades. A PIC code in comparison needs a huge number of particles to have such a low initial noise level which scales only as 1/ N part . Thus the need for a very small g pic is oviously in favour of Vlasov codes, as shown in table I.

The last point concerns the problem of the filamentation in phase space. To demonstrate the efficiency of our algorithm and the high level of complexity induced by the Vlasov equation, even in cartesian geometry, we look now at the particle dynamics in phase space. It is well known that the electromagnetic waves are parametrically unstable when they propagate through a plasma. At low intensities (a 0 1), they are clearly identified as stimulated Raman scattering (see for instance section ???). Relativistic effects add the relativistic modulational instability (RMI) and also the relativistic filamentation instability (RFI). There are two aspects of the same instability, the wave vector being, respectively, parralel and transverse to the wave vector of the pump wave (incident electromagnetic wave). Although we talk also about "filamentation", we have to consider usually here the aspect of filamentation of the laser light spot. However, as we will see, relativistic effects add also complexity in the phase space. Figures 13 show the behavior of the electron distribution function in phase space. Here we have used a phase space sampling of N x N px of 512 × 768 points and a time step of tω p = 0.01. Here the formation of smaller and smaller scales in phase space results from the physical process in duced by the relativistic parametric instability. Such a phenomena occurs even in a simple geometry (cartesian) using periodic boundary conditions and play also a major role in causal simulations. Here the source of the problem comes from the particle acceleration mechanism induced by the parametric instability, which is characterized by the tendency to generate thin filaments in phase space. At some time in the simulation, the phase space grid becomes too coarse to follow these thin filaments.

Before endind this section two remarks must be pointed out:

• Vlasov simulations presents advantages: the lack of numerical noise and the fine resolution in velocity space, provided that the dimension of velocity spaceis as low as possible. This efficiency (in terms of accuracy and computational time) is however lost when filamentation takes place in phase space which then requires an important increase of the total number of points of the phase space grid to follow this filamentation. Usual semi-Lagrangian Vlasov models are not well adapted to describe the Dirac-like distribution functions or fine phase space filaments which characterize, for instance, the particle acceleration in laser-plasma interaction or the well-known phase space filamentation in velocity of the Vlasov equation.

• Vlasov codes are powerful tools to study wave-particle interaction with interesting results for trapping and action transfer from particles and waves. Nevertheless we would point out the problem of simulating ultra high intensity lasers. Due to strong acceleration, it may happen that most of the particles are now localised in tiny regions of phase space while large regions are empty of particles. Very fine grids are needed and the number of mesh points in velocity space must be greater than the value considered in table II. For those mesh points where f is zero, the computer is spending a lot of CPU time in computing only zeros. Taking for instance N v = 1000 instead of N v = 100 (in table I) gives results which are less in favor to Vlasov codes! This is a major drawback of Vlasov methods using a uniform and fixed mesh. Since the non zero regions evolve continuously with time, there is the need to have adaptive time dependent grids with small mesh size in the non zero regions and large mesh size for the sparse regions. This would decrease N v to more interesting value for Vlasov simulation purpose. For this reason methods using an adaptive mesh must be investigated, and especially when the physical interesting phenomena occur in a localized phase-space region and do not require a fine mesh everywhere. This consideration will often make Vlasov codes preferable to PIC codes.

Adaptative method for ultra high intensity laser pulses

Hoping the reader is now convinced that Vlasov codes are powerful tools to study wave-particle interaction with interesting results for trapping and action tranfer from the plasma wave, we would point out the problem of simulating ultra high intensity lasers, i.e. a 2 0 1. In that case the whole distribution function is strongly affected and not only the particles close to the phase velocity of the plasma wave. A typical example of phase space portrait after the onset of wavebreaking during the growth of a parametric instability was shown in Fig. 13. The efficiency of Vlasov code (in term of accuracy and computational time) may be however lost when such a filamentation takes place in phase space which requires an important increase of the total number of points of the phase space grid to follow this filamentation. Usual semi-Lagrangian Vlasov models are not well adapted to describe the Dirac-like distribution functions of fine phase space filaments which characterize, for instance the particle acceleration in laser-plasma interaction or the well-known phase space filamentation in velocity of the Vlasov equation. Much progress is then expected from new adaptive time dependent grids with small mesh size in the non-zero regions and large mesh size for the sparse regions. This approach bears the promise to improve Vlasov solvers in terms of performance and accuracy.

Recently N. Besse et al (2009) [45] have developped a new method for the numerical resolution of the relativistic Vlasov-Maxwell system on a phase space grid using an adaptive semi-Lagrangian method. The authors introduce a two-dimensional phase space mesh which can be refined or derefined adaptively in time. For this purpose, hey use a technique based on wavelet multiresolution analysis, which gives a powerful and natural refinement criterion based on the local measurement of the approximation error and regularity of the distribution function. Therefore the multiscale expansion of the distribution function allows to get a sparse representation of the data and thus save memory space and CPU time. The adaptive method can be overlayed the classical semi-Lagrangian method which is based on the conservation of the distribution function along particle trajectories. The phase-space mesh can be updated using a multiresolution technique as wavelet expansion, which gives a sparse representation and a natural criterion to perform local grid refinements. Now the presented examples might be considered as a first attemps to investugate the possible use of adaptive techniques especially when the interesting physical phenomena occur in a localized phase space region and do not require a fine mesh everywhere. Thus the adaptive method can be overlaid to the classical semi-lagrangian method which is based on the conservation of the distribution function along particle trajectories. The phase space mesh can be updated using the multiresolution technique as wavelet expansion, which gives a sparse representation and a natural criterion to perform local grid refinements.

We have seen , in previous section, that Vlasov simulations provide an excellent description of the small scales of the phase space dynamics, although the filamentation mechanism or phase-space mixing are saturated by the numerical dissipation of the numerical scheme. Indeed, there exists a time when the phase space grid becomes too coarse to follow these thin filaments which get finer and finer as time goes on. Furthermore this mechanism is electromagnetic strongly amplified by relativistic effects since parametric instabilities (as RMI or RFI) tend naturally to produce thin filaments via strong particle acceleration. Thus the study of relativistic filamentation (see previous discussion of velocity space filamentation in sec. 3) provides a natural test for checking the accuracy of the adaptive scheme. Fig. 14 shows the distribution function (on the left) and the corresponding relative adaptive phase space grid (on the right). we see that the total domain of computation becomes to be filled by refinement points, thus the adaptive method becomes as competitive as a numerical method where the phase space is discretized with a fixed mesh. Physical parameters are kept identical in that simulation (in comparison to fig. 12 and13).

In this test case we are interested in the laser wake-field process, which seems to be a good way to accelerate particles to very relativistic velocity.

In the simulation the circularly polarized laser pulse has a quiver momentum

d v = 1 d v = 2 d v = 3 g pic = 10 -2
10 10 4 10 7 g pic = 10 -4 10 -1 10 2 10 5 g pic = 10 -6 10 -3 1 10 3

Table 2: PIC code versus Vlasov code : ratio N vlas /N part for N v = 1000 a 0 = 3/2 and the density of the plasma is such that n 0 /n c = 0.1. The phase space is sampled with a coarse grid N x N px = 1024 × 256 and with the possibility to refine the mesh adaptively on three more levels in each direction. In Fig. 15 the peak of the distribution function is typically associated to the wavebreaking process, which is followed by the acceleration of the electrons. We see clearly a strong acceleration of the plasma electrons in the laser pulse wake-field, up to moment value p x /mc = 30. We see that the electrons are well accelerated along a straight line as can be also reproduced by particle method. Therefore our adaptive scheme can depict the same thin filaments as PIC codes do, whereas it could be almost impossible with a fixed mesh of phase space because of the huge needs of CPU time and memory capacity. Then the advantage of the adaptivity compared to classical Eulerian numerical method constructed on a fixed mesh is that you can observe and follow in time the fine structures like these narrow localized filaments of accelareted particles with a reduced computation cost and with weak memory storage. Let us note that in Fig. 15 the whole region (white and red areas) corresponds to the finest grid we would have to compute with a non adaptive Eulerian numerical method. Several remarks must be pointed out:

• An interesting point is the existence of dynamical invariant (exact generalized canonical momenta). If we take these invariants into account as new independant variables, then we can discretize them coarsely, since there is no exchange of particles between the different populations. Providing a special choice of initialconditions, they remain exact for all time t and space x. From both an analytical and numerical point of view, the use of excat invariants allows obviously a reduction of the complexity of dealing with a full Vlasov equation with a 3D velocity space. This makes the use of an adaptive mesh scheme a very powerful method when it is coupled with such invariants (transverse canonical invariants met in laser-plasma interaction or adiabatic invariants in gyrokinetic models). The performance and accuracy obtained in the axamples presented here are very satisfying and these examples show that our adaptive scheme is able to describe accurately nonlinear effects such as the filamentation process of the parametric instabilities in the relativistic regime.

• Even though the use of an adaptive mesh in semi-lagrangian Vlasov solver will eventuelly be of benefit as well in terms of memory as in computational cost for many problems in plasma physics and in beam physics, the developement and implementation of an efficient adaptive solver is a complex problem that will still need much collaborative effort between physicists, applied mathematicians and computer scientists. Indeed, in the problems that were addressed in last sections, as in many other problems involving the Vlasov equation, the number of grid points that remain with our wavelet adaptive technique after some time becomes too high for the method to be competitive with a uniform mesh sover due to the importance overhead associated to the wavelet method which gives us the adaptive method. This difficulty becomes even more important for a parallel version of the code. There are now efficient uniform grid Vlasov solvers that scale very well on thousands of processors whereas the efficient parallelization of adaptive mesh solvers is a real challenge.

• Vlasov simulations were performed in a two-dimensional phase space where the developement of thin filaments, strongly amplified by relativistic effects, is the signature of strong acceleration in laser-plasma interaction. The phase space filamentation requires an important increase of the total number of points of the phase-space grid as the filaments get finer as time goes on. The strcture of these filaments is essentially one-dimensional even in higher dimensions. Therefore, the adaptive method could be more useful and efficient in cases where these thin filaments that need to be resolved are a very small fraction of the hyper-volume, which arises in higher dimensions because of the surface-tovolume scaling and the essentially one-dimensional structure of the filaments.

• The main drawback of our scheme is the poor local character -in phase spaceof the semi-lagrangian method and the wavelet multiscale reconstruction that we have used. Indeed the lack of space locality of the semi-lagrangian method comes from the characteristic curves tracking which prevents us tio know a priori the size of the dependance domain of each mest point and thus the number of cells that the characteristics will cross. The poor locality of interpolets comes from the size of its support. This poor locality of the scheme involves more points than it is needed. The second one is that it is difficult to construct an efficient parralel algorithm such as domain decomposition algorithm. Therefore, the main way to improve the efficiency of the adaptive method is then to increase the local character in phase space of the numerical scheme, by considering multiscale reconstruction with more compact support and by replacinf the semi-lagrangian method with more local -in space-numerical schemes as compact finite difference schemes or finite element residual schemes which are well suited for parallel domain decomposition techniques.

8 Conclusion and outlook.

The numerical simulation is a huge topic with hundreds of papers published every year. Eulerians Vlasov simulations are slowly introduced in place of the well-known Lagrangian PIC models for two main reasons: the lack of numerical noise and the very good resolution of the distribution in phase space, provided the dimension of the momentum space is as low as possible. The first property makes Vlasov codes a powerful tool to study the growth of instabilities of laser-plasma interaction even for 2D spatial problems. The very good resolution in phase space allows a precise study of wave-particle interaction including particle trapping and particle acceleration.

Another advantage of the Eulerian Vlasov code is that they incorporate easily ideas issued from the analytical results of the Vlasov equation. For example we have the possibility to reduce the phase space dimension (and therefore the computational cost) by using exact invariants (as the canonical momenta in the study of laserplasma interaction, or Water-Bag model or even adiabatic invariants in gyrokinetic modelling when well separated different time scales exist for a given problem). Such a technique of adiabatic invariants have allowed to investigate the ion temperature gradient (ITG) instabilities or (the famous trapped ion modes) by Depret et al in [46] in action-angle variable well-known in the study of the turbulence in tokamaks. We will sample consequently with less points these invariants since they appear as numerical parameters.

What are the problems we must now concentrate on? We see physical approaches in basically two directions

• The first point concern the gyrokinetic modeling. These studies are crucial for the study of the anomalous transport and plasma confinement in magnetic fusion. The fact that the gyro-averaged distribution function depends only on one parallel velocity component v makes a Vlasov approach a rival method with the particle approach. The semi-lagrangian was found to be an alternative method to solve the drift-kinetic Vlasov equation and allows to compute the distribution function with moderate dissipation. Nevertheless capturing the key elements of more realistic higher dimensionality simulations into simplified modules remains a formidable task. Significant effort in gyrokinetic modeling is the use of Gyro-Multiple Water Bag approach (see refs. for more details [47]) which introduces a special class of initial conditions allowing to reduce the full kinetic Vlasov equation into a set of hydrodynamical equations.

• A final remark concerns the use of adaptative mesh refinement approaches: for some problem encountered for instance in ultra high intensity laser pulse as in laser-Wakefield accelerator, we have to simulate large accelerated particle beams located in very thin regions of phase space while large regions are empty of particles. Similar behavior can be observed for the ion distribution function in phase space in the case of ion acceleration mechanism in overdense plasmas. Since non zero regions (and in fact fine phase space filaments or δ-Dirac-like function) evolve continuously with time, there is the need to have time-dependent adaptative grid with small mesh size in the nonzero regions and larger mesh size for the sparse regions. Much work is need but preliminary results based on wavelet multi-resolution technique are very encouraging.

• The preceding examples illustrate the challenges in modelling wave-particle interactions in laser-plasma interaction. High accuracy and resolution are required to correctly model such a trapped or accelerated particle dynamics.

The challenging nature of plasma physics in general and fusion research in particular leads to establish an interdisciplinarity research that targets the development of capabilities that "bridge" various areas of plasma physics together with computer science and applied mathematics.

Figure 1 :

 1 Figure 1: Phase space representation for nonlinear Landau damping with kλ D =0.3 and α = 0.1.We see on the three first pictures the overtaking and the subsequent breaking of a "phase space wave", followed by the formation of an BGK-like hole.

Figure 2 :

 2 Figure 2: Phase space representation for nonlinear Landau damping with kλ D = 0.25 and α = 0.05 (initial perturbation on mode k 0 , and β = 2(corresponding to the main perturbation on mode 2k 0 . The simulation shows the occuring of a pairwise coalescence leading to a single hole, similar to a KEEN wave with a frequency well below the plasma frequency. The curves show the beginning of the simulation. ( see thesis of P. Mineau, Résolution numérique du système Vlasov-Poisson, application à des problèmes de l'astrophysique et de la physique des plasmas, Orléans 1989).

Figure 3 :

 3 Figure 3: Following the results of Fig. 2, phase space representation for nonlinear Landau damping with kλ D = 0.25 and α = 0.05 (initial perturbation on mode k 0 , and β = 2(corresponding to the main perturbation on mode 2k 0 . The simulation shows the pairwise vortex merging leading to the formation of a single electron hole in phase space. Note that the size of the merging structure is small. (see thesis of P. Mineau, Résolution numérique du système Vlasov-Poisson, application à des problèmes de l'astrophysique et de la physique des plasmas, Orléans 1989).

Figure 4 :

 4 Figure 4: Successful KEEN results , with the trapping time τ B (shown by the doubleheaded arrow) being about equal to the drive duration. The result of the drive is to produced a long-lived KEEN wave. As often happens with KEEN wave drive, the final frequency is slightly reduced from the drive frequency.

Figure 5 :

 5 Figure 5: Phase space images of electron velocity distribution function during the drive-up period (t ≤ τ B ) and after at the same time (namely t/τ B = 1 2 , 1, 2, 10). The total-force potential profiles for each time frame are superposed for convenience.

  e = S 0 + S e = C e (t = 0) = const (60) which are nothing else but the conservation of quasi-particles (i.e. the total number of photons and the number of photons plus plasmons conservation).

Figure 6 :

 6 Figure 6: For a periodic SRS-F case, time evolution of the different action densities: electromagnetic pump S 0 , scattered Stokes mode S s and their sum C s on the top panel. The action is transfered back and forth between the pump and daughter waves in a classic fashion.

  . The panel shows the time evolution of the pump S 0 , the Stokes mode S s and their mutual sum C s . The corresponding curves for S e and their sum C e are shown in the bottom panel in Fig. 7 during the begining of plasma evolution (one electromagnetic cycle). Since the velocities are normalized to the light velocity c and frequencies to the plasma frequency ω p , the choice of k 0 (here equal to k s + k e ) determines the plasma box length in term of cω -1 p . The pump frequency is close to ω 0 = 2.6ω p and therefore the density as compared to the critical density is then n 0 /n c 0.15. With these parameters a good frequency match was obtained by choosing k 0 c/ω p = 2.40, i.e. a box length of L 5.23cω -1 p . We have then for the EPW k e c/ω p = 1.20 or equivalently in Debye wavelength k e λ D 0.205. The normalized pump electric field amplitude is eE 0 /m e ω p c = 0.28 .

Figure 9 :

 9 Figure 9: We present here the separatrices calculated from the simulation electric potential superimposed on the very accurate phase space representation afforded by the Vlasov code. Here color shading has been used to indicate the relative values or normalized phase space density between 10 -3 and 10 -4 . A striking agreement can be seen between the separatrices and the phase space trapping limits for the growing wave.

Figure 10 :

 10 Figure 10: For the parabolic profile and for the study of SRS-B in optical mixing, phase space plots of a given cell shown at time tω 0 = 5600. The first image clearly exhibits the formation of very neat phase space hole wave train induced by particle trapping and SRS-B. After this first frame the plasma begins the bursting or turbulent regime. In the middle frame considerable vortex merging has taken place. The last image shows the formation of new vortices resulting from the coalescence and characterized by a wavenumber close to k V M c/ω 0 ∼ 1 observed in the k-spectrum.

Figure 11 :

 11 Figure 11: Growthrate of the relativistic parametric instability Γ/ω p induced by an ultra intense circularly polarized electromagnetic wave versus the wavenumber kc/ω p . The parameters are a 0 = √ 3 and n 0 /n crit = 1.We see clearly that the most unstable mode is now located at kc/ω p = 1.4 and corresponds to a maximaum growthrate of Γ/ω p 0.409.

Figure 13 :

 13 Figure 13: Phase space representation of the electron distribution function obtained usinh an 1D semi-lagrangian full relativistic (without time splitting) Vlasov model. The parameters are a 0 = √ 3 and n 0 /n crit = 1.

Figure 14 :

 14 Figure 14: Phase space representation of the electron distribution function obtained using an 1D wavelet-MRA-based adaptive semi-lagrangian Vlasov model. The physical example correspond to the relativistic parametric instability. The parameters are a 0 = √ 3 and n 0 /n crit = 1. We have represented the electron distribution function (at left) and the adaptive numerical grid (at right).

Figure 15 :

 15 Figure 15: Phase space representation of the electron distribution function obtained using an 1D wavelet-MRA-based adaptive semi-lagrangian Vlasov model. The physical example correspond to the Laser Wave Field Accelerator process. We have represented the electron distribution function (at left) and the adaptive numerical grid (at right).

and first shown in the thesis of P. Mineau, Résolution numérique du système Vlasov-Poisson, application à des problèmes de l'astrophysique et de la physique des plasmas, Orléans 1989.

Acknowledgments:

The authors acknowledge the IDRIS (Institut du Développement et des Ressources pour l'Informatique Scientifique, CNRS, Orsay, France) for computer time allocation on their supercomputers.