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 individual 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.

As a matter of fact it was quickly realized that a numerical solution of the Vlasov equation itself was more appropriate. A first solution through a spectral (Fourier) method in space and Hermite functions in velocity space was proposed by Feix and Grant in [1]. More recently Schumer-Holloway in [2], S. Le Bourdiec et al [3] 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. The Hermite basis is a natural choice for Maxwellian-like distribution function because the lowest order Hermite function is a Gaussian one. However the Hermite development exhibits a slow convergence.

Nevertheless recent developments using the splitting have attempted to address this convergence problem. A first such an approach of the Vlasov-Poisson equation was given by Cheng and Knorr in ref. [4] consisting in splitting up the Vlasov equation into two partial derivative equations, one in the x, t, the other in p, t; where the symbols have their conventional meaning. The method is the basis for modern Vlasov codes.

Considering the dimensionless form of the one-dimensional (1D) Vlasov-Poisson equations:

∂f ∂t + p ∂f ∂x + E(x, t) ∂f ∂p = 0 (1)
∂E ∂x = +∞ -∞ f (x, p, t) dp -1 (2)
The method introduced two ideas: (i) the time splitting which consists in splitting up the free-streaming term and the acceleration term in (1). The integration of Eqs. ( 1) and ( 2) is reduced to the following shifting sequence of the distribution function:

1. Compute f * (x, p) = f n (x -p∆t/2, p), f n (x, p) being the distribution function at time t n = n∆t.

2. Compute the electric field E * by substituting f * into Eq. (2).

3. Compute f * * (x, p) = f * (x, p -E * ∆t) 4. Repeat again step 1 to obtain f n+1 (x, p) = f * * (x -p∆t/2, p) (ii) the reconstruction of f * , f * * and f n+1 on a mesh in phase space. Cubic spline interpolation have been extensively used (see for instance Ghizzo et al in [5] and references therein). Fourier interpolation have also been developed (see for more details Ghizzo et al in [6]). Another scheme for the Vlasov equation is the flux balance method (FBM) (see ref. [7] and more recently the positive and flux conservative PFC method (of Filbet et al in ref. [8]). The basic idea is to compute the average of the Vlasov equation solution in each cell of the phase space grid by a conservative method. Unlike classical Eulerian algorithms such as finite difference or finite volume schemes, these methods are not restricted by the Courant-Friedrich-Levy (CFL) condition on the time step. The PFC method enforces the conservation of the global and controls numerical oscillations of f . For the Vlasov-Poisson system, this scheme is equivalent to an integration of the Vlasov equation along the characteristics and is correct up to the second order in time step. For a recent extensive review on the analytical and numerical methods to solve the Vlasov equation see the review paper of M.R. Feix et al in ref. [9].

During the '80s, PIC codes were able to perform 2D simulations, getting useful information out of them and direct Vlasov solvers lost their interest for most people due to their high numerical cost. However the Vlasov equation coupled with the Poisson's equation or Maxwell system often contains filamentation, which has been one of the reasons why Vlasov solvers have been poorly considered, compared to PIC solvers which are insensive to this problem. The different attempts on the numerical solution of the Vlasov equation (especially with spectral methods) were not able to solve this problem of filamentation in velocity space unless using a huge number of Hermite polynomials or Fourier modes. The source of this problem lies in the treatment of the free streaming term of the Vlasov equation. A Fourier transform on the x-axis gives:

∂ f ∂t -ikp f (k, p, t) = 0 (3)
leading to a solution in the usual form:

f (k, p, t) = f (k, p, 0) e ikpt (4)
We see in (4) that, as a function of momentum p, f (k, p, t) oscillates at the frequency kt. When this frequency reaches the inverse of the size of the velocity cell ∆p, we cannot follow the exact solution of f . Indeed the distribution function f (k, p, t) is constant along the characteristic curves, which become close, so that phase space regions where f (x, p, t) has different values come close together and steep gradients are thus generated. Furthermore in laser-plasma interaction, this mechanism is strongly amplified by relativistic effects since electron relativistic parametric instabilities have already this tendency to generate thin filaments in phase space (in particular during particle acceleration processes, see for instance ref. [START_REF] Ghizzo | [END_REF]). Even with the self consistent longitudinal field, the whole Vlasov equation exhibits small scale oscillations in phase space, for instance due to the trapping of particle trajectories around the phase velocity of the wave phenomenon under consideration.

The Cheng and Knorr method in ref. [4] proved to work very well in the electrostatic case and in the case of the weak relativistic regime of the laser-plasma interaction. Direct solution of the Vlasov equation itself on the phase space grid have been found to be a powerful tool to study in detail the particle dynamics in stimulated Raman scattering (referred here as SRS), both in a conveniently short but idealized periodic plasma (see ref. [10]) and in the more realistic open system (see refs. [11,12,13]) with external source. These codes allow precise measures of action transfer between photons, plasmons and electrons in Raman scattering processes. These Vlasov simulations of plasma have revealed a rich variety of phenomena associated with the trapped particle dynamics induced by SRS or by a beatwave experiment. However these phenomena become very burdensome to treat numerically when electromagnetic frequencies are large compared with the plasma frequency ω p , as met in real experiment (ω 0 /ω p ≃ 33 in ref. [14] or ω 0 /ω p ≃ 100 in ref. [15] where ω 0 is the usual pump laser frequency and ω p the electron plasma frequency). We have demonstrated in refs. [16,17] an interesting possibility to extend such detailed modeling to frequency ratios greater than the current practical maximum of 10 or so for Vlasov or PIC codes, by using mode equations for the electromagnetic waves rather than the full electromagnetic Vlasov codes. The resulting hybrid version which we have called the Hilbert-Vlasov code thus provides a saving of order (ω 0 /ω p ) 2 in computer time for forward Raman-like processes, as compared with a direct attack by the full electromagnetic Maxwell-Vlasov scheme with the highest space-time resolution required.

The advantage is the semi-lagrangian Vlasov scheme can use large time steps than explicit Eulerian ones, is noiseless in comparison to PIC codes, but the price to pay is to reconstruct a regular grid using interpolations. However the Cheng and Knorr scheme could not be applied as easily to more complicated Vlasov problem, for instance the gyrokinetic Vlasov equation used in Tokamak plasmas or the full relativistic Vlasov equation in which the Lorentz factor couples for instance the different components of momentum p. The original algorithm of Cheng and Knorr has been cast into the more general framework of semi-lagrangian method by Sonnendrucker et al in ref. [18] and Nakamura and Yabe in ref. [19]. In particular it would appear that direct integration using 2D full advection (without time splitting) is necessary in the strong regime of the laser-plasma interaction (referred here as LPI) for keeping accurate solutions of the Vlasov-Maxwell system. The generalization to the relativistic case is presented in ref. [20] We present here 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 chapter 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 in LPI.

The structure of the paper is organized as follows. Part 0.2 is devoted to phase space properties and the difference between the conservative and advective forms of the Vlasov equation. In part 0.3 we used the exact invariance of generalized transverse momenta to built a Vlasov-Maxwell model which was found to be pertinent in LPI simulations. Then we show in section 0.4 how the existence of such invariants allows us to reduce the computational CPU time and points out the efficiency of one-dimensional Vlasov codes in the study of "kinetic" regime recently observed in LPI. Section 0.4.1 is devoted in Manley-Rowe partition between photons or plasmons and its connection with electromagnetic Vlasov models and in particular on the importance of resonant wave-particle interaction and the use of full kinetic (Vlasov) model. Section 0.4.2 gives another example of kinetic regime met in LPI concerning backward-SRS in optical mixing and the importance of a precise description of trapped particle dynamics at low laser intensities.

Intensities above 10 10 W cm -2 are now reached with recently developed chirped pulsed laser and relativistic plasma wave-breaking can now be experimentally investigated as well as the self-induced transparency (i.e. laser propagation through a so-called overdense plasma layer in which the laser propagation is classically forbidden). We will come back to these effects and the possibility to use Vlasov code in the relativistic regime in part 0.5. Sec. 0.5.1 is devoted to the usual hole boring scenario which is used here to check the model and to show the difference met in the self-induced transparency mechanism leading to a strong (anomalous) absorption which is treated in sec. 0.5.2. Sections. 0.5.3 and 0.5.4 concern the excitation of kinetic nonlinear low-frequency waves which may be accurately described by such a Vlasov code and which correspond to curious electron phase space anomaly driven by a ponderomotive force. Finally our conclusions are presented in part. 0.6.

The Vlasov model.

The traditional approach to plasma kinetic theory, in the case of electrostatic systems, starts from the Vlasov-Poisson system (described by Eqs. ( 1) and ( 2)) coupling the particle distribution function f (x, p, t) with the electric field E (x, t).

To simplify the notation and without any loss of generality we supposed an onedimensional space. 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 = dxdp of a close contour defined at time t 0 . At a later time t the particles on the initial contour have moved and define a new contour, the volume of which is dxdp = dS. Using the Hamilton's equations, i.e.

ẋ = ∂H ∂p and ṗ = - ∂H ∂x (5)
we have then dS = dS 0 . This is this property which is used to built the Water-Bag (WB) model. Introduced initially by De Packh in ref. [21], then by Feix, Hohl and Bertrand in refs. [START_REF] Feix | A minimum energy principles in plasma physics and stellar dynamics[END_REF]23], the water-bag model was shown to bring the bridge between fluid and kinetic description of a collisionless plasma, allowing to keep the kinetic aspect of the problem but with the complexity as the fluid model. The extension to gyrokinetic modeling is recently realized in refs. [24,25]. On the other hand, in laser-plasma interaction the WB model allows us to consider the effect of a finite temperature on the growthrate of the Stimulated Raman Scattering (SRS)-like instability in the relativistic regime, in particular the fact that a hot electron temperature may strongly reduce the growthrate of the such an instability.

The second fundamental property is that particles which at time t 0 , are inside the contour C 0 are exactly the same which, at time t are inside the contour C. We can then define phase space densities 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 have then f (x, p, t) = f (x 0 , p 0 , t 0 ). We deduce

Df Dt = 0 (6)
in the limit where dS 0 and dS tend to zero. The kinetic limit reads as the conservation law for the distribution function f , allowing to write the Vlasov equation in a conservative form:

∂f ∂t + ∂ ∂x ( ẋf ) + ∂ ∂p (E(x, t)f ) = 0 (7) As ∂ x (p) + ∂ p (E) = 0 (8)
(in its dimensionless expression) and since the plasma is a Hamiltonian system, the conservation law (7) may also be written as the transport equation ( 1). An important concept in plasma physics is that we focus on collective behaviors of the particles, which are not very sensitive to the grain character (particle discreteness). In 1960 in ref. [26] Rostoker and Rosenbluth introduced the "virtual" dichotomy experiment consisting in "fragmenting" particle of mass m e and charge e into many smaller particles, keeping the total charge (per volume unit) in the system n 0 e or the ratio e/m e constant. In particular, the plasma frequency ω p and the electron Debye length λ D are invariant in this transformation. In this Vlasov limit, one discards the granular nature of the plasma. This f conservation property is used in the time-splitting property in semi-lagrangian scheme for solving the Vlasov equation.

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

∂f ∂t + div X (U(X,t)f ) = 0 (9)
where X stands for the phase space coordinates and U is a divergence-free advection field. Splitting the components of X into two sets X 1 and X 2 , Eq. ( 9) 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 (10)
Note that since U is divergence free, Eq. ( 9) can also be written in an advective (transport) form:

∂f ∂t + U (X, t) .∇ X f = 0 (11)
which is equivalent in 1D to Eq. ( 1). Moreover, it is well-known (see for instance ref. [START_REF] Marchuk | Meth. of Num. Math[END_REF]) that solving separately

∂f ∂t + div X 1 (U 1 (X 1 , X 2 ,t)f ) = 0 (12) ∂f ∂t + div X 2 (U 2 (X 1 , X 2 ,t)f ) = 0 (13)
keeps the second order accuracy for the whole Eqs. (9) and equivalently (10) 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 (14)
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 (15)
We can then introduce the characteristics of Eq. (11), which are the solutions of the dynamical system:

DX Dt = U (X(t), t) (16)
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. ( 16) we have then:

D Dt (f (X(t), t) = ∂f ∂t + Ẋ (X, t) .∇ X f = ∂f ∂t + U (X(t), t) .∇ X f = 0 (17)
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) (18)
for any time t, s and phase space coordinate x. Note that if the conditions (14) are not fullfilled, the resolution of Eq. ( 10) 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 ) (19)
and

∂f ∂t + U 2 .∇ X 2 f = -f div X 2 (U 2 ) (20)
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. [20] 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 in ref. [START_REF] Ghizzo | [END_REF].

Vlasov models for laser-plasma interaction

Vlasov models have long been used to study short pulse high intensity laser-plasma interaction where collisions can be ignored. This is the case for parametric instabilities, beat-wave, Raman, Brillouin scattering, particle acceleration mechanism or self-induced transparency for overdense plasmas. To describe the behavior of an 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 et al in ref. [20], 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 p = 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 (21)
where the Lorentz factor has the conventional form:

Γ 2 = 1 + p 2 m 2 e c 2 (22)
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) (23)
with

Γ = 1 + (P c -eA) 2 m 2 e c 2 1 2 (24)
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 (25)
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 (26)
and in the transverse direction

DP c⊥ Dt = -∇ ⊥ H = 0 (27)
This last equation ( 27) 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))] (28)
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. The family of fluid closures, based on the Water-Bag reduction of the Vlasov equation (see refs. [START_REF] Feix | A minimum energy principles in plasma physics and stellar dynamics[END_REF]23]), constitutes another example of reduction technique. This WB model characterized by a special class of initial conditions, is strictly equivalent to the Vlasov equation for all wavelengths and all degrees of nonlinearity. 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) (29)
we obtain for the reduced Vlasov equation, using the Liouville's theorem:

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

∂f ∂t + p x m e γ ∂f ∂x + eE x - e  m e γ ∂A 2 ⊥ ∂x ∂f ∂p x = 0 (31)
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 (32)
Eq. ( 31) 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 (33)
and by the vector potential equation

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

J ⊥ (x, t) = - e 2 A ⊥ m e +∞ -∞ f γ dp x (35)
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 in ref. [29].

Parametric SRS instability in "kinetic" regime

In a plasma, the parametric instability refers usually to the decay of a large amplitude wave to other waves or modes. The simplest parametric instabilities involve an interaction among three waves: the large-amplitude pump wave and two decay (or daughter) waves. Stimulated Raman Scattering (SRS) is an example of such an instability involving the stimulated decay of an incident electromagnetic wave (so called "pump wave") into an another electromagnetic scattered wave (the Stokes wave) and a high-frequency electron plasma wave (EPW). Unlike regular SRS the phenomenon in plasmas is related to stimulated scattering from optical photons requiring frequency/wavenumber matching between the three considered waves. It should more properly be called SPWS for Stimulated Plasma Wave Scattering since an EPW may be excited by resonance. The driven version of this process is the Plasma Beat-Wave process which involves the beat of two injected laser pulses into the plasma.

It is well known that de-tuning or the loss of resonance (matching conditions cease to be valid) may lead to the saturation of the instability. Let us consider the formation of an EPW of wavevector k e and we note by λ D the (electron) Debye wavelength. In the weak Landau regime (k e λ D 0.2), another mechanism which may be responsible for the saturation of SRS-B (SRS in the backward direction, sometimes called BSRS) is the Langmuir Decay Instability (LDI), but the evidence that LDI limits the growth of SRS-B was only gathered in that regime. In the regime where k e λ D 0.3 (i.e. at high plasma temperature or low density), referred here as the kinetic regime of the instability, the saturation proves to be more complex. Different explanations for this regime were proposed. Vu et al in ref. [30] propose a classic saturation of SRS-B by de-tuning due to the nonlinear frequency shift of the EPW while Estabrook et al (see for instance ref. [31]) have proposed that SRS saturates as a result of increased nonlinear Landau damping induced by the production of hot electrons during the particle trapping process. Recently Brunner and Valeo in ref. [32] have suggested the breakup of the EPW through the trappedparticle instability using Vlasov simulations. Another mechanism might explain the saturation of SRS by the excitation of (strong) Langmuir turbulence (see for instance ref. [34]). These works show the importance of kinetic processes: particle trapping or hot electron generation which may modify the distribution function in momentum space.

The effect of kinetic processes on Langmuir wave turbulence was studied in ref. [33] using a reduced-description Particle-In-Cell (RPIC) code. SRS-B was observed to occur intermittently, in spatial-temporal pulses of activity which may extend over several hundred Langmuir wavelengths. This bursting behavior of SRS-B appears as a common feature of different nonlinear kinetic simulations. SRS-B reflectivity presents such a bursting behavior, building up and breaking down in cycles. Thus in these models electron trapping seems to play a major role, so the nonlinear behavior of the plasma is dominated by kinetic effects. A possible transition from fluid to kinetic plasma behavior was recently identified by Kline et al and presented in ref. [35]. It becomes clear that in that new kinetic regime, trapping effects and the associated nonlinear shift in frequency are expected to dominate since LDI is too heavily Landau damped to compete. Not only electron trapping modifies the velocity distribution which may result in an effective beam population of low density, but leads to a different plasma dynamics involving new low-frequency plasma waves of acoustic-like nature different from the conventional Langmuir high-frequency waves.

Another clear indication of kinetic behavior and the modification of the particle dynamics was the recent observations shown in ref. [37] at the Trident laser facility of SRS-B experiments, associated which what was termed Stimulated Electron Acoustic Scattering (SEAS), a novel scattering apparently involving a so-called Trapped Electron Acoustic Wave (TEAW) whose phase velocity v ϕ = 1.3v th is between an EPW and an Ion Acoustic Wave (IAW) (and v th being the electron thermal velocity). This new wave differs from the classic low-frequency electron acoustic mode introduced by S.P. Gary et al in ref. [61] which may take place in a two-temperature electron plasma and which plays the same role for electron plasma than IAWs in a plasma of both electron and ions. Here these TEAW are characterized by only one temperature electron component and are driven in a nonlinear way till particle trapping occurs. We refers here these new plasma wave ( usually denoted in plasma community by the same EAW denomination) as (Trapped Electron Acoustic Wave) TEAW in order to make the difference with the standard electron acoustic modes. The experimental results showing a small narrow signal with an uncharacteristic frequency shift of about 0.37ω p which could not be attributed to the usual plasma SRS (also seen but at much higher level ∼ 3000 times stronger) with its frequency shift about ω p . Such TEAWs are justified as the small-amplitude limit of waves which depend essentially on trapped electrons as theoretically predicted by Krapchev and Ram in ref. [38] and Holloway and Dorning in ref. [39] and are quite distinct from usual IAWs, which are linear waves. TEAWs are undamped electron wave absent from fluid description of the plasmas. TEAWs require a non Maxwellian electron velocity distribution, flattened at the phase velocity of the wave, without which it would be critically damped.

Let us consider the well-known linear analysis of a small perturbation in a homogeneous plasmas (with a normalized homogeneous equilibrium distribution function F 0 (v) for the electrons and a motionless neutralizing ion background) through the introduction of the plasma dielectric function ε D (ω, k) where:

ε D (ω, k) = 1 - ω 2 p k 2 +∞ -∞ ∂F 0 ∂v dv v -v ϕ = 0 (36)
and where ω p is the usual plasma frequency and v ϕ = ω/k the phase velocity and v = p/m e the particle velocity. The integral in the Landau dispersion relation may be written in the following form:

+∞ -∞ ∂F 0 ∂v dv v -v ϕ = P P +∞ -∞ ∂F 0 ∂v dv v -v ϕ + iπ ∂F 0 ∂v v=vϕ (37)
and where PP means integral in principal part. This dispersion relation leads to a branch of waves which are essentially Langmuir waves. However for TEAWs, which are another branch at low frequency and long wavelength, the trapped particle distribution effectively makes the distribution flat at velocity close to v ≃ v ϕ (i.e.

(∂F 0 /∂v) v=vϕ = 0). Recently a direct Vlasov solver was used by N.S. Sircombe et al (in ref. [START_REF] Sircombe | [END_REF]) to address fundamental aspects of the existence and stability of these TEAWs. The implications for laser-plasma interaction (relevant to the Trident experiment) were examined in detail, including LDI and SEAS. The concept of these TEAWs are a novel version of a Bernstein-Greene-Kruskal (BGK) nonlinear waves first introduced in ref. [40], 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 low-frequency kinetic waves.

The TEAWs are very similar to the kinetic electrostatic electron nonlinear (KEEN) waves. The KEEN waves have been studied in periodic Vlasov-Poisson simulations in ref. [START_REF] Afeyan | proceeding of the third International Conference on Inertial Fusion Sciences and Applications M034[END_REF] ponderomotively driven for a short period in time to form self-consistently. 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 IAW and Langmuir wave) were all legitimately sustainable as well. In the case of KEEN waves, no initial flattened electron velocity distribution function need to be invoked. KEEN waves are quite non resonant easy to drive up over a significant range of frequencies for a given wave number.

Recently, by means of PIC simulations the stability of TEAWs against decay was investigated in ref. [43]. In particular the possibility of decay of TEAWs to a longer wavelength mode was shown. Such a behavior is a typical feature met in the study of the marginal stability of BGK waves. Marginal stability analysis indicates that such BGK equilibria are indeed unstable when the system contains at least two phase space holes. The tendency to BGK structures to merge was studied using Vlasov simulations in [44]. In particular those electron holes are closely related to the "clumps" of Dupree that appear in weak turbulence theory (see for instance ref. [45]), although BGK instabilities give rise to strongly nonlinear effects. Recently the influence of low-frequency nonlinear BGK-type waves induced by trapped electrons was observed in a more complex situation in ref. [46] in SRS-B in optical mixing. In particular the saturation of Raman backscattering begins with phase space hole merging followed by a transition regime to lower wavevectors resembling weak turbulence.

Resonant wave-particle Vlasov simulations for SRS:

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 (38)
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. Thus as a first approach of the problem, we replace Eq. ( 31) by its corresponding fluid reduction, assuming the distribution function is close to a Dirac one. 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 (32) 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 (39)
∂ 2 A ⊥ ∂t 2 + ω 2 p A ⊥ -c 2 ∂ 2 A ⊥ ∂x 2 = eA ⊥ m e ∂E x ∂x (40)
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. ( 39) and ( 40) 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. (41)
E x (x, t) = 1 2 εE e (x 1 , t 1 )e i(kex 0 -ωet 0 ) + c.c. ( 42
)
together with the following expansion where ε → 0:

∂ t = ∂ t 0 + ε∂ t 1 and ∂ x = ∂ x 0 + ε∂ x 1 (43)
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 (44)
∂ ∂t + v gs ∂ ∂x a s = Λa 0 a * e ( 45
)
∂ ∂t + v ge ∂ ∂x + η a e = Λa 0 a * s ( 46
)
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. ( 44) to ( 46) over a box of length L yields to the usual Manley-Rowe relations:

D Dt (N 0 + N s ) = F 0 + F s (47) D Dt (N 0 + N e ) + 2ηN e = F 0 (48)
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. ( 47) and ( 48) become when η = 0:

C s = S 0 + S s = C s (t = 0) = const (49)
C e = S 0 + S e = C e (t = 0) = const (50)
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).

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. 1. The top 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 0 , S e and their sum C e are shown in the bottom panel in Fig. 1 till a time of tω p = 4000. 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 plasma temperature is here T e = 15keV and 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.14 which gives a quiver momentum of a osc = 0.0538. This corresponds to an intensity of Iλ 2 0 = 3, 96×10 15 W cm -2 µm 2 (λ 0 being the pump wavelength given in µm). We have used a grid sampling of N x N px of 128 × 512 in phase space and a time step of ∆tω p = 0.01. Here no hot electron population was added initially. The well conservation of the quantity S 0 + S e shows that the regime of the SRS-instability is fluid and no electron trapping is observed in phase space. 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 now 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 waveparticle interaction around the phase velocity in phase space. We choose to keep identical the others physical parameters in a first step. However the phase space sampling is higher in that case with N x N px of 512 × 512 to allow a fine description of the dynamics of trapped particles. The time behavior of the action densities from the Vlasov simulation are shown in Fig. 2 together with the relevant action sum. 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. ( 49), the first exact invariant C s = S 0 + S s is well conserved, while the pump-plus-plasma wave sum C e strongly decreases (but not shown here). 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. [10] that this integral is equal to the action density S kin may be estimated by:

S kin = W kin ω e = m e c 2 ω e L L 0 dx +∞ p low (γ -1) f (x, p x , t) dp x ( 51
)
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 (see in particular the bottom panel of Fig. 2). Action behaviors in time show that the instability arises in a weak regime (in particular the pump decreases slightly during the parametric instability). Fig. 3 shows the xp x phase space representation of f afforded by our Vlasov code, at four different times during the evolution when the hot electron population is introduced. Our idea was to follow testparticle 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 e 2 /2m e γ ∂A 2 ⊥ /∂x ≪ eE x in the Vlasov equation (31), 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) (52)
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:

H w = γ ϕ eΦ lab (x -v ϕ t) + m e c 2 (γ -1) -v ϕ p x (53)
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 (54)
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 (55)
where u ϕ = γ ϕ β ϕ and β ϕ = v ϕ /c and ∆Φ = Φ 0,lab -Φ(x).

In Fig. 3 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. The distribution of particles in xp x space settles to a regime where a fraction of the initial population remains together (and forms the trapped population) and oscillates in the well located inside both separatrices.

After a few trapping and detrapping, one may assume that the distribution of particle in the domain swept by the wave-pulsating separatrices are quite uniformly distributed in the phase space in the asymptotic limit. This diffusion process is indeed observed on a long time behavior at time tω p = 3800 in Fig. 3.

To analyze in detail the physical processes, we chose now to increase the laser field amplitude to eE 0 /m e ω p c = 0.28 keeping the other parameters identical in order to increase trapping effects. The temporal behavior of the action densities and Manley-Rowe partition is shown in Fig. 4 for the electromagnetic components (on top panel) and for the plasma mode and particles (bottom panel in Fig. 4). The corresponding phase space representation of the distribution function is plotted in Figs. 5 and6. Now very strong pump depletion is observed. Fig. 4 seems to indicate that a strong trapping process is now effective. Again as expected by Eq. ( 49) the electromagnetic action sum C s is well conserved. A striking agreement can be seen between the separatrices and the phase space trapping limits for growing plasma wave. At the frame tω p = 160, the spiral structure inside the trapping limits reflects the monotonic increase in wave energy. After the plasma wave momentum is reached (for tω p = 200) the marginally trapped particles now begin to detrap.

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. Fig. 6 shows the plasma behavior after several periods of oscillations. Hamiltonian stochasticity may be obtained even for a simple plasma wave but with slowly time-periodically oscillating amplitude. Vlasov simulations show that particle dynamics is diffusive in Fig. 6 and that the phase space hole due to trapping is now filled up. This is the case on the top panel in Fig. 6 which shows the behavior of plasma in the phase space at time ω p = 900 and ω p = 920. This 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. Chaotic diffusion, occurs over the entire region of phase space where the separatrix crossing occurs. In ref.

[47] D.L. Bruhwiler and J.R. Cary have shown how separatrix-crossing theory may be applied to one-dimensional accelerating structures with slow spatial variation. In particular a two-step behavior was predicted. First the system evolves dramatically during the first few wave-particle interactions as resonant particles strongly scatter with each interaction. On the long time scale, the system slowly spreads because the separatrix-crossing region (the so-called separatrix-crossing map in ref. [47]) becomes now chaotic. Such a regime is observed here and is beyond the realm of quasilinear theory and makes the Vlasov code a good candidate to make numerical comparison with an exact N-body evolution with symplectic integrator (see for instance refs. [48,49,50] for more details).

It must be pointed out that, although the second Manley-Rowe invariant C e = S 0 + S e + S kin is well conserved; high frequency fluctuations were observed on a long time in the bottom panel in Fig. 4. This anomaly is indeed due to the occurring of another plasma mode (with a negative wave number k ′ e < 0 and k ′ e = -2 △ k < 0) which may be excited by nonlinear coupling and results here of the excitation of an anti-Stokes mode (k 0 = k ′ e + k as = -2 △ k + 6 △ k) which may be observed in the bottom panel in Fig. 6 at time tω p = 1500 and tω p = 1520.

Numerical Vlasov simulations of SRS-B in optical mixing:

We present now results on the so-called kinetic regime of Stimulated Raman back-Scattering (SRS-B) dominated by weak plasma turbulence and nonlinear shift of the EPW introduced in the beginning of Sec. 4. Vlasov-Maxwell simulations performed in the case of an open system have shown the importance of weak turbulence which is characterized by an inverse cascade process, on the nonlinear evolution of SRS-B and its saturation. The observed behavior (detailed mechanisms have been reported in ref. [46]) corresponds to the combined action of the frequency shift of the amplified EPW induced by nonlinear trapping effects, compensated for by a retuning of the wavenumber as well, so as to maintain the parametric resonance over a long time.

The simulations show that the main process that limits the continuation of SRS-B instability by Langmuir wave frequency wavevector chirping and retuning is provided by pairwise vortex merging of phase space vortices or phase space "holes".

We begin with a simulation of SRS-B in 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 .

We recall the frequency and wavenumber matching conditions between waves in the plasma for the SRS-B instability are ω 0 = ω s + ω e and k 0 = -k s + k e . Since the case of a backward stimulated Raman scattering is investigated here, the Stokes mode is a backscattered electromagnetic wave which propagates in the opposite direction of the pump wave (with a wave vector -k s and with the choice k s ≥ 0). Normalizing frequencies to the pump frequency ω 0 (because the density profile is here inhomogeneous), we have ω s = 0.666ω 0 and ω e = 0.333ω 0 and k 0 c/ω 0 ≃ 0.9578 (for the pump wave) and k s c/ω 0 ≃ 0.6016 (for the seed). The EPW has then a wave vector close to k e c/ω 0 ≃ 1.5594 which corresponds to a value of k e λ De ≃ 0.34, well into the kinetic regime. Matching conditions of SRS-B occurs right on the top of the parabolic profile at a density of n 0,max ≃ 0.0825n c , where both electromagnetic waves (pump and probe) met to start the SRS-B interaction.

The behavior of the plasma density is shown in Fig. 7 at four different times during the time evolution. The probe beam is injected at the x = L side and is propagating here in the left-hand direction. Its amplification breaks the symmetry of the interaction with respect the top of the parabolic profile. We have also checked that the beginning of the interaction is indeed in well agreement with the usual envelope model of Eqs. (44) to (46) used to describe the three-wave parametric SRS-B instability. However at time tω 0 ≃ 5600 (bottom left panel in Fig. 7) a perturbed behavior is observed resulting from strong kinetic effects. At tω 0 = 5600 there are, in the right-hand region of the top of the parabolic profile (i.e. for 2000 < xω 0 /c < 3000), strong fluctuations in density showing that plasma behavior has lost its coherence in phase space (this feature is due to the vortex coalescence in this region).

Fig. 8 shows the longitudinal electric field eE x /m e ω 0 c and its corresponding k-spectrum at two different times tω 0 = 4200 and tω 0 = 5600. We clearly distinguish the change of the regime of the instability. At time tω 0 = 4200 one observes rather coherent shape of the field and this phase corresponds to weak kinetic effects (the beginning of particle trapping). The k-spectrum exhibits a numerical value of k e c/ω 0 ≃ 1.560 in well agreement with the expected value of k e c/ω 0 ≃ 1.5594. After this phase ends, we recover the bursting behavior of SRS-B. At time tω 0 = 5600, this is the end of the fluid-like phase. However at the left side of the box, the plasma is still unperturbed keeping its coherent shape for the electron plasma field. This region exhibits very regular trapping structures (as can be seen in the top panel in Fig. 9) and contributes to the value of the SRS-B induced EPW wave vector k e c/ω 0 ≃ 1.56 in the spectrum. The k-spectrum exhibits two main characteristics at that time:

(i) First a small wavenumber shift is observed till a higher value of ke c/ω 0 ≃ 1.58 in comparison with the initial value of k e c/ω 0 ≃ 1.5594. A nonlinear shift in frequency ω e induced by trapping effects and predicted by G.J. Morales and T.M. O'Neil in ref. [51] was also observed in the simulation. Note that in ref. [46] in which simulations were performed over a long time, the corresponding down-shifted limit value is close to ωe ≃ 0.25ω 0 for the EPW rather than the original value of ω e ≃ 0.333ω 0 . While the frequency shift was seen in ref. [30], the saturation of the instability by detuning as predicted in ref. [30] is not observed in optical mixing. This is because there is the corresponding shift in wave vectors (as can be seen for instance in the Fourier k-spectrum in Fig. 8); which means that matching conditions are maintained in spite of the shift and become here ω 0 = ωs + ωe (here 1 = 0.75+0.25 in ω 0 units) and for the wave vectors k 0 = -ks + ke (here 0.9578 = -0.6078 + 1.65 in ω 0 c -1 units). The pump wave frequency ω 0 and the wave vector k 0 of course remain unmodified as a constant source of energy. Such a behavior was also observed by H.X. Vu et al in ref. [52].

(ii) We observe also in the k-spectrum of Fig. 8, at time tω 0 = 5600, the presence of a rather broad (but at lower level) spectrum around the value k V M c/ω 0 ≃ 1.14 that is here the consequence of the vortex merging phenomena in the plasma.

To aid in understanding what is happening, we present in Fig. 9 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 . 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. 9. 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. 9 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. [46]). The vortex merging of EPWs of shifted frequencies leads to the formation of lower wavenumber BGK-type structures 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 in [46] 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. Such modes behave as BGK modes and previous simulations in ref. [6] of BGK wave train using Vlasov-Poisson and previously in [53] using a Water-Bag model may also subject to such instabilities. Note such a vortex merging was also observed by Valentini et al in [43] in the case of the study of TEAWs in the nonlinear regime and using periodic and electrostatic PIC simulations.

To conclude this point we focus now on the self-sustaining feature of these structures obtained through the hole merging. In order to suppress any super-reflection in the ambipolar region, we have chosen to let the plasma evolve without laser field till a time tω 0 = 50000 and to begin the simultaneous injection of both pump plus probe at that time only. The ambipolar field is thus reduced by a factor ten during this first step of the simulation allowing to suppress any possible effect on the boundaries of the plasma. Two simulations were carried out in order to analyze the self-sustaining character of the phase space structures resulting from the LMO wave coalescence. The first simulation corresponds to a continuous injection of the pump plus probe laser beams during all the time of the simulation i.e. from tω 0 = 50000 till the end of the simulation. The result is shown on the top panel of Fig. 10. Once these phase space structures are well established (say at tω 0 = 65000), both laser beam injection are then stopped in a second step to check the BGK-character of the waves. While in the first simulation (top panel in Fig. 10) we have kept a continuous injection of pump plus probe beams during all the evolution time of the simulation, a second simulation ( the corresponding result is shown in the bottom panel in Fig. 10) is then performed turning off the pump and probe beams at tω 0 = 65000.

The top panel in Fig. 10 shows the plasma behavior (for a given cell) at time tω 0 = 75000 in the case of a simulation performed with the identical physical parameters of Fig. 9 i.e. for a pump quiver momentum of a osc = 0.025, an electron temperature of T e = 2keV and a probe seed of intensity I seed = 10 -2 I pump . Here ions are mobile with a mass of m i = 1836m e and an ion temperature of T i = 0.25T e . One begins to see accelerated electrons which were reflected in the ambipolar sheath field and are passing the EPW region but now moving leftward towards the laser. The bottom panel corresponds now to the case where the pump wave is turned off at time tω 0 = 65000 and we observe that the LMO waves, continuously generated by the SRS-B instability have now disappeared, while on the other hand the BGK-like waves are self-sustained and arrive to persist over a long time after the pump was turned off. These simulations show clearly that phase space hole merging leads to the occurring of BGK-like (but here non stationary) structures.

LPI in the relativistic regime

In this section we focus now on the laser-plasma interaction in the relativistic regime. Recent developments in laser technology have made possible the generation of short laser pulses at intensities above the limit of 1.368 × 10 18 W cm -2 (the corresponding quiver momentum being a osc = 1 for a linear polarization of the pump electromagnetic wave) and therefore have opened possibilities to explore laser-plasma interaction in new regimes. We have seen that electromagnetic waves are parametrically unstable when they propagate through a plasma. At lower intensities (a osc ≪ 1), the instabilities are clearly identified as SRS for purely electronic parametric instabilities. When electrons quiver at relativistic velocities, plasma is now subject to new instabilities not evident at lower intensities. Not only relativistic effects modify strongly the features of SRS but lead to the relativistic modulational instability (RMI) and the relativistic filamentation instability (RFI). There are two aspects of the same physical process, the wave vector being respectively, parallel or perpendicular to the pump wave vector direction. Several remarks must be here pointed out:

(i) First in this strongly relativistic regime of laser-plasma interaction, an unadapted use of the time-splitting scheme of Knorr and Cheng (see ref. [20] for more details) may lead to a bad density conservation and triggers to a numerical instability when canonical invariants are used in the transverse directions. In order to avoid these difficulties, a direct 2D advection solver must be applied in this regime.

(ii) The semi-lagrangian Vlasov code was used to investigate the SRS parametric instability in the relativistic case in ref. [54] in the 1D case then in ref. [55] in a 2D plasma. The most striking point is that we did not introduced any initial density perturbation in our simulation. The instability starts just from the round-off errors, which for a 64-bit computer is of order of 10 -15 . Such a result clearly demonstrates that the low noise character of the Vlasov code allows an extremely accurate study of the growth of such SRS-like instabilities over a large number of decades.

(iii) High longitudinal temperature of the bulk of the plasma modifies the SRS instability even in the relativistic regime. Some theoretical and numerical studies have shown the suppression of SRS in the relativistic regime. SRS shows new features in the relativistic regime since several branches of instability (SRS-F, SRS-B and RMI -the relativistic modulational instability -may merge and lead to a same instability). A water bag model was used in ref. [56] by Z.M. Sheng et al showing that high (relativistic) temperatures, in the longitudinal direction or the laser propagation direction, significantly reduce the growthrate of the SRS instability and that the maximum growth tends to shift towards small k values. In refs. [57] the authors, using PIC simulations, have shown that for sufficiently a hot and dense plasma, the SRS instability becomes non-resonant, turning into Compton-like scattering rather than Raman-type.

(iv) Periodic systems, while excellent for untangling the basic physics, are quite unrealistic. Boundary conditions and in particular the effective laser propagation are important as in the self-induced transparency (SIT), i.e. the laser propagation through an overdense plasma in which the propagation is classically forbidden. Note that for such an overdense plasma the plasma density n 0 is then bigger than the critical density n c = n 0 ω 2 0 /ω 2 p . SIT differs radically from the conventional hole boring scenario or forward motion of the critical surface caused by the action of the ponderomotive pressure. SIT is characterized by an anomalous longitudinal electron heating. We may again speak about a new kinetic regime for which particle trapping seems to play again a major role as in the case of SRS for underdense plasmas. The situation is however somewhat different in the case of overdense plasmas. The interaction of such laser beams with moderate underdense to slightly overdense plasma is of particular interest in the "fast ignitor fusion concept" proposed by Tabak in ref. [58] and is expected to display new physical phenomena.

(v) The last comment concerns the formation of the KEEN waves. The situation is however somewhat more complex in the regime where the driving lasers are at relativistic intensities. In particular the plasma is now subject to strong electron parametric instabilities as SRS which might mask the formation and growth of these KEEN waves. Indeed the strong heating which accompanies these SRS instabilities (and which is responsible also of their suppression on a long time behavior) allows the growth of these modes. KEEN waves may be seen as a transition process from the SRS instability to a Compton-like instability in regime dominated by hot (longitudinal) electron temperature.

Conventional hole boring scenario

Before investigating new kinetic effects associated with the anomalous absorption met for moderate overdense plasma, it is interesting to come back to the standard model or hole boring model. The mechanism involves the ponderomotive force associated with the high intensity laser pulse, which in the 1D picture, considered here, pushes forward particles. The laser pulse, acting as a piston, drives a collisionless ion shock wave into the plasma slab which reflects ions thus accelerated up to velocities twice the wave front velocity (see for instance the top panel of Fig. 11). This mechanism was first studied by Denavit in [59] using PIC simulations. Of course such a mechanism is accompanied by the occurring of ion kinetic effects, however a weak electron heating is usually observed in this regime. Recently the case of a circularly polarized wave was investigated in ref. [60] showing that one obtains moderate ion energies and weak electron heating.

We focus now on numerical results given by our Vlasov code. Here only circularly polarized laser pulses were considered to simplify the presentation of the physical processes. Since an homogeneous plasma slab surrounded by vacuum regions on both sides of the slab is considered, we use again a standard normalization. The time is given in inverse of the plasma frequency ω -1 p , the lengths in cω -1 p unit and the electron momenta in m e c (momenta for ions being given in m i c).

Fig. 11 shows the presence of the expected electrostatic collisionless ion shock in the ion phase space in the case of the standard hole boring model. The corresponding electron phase space is shown in the bottom panel at the same time tω p = 900. The simulation was performed with the physical parameters n 0 = 4n c and a osc = 5. The electron temperature is chosen to T e = 100keV . Ions are mobile with a temperature of T i = 20keV and we have chosen a ratio of the ion (rest) mass to the electron mass to m i /m e = 1836.

Two remarks must be pointed out at this step. (i) First the wave front velocity v F was found in well agreement with the theoretical value predicted by Denavit in [59] i.e.

v F c = n c m e n 0 m i a 2 osc ( 56
)
Here the numerical value is estimated to v F ≃ 0.054c in well agreement with the value obtained by Eq. ( 56) leading to an estimation of v F ≃ 0.058c.

(ii) The observed ion shock velocity in the phase space in Fig. 11 (top panel) is close to v shock ≃ 2v F ≃ 0.11c (which corresponds to a value of the Mach number M close to M = c s /2v F ∼ 4 -5 where the ion acoustic velocity c s = k B T e /m i ≃ 0.025c). Such a behavior has been already predicted in [60] by using PIC simulations. This mechanism action on ion dynamics is accompagnied by an electron heating located in the region where the ion shock takes place. Thus the penetration of the laser pulse is associated with the ponderomotive pressure leading to a quasiinstantaneous electron heating and the formation of high Mach number ion shock at twice of the velocity of the wave front. As we will see after, SIT is a different physical process since it leads to a stronger electron heating till relativistic temperature of several MeV which is not reached here.

Enhanced electron kinetic regime met in SIT

It is generally recognized that PIC codes present a high level of numerical noise and can overestimate the plasma heating due to their numerical heating. The idea to use a Vlasov code to describe SIT for high intensity electromagnetic laser pulse becomes then natural (for electrons the gain is immediate but for ions the advantage is somewhat less evident).

From our Vlasov simulation, a new heating mechanism has emerged. The penetration by SIT of a pump wave is characterized by a beatwave mechanism between the incoming pump wave (ω 0 , k 0 ) and its reflected electromagnetic wave (ω r , k r ) at the wave front. The reflected wave frequency is lowered by a relativistic Doppler shift which arises at the wave front, which acts here as a moving mirror. We can indeed obtain a good estimation of the characteristics of this reflected electromagnetic wave in terms of frequency and wavevector by considering the relativistic Doppler shift induced by this moving wave front. An estimation of this frequency ω r and wave vector k r can be found in ref. [64] leading to the following expressions:

ω r = ω 0 γ 2 F 1 + β 2 F -2 β F β 0 (57)
k r = k 0 γ 2 F 1 + β 2 F -2β F β 0 (58)
where

β 0 = ω 0 /k 0 c . Here γ F = 1 -β 2 F -1/2
is the Lorentz associated with

β F = v F /c
. It is the beating between these two electromagnetic waves which now determines the characteristics in frequency and wave vector of the plasma response.

The beatwave drives a low-frequency electron acoustic wave (EAW) (ω eaw , k eaw ).

Matching conditions required for such a resonance may be written in the following form:

ω 0 = ω r + ω eaw and k 0 = -k r + k eaw (59)
The electromagnetic waves (pump plus reflected wave) satisfy the standard (linear) dispersion relation for circularly polarized waves in a hot plasma:

ω 2 0,r = ω 2 p γ -1 + k 2 0,r c 2 (60)
The low frequency EAW which is driven by the beating are similar to the acoustic modes introduced by S.P. Gary in ref. [61] and may occur in a plasma that comprises two electron components distinguished by widely disparate temperatures.

Here the EAWs are driven in a nonlinear way and are characterized by the formation of trapping structures (here "bumps" in phase space rather than the usual "holes" met in TEAWs as can be seen in the middle panel of Fig. 14). Notice also that in the case of TEAWs or KEEN waves no two electron components and well-separated temperatures need to be invoked. Furthermore these EAWs, even in the nonlinear step, are no self-sustaining property i.e. cannot survive after the pump wave was turned off.

The "forced" electron mode (in the sense that it disappears once the laser is turned off) obeys to the condition:

ω eaw = k eaw v F = (k 0 + k r ) v F ∼ k eam v th n cold n hot (61)
where we assume that an electron component is hot (of density n hot ) and an another cool (with the subscript cold in density). Eq. ( 61) can be easily obtained via a Lorentz transform in frequency. In the frame of the wave front, the frequency of this beat mode is then given by ω ′ eaw = γ F (ω eaw -k eaw v F ) = 0 since in this frame both pump and reflected light have the same frequency.

In this regime, the expansion of ions into vacuum is too slow for the electron motion to destabilize the electron modes. This process was observed in ref. [62] for short systems and fast electron scales. We have referred to this unusual beatwave process involving (nonlinear) EAWs as the Doppler BeatWave Heating (DBWH) in [62], to emphasize both requirements, the Doppler shift and the anomalous electron longitudinal heating. This is a novel version of the beat-wave heating scenario first introduced by Cohen et al in ref. [63] (involving a Langmuir wave in the Cohen's scenario) and which can be applied in an overdense plasma here. The underlying physics of the observed anomalous electron heating is essentially the generation of these EAWs in a nonlinear regime. To be an efficient heating process, the plasma mode must not be suppressed by the high Landau damping of any SRS-like processes in high temperature. Here the plasma response is well adapted since EAWs are weakly affected by Landau's damping and have weak phase velocity located in the bulk of the distribution function. This last point allows to "trap" particles in high quantity. It must be pointed out that the frequency range in which the beat mode can be excited is not limited to (nonlinear) EAWs. Ion Acoustic Waves (IAWs) may also be excited. This excitation however arises when the corresponding ion acoustic velocity c s is very close to the wave front velocity v F . Fig. 12 displays the electromagnetic and electrostatic spectra in frequencies for the following physical parameters n 0 /n c = 2 and a quiver momentum of a osc = √ 5 corresponding to an irradiation of Iλ 2 0 = 1.37 × 10 19 W cm -2 µm 2 . We have chosen to increase the electron temperature to T e = 400keV to minimize SRS-like relativistic instabilities. For these physical parameters, the corresponding value of the ratio of relative electron density to the critical density is then n 0 γ -1 /n c ≃ 0.32, which indicates that the plasma now becomes transparent for the electromagnetic pump wave. The upper panel in Fig. 12 corresponds respectively to the reflected light spectrum. The ω-spectrum of the electron density is shown in the bottom panel. Both spectra are calculated in the time domain of 0, 600ω -1 p . The top left panel exhibits a dominant peak at frequency ω r ≃ 0.51ω p in well agreement with the theoretical value predicted by Eq. ( 57) giving ω r = 0.526ω p . Using Eq. ( 58) the corresponding value of the wave vector is then k r = 0.511k 0 ≃ 0.29ω p /c. Using n 0 /n c = 2 we obtain a pump wave frequency of ω 0 = 1/ n 0 /n c = 0.707ω p and k 0 c/ω p = 0.58. Using γ -1 ≃ 0.16 the (linear) dispersion relation (60) leads to a value of the pump frequency of ω 0 = 0.704ω p in well agreement with the numerical value obtained in spectra. The wave front velocity used in Eqs. ( 57) and ( 58) (and measured directly in phase space plots in Figs 14) is here close of v F = 0.22c. In ref. [62] this wave front velocity was found to depend of γ -1 i.e. the parameters a osc and T e and finally of the ratio n 0 /n c . Thus for a given value of γ -1 the wave front velocity v F decreases in a linear way as a function of n 0 /n c untill the plasma layer becomes opaque. For moderate overdense plasmas v F can be found bigger than the value found by Eq. ( 56) met in the hole boring scenario.

The DBWH arises when matching conditions (59) are satisfied. Here these conditions lead to ω 0 = ω r +ω eaw (i.e. 0.707=0.51+ 0.197 in ω p units) and k 0 = -k r +k eaw (0.58=-0.30 + 0.88 in ω p /c units). The electrostatic spectrum (the bottom panel in Fig. 12) illustrates the occuring of the EAW into the plasma at frequency of ω eaw ≃ 0.19ω p but another peak located around a higher value of ω K ≃ 0.59ω p , but nevertheless weaker than the usual plasma frequency ω p . This second peak may be attributed to the excitation of a KEEN wave by nonlinear coupling at the same wavenumber k eaw .

The beating induced by the Doppler shift of the pump wave at the moving wave front leads to a strong electron heating. Fig. 13 displays the mean distribution function versus p x /m e c at three different times during the interaction. The top panel shows the plasma state in the beginning phase in time at tω p ≃ 92 when the laser pulse has not reached the plasma slab (the plasma is then very close to the initial plasma state). The interaction of the laser pulse with the plasma layer is indeed characterized by a two-step process. The middle panel at time tω p ≃ 384 corresponds to the DBWH process: we see clearly a two-temperature plasma with the existence of a cold return current. While the laser pulse has crossed all the plasma layer at time tω p ≃ 1538 (bottom panel in Fig. 13), the plasma heating continues which seems indicate that another new process is now present. Such two-step plasma behavior was already observed in PIC simulations in ref. [64] for moderate overdense plasmas.

Figs. 14 show the behavior of the plasma, at three different times during the first phase of the interaction i.e. when DBWH takes place. Fig. 16 shows what happens on a longer time in the second phase of the instability. As the pump wave penetrates by SIT inside the plasma, new phase space vortices are formed and their creation was associated with the displacement of the wave front. The wave front is now moving at the velocity v F ≃ 0.22c. Till the time tω p ≃ 480 in Fig. 16, which corresponds to the fact that the electromagnetic pulse has reached the other side of the plasma, the KEEN wave excitation does not arise and we just observe in Fig. 15 the formation of the EAW characterized here by "bumps" in phase space and by a two-temperature electron population. We observe that the cold population is not located at the zero value in momenta indicating that a strong return current is present inside the plasma. At tω p ≃ 481 in Fig. 16, the plasma exhibits the formation of seven vortices for a plasma length of 72.5c/ω p or 6.4λ 0 , which is close to the expected value of k eaw c/ω p = 0.88 (the value measured in the electron phase space being close to k eaw c/ω p = 0.86). A simulation with mobile ions with a mass of m i = 1836m e was also performed for the same physical parameters. As expected, ion dynamics plays a minor role and no ion collisionless shock has been observed in this regime of the laser-plasma interaction, showing that the DBWH strongly differs from the conventional hole boring mechanism.

Second phase of SIT: Stimulated KEEN Slab Scattering

Once the plasma layer crossing by the wave front is achieved, the beatwave normally ceases. However the behavior of the electron distribution function in phase space exhibits, even in this limit, trapping structures but at a different wave number (see Figs. 16). Indeed after the pump wave has crossed all the plasma layer, the plasma heating does not stop. Consequently another physical mechanism is now effective when the plasma is almost transparent to the pump light. The simulation plasma is then in the form of a relatively sharp-edged slab and this may support (more and less) an electromagnetic cavity-like radiating mode structure whose wavevector k s is determined by the slab width L slab and whose frequency is close to the (relativistic) plasma frequency ω s ≃ ω p γ -1 1/2 . Indeed in that regime the change in physical parameters in term of wave number and frequency of the pump do not modify the value obtained for ω s which was always observed at the local plasma frequency both for the backward and forward direction of propagation of the scattered mode.

The second point is that now the plasma temperature is high enough to allow the growth of KEEN waves in the relativistic regime. The top panel of Fig. 15 shows the forward electromagnetic spectrum and its corresponding electrostatic spectrum is plotted in the bottom panel. The spectra have been obtained in the time interval 600, 1500ω -1 p . As expected the electromagnetic spectrum exhibits a central peak at the pump frequency ω 0 = 0.707ω p and a second peak at a frequency close to the local plasma frequency ω s ≃ ω p γ -1 1/2 ≃ 0.40ω p .

In systematic simulation studies of this scheme (we have varied the ratio electron density to the critical density) our investigation reveals that this electromag-netic scattered wave was driven near critical, i.e. with a frequency always given by ω p γ -1 1/2 and with a wavevector k s close to zero (in fact k s was found to be a multiple of π/L slab and k s = nπ/L slab ≪ 1). This feature is the signature of the Stimulated Slab KEEN Scattering (SSKS) instability. In ref. [65] we have shown that, for somewhat underdense plasma (n 0 /n c 0.25), the simulation resonance was in fact determined by the beating of the pump with a new radiating pseudo-cavity electromagnetic mode for the slab at a frequency close to ω p with relatively low loss. The resonant conditions of this SRS-like instability involving a (ω K , k K ) KEEN wave are then:

ω 0 = ω s + ω K and k 0 = ±k s + k K ∼ k K (62)
Here we have for frequencies 0.707 = 0.40 + 0.307 in ω p units (i.e. ω s = 0.40ω p and ω K ≃ 0.30ω p i.e. well below the value of the local plasma frequency) and for the wave vectors

k 0 = k s + k K ≃ k K with k 0 c/ω p = 0.58, k s c/ω p ∼ 0.1 and k K c/ω p = 0.48.
The low-frequency KEEN wave can be seen in the bottom panel in Fig. 15: a peak located at ω K ≃ 0.30ω p is clearly visible (and with its two first harmonic ω 2K and ω 3K ).

As tω p = 865 (in the middle panel in Fig. 16) as SKSS is well established in simulation, new trapping structures are formed in phase space. At that time we see clearly the formation, for positive momentum, of six or seven vortices for a plasma slab of 74c/ω p length (i.e. somewhat 7λ 0 with λ 0 being the pump wavelength into the plasma). This corresponds to a wavelength for the KEEN wave of λ K ω p /c ≃ 12.25 or equivalently to a wavevector of k K c/ω p ≃ 0.50 in well agreement with the expected value. For negative values of the momentum, one observes the particle recirculation from the right-hand plasma edge, with clear phase mixing convection as the electrons return against the incoming pump wave. Trapping structures initially generated for positive values of the momentum (since the corresponding wavevector of this wave k K is positive) are now destroyed for negative velocities (since the parametric SKSS instability is ineffective for negative values of k K ).

Plasma Kinetic Waves

The lack of low-level narrow-band resonance in the electrostatic oscillations precludes the usual stimulated scattering scenario as met for instance in SRS-B in sec. 0.4 i.e. a three-wave resonance amplification from a low-level seed. In the usual analysis the initial frequency of the scattered electromagnetic wave and the driven electrostatic response are necessarily chosen to satisfy three-wave resonance. When this electrostatic resonance is very broad, the resonance mechanism of this electrostatic wave amplification is then lost. In fact it is an electromagnetic daughter wave that sets the resonance via a pseudo-cavity mode for instance in the SKSS instability. It is this radiating pseudo-cavity plasma slab structure that determines the frequency ω s (just above the local plasma frequency) and the wave number k s ∼ π/L slab , which can be combined with the pump wave to produce a well-defined difference ponderomotive force drive which can thus drive up the (non resonant) KEEN waves. The situation resembles the case of highly Landau damped SRS which is in effect stimulated Compton-like scattering (SCS) acting directly on electrons. Thus the formation of these KEEN waves appears as an intermediate regime between SRS and SCS provided that the plasma temperature is high enough. Contrary to SCS, the plasma keeps a collective coherent response in the form of these KEEN waves. However one of the most fundamental aspect is their kinetic nature.

We have performed a numerical experiment in the case of a periodic system: a beatwave scenario between both electromagnetic pump (ω 0 , k 0 ) and a pseudocavity (whose features are in terms of frequency and wavevector ω s ≃ ω p γ -1 1/2 and k s ∼ 2π/L slab = △k ; △k being the fundamental wave number). The beating between both electromagnetic modes is used to excite a KEEN wave into the plasma. Since the simulation box is periodic, we assume perfect k-matching: the pump wave vector k 0 is chosen as k 0 = 12 △ k. The scattered mode k s = -△ k is chosen as small as possible (as in the SKSS scenario for an open system). The KEEN wavevector is then k K = 13 △ k. Different initial conditions were tested: the choice of a twotemperature Maxwellian (i.e. with a cold electron component at the temperature T cold = 15keV and a second hot electron population at T hot = 150keV ) allows a better stabilization of the KEEN wave but are not necessary to excite the KEEN wave.

The physical parameters may be found in ref. [66]. We use again action transfer diagnostics as used in sec. 0.4.1. However when KEEN waves are involved, the longitudinal field is strongly Landau damped and it is now necessary to take into account the potential vector in the calculation of the separatrices. In fact the ∆Φ = Φ 0 -Φ lab (x) term used in eq. ( 55), and which allows the determination of both lower and upper separatrix limits in phase space, is now replaced by the following expression:

∆Φ = Φ 0 -Φ lab (x) + m e c 2 2e eA ⊥ m e c 2 (63)
and where Φ 0 is the generalized potential at an X-point where the two separatrices join. Fig. 17 shows the electromagnetic action densities versus time i.e. on the top panel, the pump action density S 0 (in solid line), the action of the corresponding scattered mode S s (in dashed line) and their sum C s = S 0 + S s (in dotted line ). The bottom panel shows the time evolution of the pump S 0 (again in solid line) but now with the plasma mode action S e (here negligible) and the kinetic part S kin evaluated from Eq. ( 51) and the sum C e = S 0 + S e + S kin here close to the value S 0 + S kin since S e remains at a very small level. Fig. 17 shows the action S kin transferred to trapped particles is now higher in comparison to the SRS-F case investigated in sec. 0.4.1. Note also since the electric action S e obtained from the data of the longitudinal electric field is here negligible, there are in fact only three actors in the system and that their sum C e = S 0 + S kin is well conserved in the Vlasov simulation.

The corresponding phase space of the electron distribution function is plotted in Fig. 18 showing the formation of electron holes in phase space. From these results, it is become clear that a simple three-wave oscillator model based on a fluid description for the envelope equations is not adequate for studying such a kinetic instability and that particle trapping effects must be correctly included. Thus the nature of these KEEN waves is purely kinetic and require, like the BGK waves, a kinetic and nonlinear treatment.

Conclusion and outlook.

The numerical simulation is a huge topic with hundreds of papers published every year. Eulerian 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 LPI even for 2D problems. The very good resolution in phase space allows a precise study of wave-particle interaction including particle trapping and particle acceleration. We are investigating one such curious electron phase space anomaly, namely the KEEN waves, a version of the classic BGK concept but unrelated to usual electron plasma waves. 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 LPI, or Water-Bag model or even adiabatic invariants 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 trapped ion modes in ref. [67] 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? (i) 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 gyroaveraged 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 driftkinetic Vlasov equation (see refs. [6] and [68]) 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. [24] and [25] for more details) which introduces a special class of initial conditions allowing to reduce the full kinetic Vlasov equation into a set of hydrodynamical equations.

(ii) 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 modeling wave-particle interactions in LPI. 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 interdisciplinary research that targets the development of capabilities that "bridge" various areas of plasma physics together with computer science and applied mathematics. Figure 3: For a periodic SRS-F case, phase space representation of the electron distribution function corresponding to the situation shown in Fig. 2 when a small hot electron population is introduced initially. SRS-F enters now in a weak "kinetic" regime characterized by a small trapped electron population localized inside the wave-pulsating separatrices calculated using the data of the electrostatic potential given by the code.

Figure 4: From a spatially periodic semi-lagrangian Vlasov simulation of SRS-F: time evolution of the pump action S 0 (in solid line), the corresponding scattered wave action S s (in dashed line) and their action sum C s = S 0 +S s (top panel). In the bottom panel: pump action density S 0 plotted again in solid line, but now with the electron plasma wave action S e , the transfer action S kin to trapped particles, together with their sum C e = S 0 + S e + S kin where S kin accounts for particle accelerated above the lower separatrix (see Figs. 5 and6). A stronger wave particle-interaction is now observed by increasing the pump quiver momentum to a osc = 0.107; the other parameters being the same than those used in Fig. 2. 6: Continuing on from the previous set of phase space frames in Fig. 5, we show now the plasma behavior after several periods of oscillations. Hamiltonian stochasticity may be obtained even for an unique EPW mode with slowly timeperiodically oscillating amplitude. In the top panel at tω p = 900 and tω p = 920 the phase space hole due to trapping is now filled up by electrons. On a long time scale at tω p = 1500 the system slowly spreads in the separatrix-crossing region. The parameters are here a quiver momentum of a osc = 0.107 and an electron temperature of T e = 15keV with a small (5%) population of hot electrons (at 100keV ) introduced to enhance the wave-particle interaction. Figure 10: Phase space representation of the electron distribution function in the case of a parabolic profile using identical parameters than those used in Fig. 9 but now with mobile ions. We have thus a osc = 0.025 for the quiver momentum, T e = 2keV for the electron temperature and T i = 0.25T e for the ion temperature and m i = 1836m e for the ion mass. We have chosen to let the plasma evolve without laser field till a time tω 0 = 5 × 10 4 and to begin the simultaneous pump and seed injection at that time. Top panel: the pump and probe laser beams were continuously injected: we see the formation of LMO phase space vortices for positive velocity and accelerated electron which were reflected in the ambipolar sheath field and are now passing the EPW region but now moving leftward towards the laser. Such holes results from previous pairwise vortex merging. Bottom panel: the pump and probe laser beams were turned off at tω 0 = 6.5 × 10 4 and we observe that LMO waves have now disappeared while the BGK-like structures with negative phase space structures persist in time.

Figure 11: Vlasov semi-lagrangian simulation performed for hole boring scenario for an uniform short plasma slab (with an electron density of n 0 = 4n c and a quiver momentum of a osc = 5). The ion behavior in phase space is plotted in the top panel. The laser pulse drives a collisionless ion shock wave into the plasma slab which reflects thus accelerated up to velocities twice the wave front velocity. Such a mechanism is accompanied by an heating in the electron distribution function, plotted here at the same time in the bottom panel. This electron heating however remains weaker than that observed in SIT (see for instance Fig. 14). bottom) we have used the separatrix limit to determine S kin (by computing the relativistic kinetic energy above the lower separatrix). For the whole phase space representation f has 13 vortices. We have just made a zoom and have not represented all the phase space. We have also superimposed the separatrices on the curves.

Figure 1 :

 1 Figure 1: For a periodic SRS-F case, time evolution of the different action densities: electromagnetic pump S 0 (in solid line), scattered Stokes mode S s (in dashed line) and their sum C s on the top panel. In the bottom panel we have plotted the EPW action density S e (in dashed line) superimposed on the pump S 0 (in solid line). Without adding a hot electron population we observe a fluid-like behavior with wellconservation of Manley-Rowe invariant C s (top panel) and C e (bottom panel). The physical parameters are T e = 15keV and a osc = 0.0538.

Figure 2 :

 2 Figure 2: 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. The well-conservation of the pump plus Stokes action density is shown on the top panel while it is now necessary to take into account the transfered action to trapped particles S kin to obtain a good conservation of the second invariant C e = S 0 + S e + S kin of the Manley-Rowe partition. Notice the weak depletion of the pump now. Physical parameters are identical to those used in Fig.1.

Figure 5 :

 5 Figure 5: Phase space (p x /m e c vs xω p /c) electron distribution function showing the evolution of the trapped electron generation in the case of SRS-F and corresponding to action evolution of Fig. 4. The separatrices obtained by Eqs. (55) are also drawn. We observe the trapping spiral structures associated with a growing plasma wave during the first pump depletion.

Figure

  Figure6: Continuing on from the previous set of phase space frames in Fig.5, we show now the plasma behavior after several periods of oscillations. Hamiltonian stochasticity may be obtained even for an unique EPW mode with slowly timeperiodically oscillating amplitude. In the top panel at tω p = 900 and tω p = 920 the phase space hole due to trapping is now filled up by electrons. On a long time scale at tω p = 1500 the system slowly spreads in the separatrix-crossing region. The parameters are here a quiver momentum of a osc = 0.107 and an electron temperature of T e = 15keV with a small (5%) population of hot electrons (at 100keV ) introduced to enhance the wave-particle interaction.

Figure 7 :

 7 Figure 7: For the parabolic profile (maximum density of n 0 ≃ 0.0825n c , a osc = 0.025and T e = 2keV ) the typical density behavior at four different times. At tω p = 3500 the probe amplification breaks the symmetry of the interaction with respect to the top of the parabolic profile. At tω p = 6300 the lose of coherence of the shape for the electron density as the result of the pairwise vortex merging in that region.

Figure 8 :

 8 Figure 8: For SRS-B in optical mixing and for the parabolic plasma simulation of Fig. 7, we show here the longitudinal electric field vs xω 0 /c at two different times, together with their k-spectra. Notice the increase in the wave number k e at tω 0 = 5600 (bottom right panel) from the initial k e value to a larger value of k e . This is due (see text) to a Langmuir Morales-O'Neil (LMO) shift with intensity.

Figure 9 :

 9 Figure 9: 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 of Fig. 8 (bottom right panel).

Figure 12 :

 12 Figure 12: For the study of self-induced transparency scenario, display of the reflected electromagnetic (top panel) and of the electron density (bottom panel) spectra in frequency illustrating the DBWH mechanism. The top panel shows that the dominant part of the backscattered light is induced by the Doppler shift due to the moving wave front. The bottom panel exhibits the beatwave frequency close here to ω eaw = 0.20ω p (here an electron acoustic wave) and a second mode at ω K ≃ 0.59ω p corresponding to a KEEN wave. The physical parameters are n 0 = 2n c , a osc = √ 5 and T e = 400keV ; ions being fixed.

Figure 13 :

 13 Figure 13: For a hot plasma of electron temperature of T e = 400keV and fixed ions, the mean distribution function versus p x /m e c at three different times. The top panel show the initial plasma state when the laser has not reached the plasma slab. At time tω p ≃ 384 we see the increase of the electron temperature induced by the Doppler beatwave heating. After this phase the plasma heating does not cease once the plasma layer crossing by the wave front is achieved (at time tω p ≃ 1538). A stimulated slab KEEN scattering instability takes place (bottom panel) leading to a strong plasma heating.

Figure 14 : 5 .

 145 Figure 14: For SIT scenario and an uniform and hot plasma with an electron temperature of T e = 400keV , we show the representation of the electron distribution function in phase space for a numerical solution carried out with n 0 = 2n c and a osc = √ 5. The panels show the vortices corresponding to EAW modes induced by the Doppler shift of the incoming pump wave after reflection. Notice the feature of a two-temperature component of the distribution function, feature of the excitation of an EAW.

Figure 15 :

 15 Figure 15: For the self-induced transparency (SIT) scenario of an uniform plasma slab of n 0 = 2n c , display of the (forward) electromagnetic ω-spectrum (on top panel) and of the electron density (on bottom panel) calculated in the time interval 600, 1500ω -1 pfor the same simulation presented in Fig.12. We focus here on the second phase of the instability dominated by SKSS i.e. when DBWH is ineffective. The low-frequency KEEN wave can be seen in the bottom panel with a peak located around ω K ≃ 0.30ω p well below the local plasma frequency.

Figure 16 :

 16 Figure16: Continuing on from set of phase space frame in Fig.14and for the selfinduced transparency (SIT) scenario of an uniform plasma slab of n 0 = 2n c , we show here the electron distribution function in phase space. The top panel corresponds to the formation of an EAW induced by DBWH. In both middle and bottom frames, the plasma slab crossing by the wave front is now achieved and new trapping structures are now observed. KEEN vortices initially generated for positive values of the momentum were destroyed after reflection on the ambipolat sheath field (SKSS being not possible for a negative value of k K ). The corresponding spectra in frequency are shown in Fig.15

Figure 17 :

 17 Figure 17: For a periodic plasma and a beatwave scenario between the pump wave and an initially imposed backscattered electromagnetic mode (playing the role of the pseudo-cavity in SKSS), we plot the action densities versus time: the electromagnetic action density S 0 (in solid line), the backscattered action density S s (in dashed line) and their sum C s = S 0 + S s on the top panel. The second Manley -Rowe invariant C e = S 0 + S e + S kin is plotted on the bottom panel. Notice that the plasma mode S e (here corresponding to the action of the longitudinal electric field associated to the KEEN wave) is negligible and the mot important contribution comes from S kin which seems to be a feature of the KEEN wave growth.

Figure 18 :

 18 Figure 18: Phase space representation of the electron distribution function for the same simulation presented in Fig. 17. For the calculation of the action transfer to trapped particles S kin used in the Manley-Rowe partition of Fig. 17(bottom) we have used the separatrix limit to determine S kin (by computing the relativistic kinetic energy above the lower separatrix). For the whole phase space representation f has 13 vortices. We have just made a zoom and have not represented all the phase space. We have also superimposed the separatrices on the curves.

  Figure 18: Phase space representation of the electron distribution function for the same simulation presented in Fig. 17. For the calculation of the action transfer to trapped particles S kin used in the Manley-Rowe partition of Fig. 17(bottom) we have used the separatrix limit to determine S kin (by computing the relativistic kinetic energy above the lower separatrix). For the whole phase space representation f has 13 vortices. We have just made a zoom and have not represented all the phase space. We have also superimposed the separatrices on the curves.

Acknowledgments:

Contributions from Drs P. Bertrand, T.W. Johnston, E. Sonnendrucker via the CALcul and VIsualization INRIA Project are warmly acknowledged. The work is supported by the ANR project NT05-1-41542 VALSOV and The the IDRIS center (Institut du Développement et des Ressources pour l'Informatique Scientifique, CNRS, Orsay, France) for computer time allocation on their supercomputers.