**2. The relevant equations of the Eulerian Vlasov code and the numerical scheme**

We study this current problem by using an Eulerian Vlasov code for the numerical solution of the one-dimensional (1D) relativistic Vlasov-Maxwell equations. The relevant equations for the Eulerian Vlasov formulation are those previously presented in references [12-15] for instance. We present here these equations in order to fix the notation. Time *t* is normalized to the inverse plasma frequency *ω<sup>p</sup>* −1 , length is normalized to *l* <sup>0</sup> =*cω<sup>p</sup>* −1 , velocity and momentum are normalized respectively to the velocity of light *c* and to *Mec*, where *Me* is the electron mass and *c* is the velocity of light. We have the following Vlasov equations for the electrons and the ions distribution functions *f <sup>e</sup>*,*<sup>i</sup>* (*x*, *pxe*,*<sup>i</sup>* , *t*):

$$\frac{\partial f\_{e,i}}{\partial t} + \frac{p\_{xe,i}}{\mu\_{e,i}\mathcal{V}\_{e,i}} \frac{\partial f\_{e,i}}{\partial \mathbf{x}} + (\mp E\_x - \frac{1}{2\mu\_{e,i}\mathcal{V}\_{e,i}} \frac{\partial a\_{\perp}^2}{\partial \mathbf{x}}) \frac{\partial f\_{e,i}}{\partial p\_{xe,i}} = 0 \tag{1}$$

where *γe*,*<sup>i</sup>* =(1 + (*pxe*,*<sup>i</sup>* / *μe*,*<sup>i</sup>* )<sup>2</sup> + (*a*<sup>⊥</sup> / *μe*,*<sup>i</sup>* )2)1/2 .

box, SRBS finally grows after KEEN waves reach that area and change the local distribution function by softening it. SRBS eventually swallows up the KEEN wave and dominates since its growth rate is far higher. This new scenario confirms that nonlinear trapping evolution of SRBS is not just a question of EPWs but also of KEEN waves and SKEENS and that SRFS wavenumber perturbations can initiate the latter if a standing wave already exists much before SRBS can occur. The later evolution of all these processes is very complicated still involving positive slope electron distribution functions which will then accelerate the self-destruction of these modes and render the picture even more transient, intermittent and chaotic. We stop the simulations short of that eventuality where even Brillouin scatter begins to occur and dynamics and predictions become more challenging to track requiring ion acoustic waves and fluid

Nonlocal and collective kinetic effects are involved in this physics. A direct Vlasov solver, capable of resolving these kinetic processes, is used here to address some aspects of the scattering properties of SRS and SKEENS. The code evolves relativistically both electrons and ions. In Vlasov codes used to simulate these problems, noise and other numerical fluctuations are very low for the SRBS to grow from, therefore it is usual in several simulations to stimulate artificially the counter-propagating daughter light wave at a low level as an injected seed, in order to enhance the SRBS growth and to allow the saturation phase to be reached rapidly. This saturation results from the competition of non-linear effects which include frequency shift, pump depletion, damping reduction, trapped particle instability, spatiotemporal chaos, among others. A detailed study of the resulting distribution function obtained at saturation has been presented, for instance, in Strozzi *et al*, 2007, which showed that the stage following saturation involved the transformation of Raman Langmuir waves into a set of beam acoustic modes or BAM (see also [24,25]), and an EAW appears at this stage with a weak reflected light that phase-matches for scattering off this mode, a process called electron acoustic scatter (EAS). SEAS has been experimentally observed in [26,27], and has been reported in simulations of plasmas overdense to SRBS and at relativistic pump intensities [11,28], and has been also observed in underdense Vlasov simulations [29]. Distinguishing between BAMs and EAWs is discussed in the literature, see for instance [22,29]. We will avoid this discussion, because they play secondary roles in the results presented here. For the parameters we are using, which involves strongly damped electron plasma waves, our simulations are dominated by SRFS, SKEENS and SRBS, in that order. It is the purpose of the present work to study these three processes and their mutual interactions, SRFS, SKEENS and SRBS which arise during the SRS dynamics process, and which appear in the early stage which precedes the saturation of the SRBS. To avoid any interference from artificially distorted distribution functions or imposed seeding, we start the code from an initial Maxwellian distribution, and the system evolves under the influence of a pump light wave which provides fluctuations from which SRS develops, without any additional imposed initial perturbation except for a standing wave at the resonant wavenumber of SRFS. This then develops SRFS but also drives SKEENS at the backscattering portion beating with the pump electrostatic ponderomotive field which seeds a KEEN wave directly. We do not seed the counter-propagating daughter light wave to stimulate the growth of the SRBS. We identify in the early phase of the Raman interaction a backscattered light that phase-matches for scattering off a KEEN wave, and which precedes

saturation of SRS as well via Langmuir decay instability, etc.

254 Computational and Numerical Simulations

(the upper sign in Eq.(1) is for the electron equation and the lower sign for the ion equation, and subscripts *e* or *i* denote electrons or ions, respectively). In our normalized units *μe* =1 and *μi* =*Mi* / *Me* is the ratio of ion to electron masses.

$$E\_{\pm} = -\frac{\partial \phi}{\partial \mathbf{x}} \text{and } \vec{E}\_{\perp} = -\frac{\partial \vec{a}\_{\perp}}{\partial t} \tag{2}$$

and *ϕ* is given by Poisson's equation, which is given here by:

$$\frac{\partial^2 \phi}{\partial x^2} = \int f\_e(\mathbf{x}, p\_{xe}) dp\_{xe} - \int f\_i(\mathbf{x}, p\_{xi}) dp\_{xi} \tag{3}$$

The transverse electromagnetic fields *Ey*, *Bz* for the linearly polarized wave obey Maxwell's equations. Defining *<sup>E</sup>* <sup>±</sup> <sup>=</sup>*Ey* <sup>±</sup> *Bz*, we have:

$$(\frac{\partial}{\partial t} \pm \frac{\partial}{\partial x})E^{\pm} = -J\_y. \tag{4}$$

where (*x*(*tn*), *pxe*,*<sup>i</sup>*

,

Equations (7) and (8) can be rewritten as:

Which are implicit equations for *Δxe*,*<sup>i</sup>*

vector in Eq.(12) and *Ve*,*<sup>i</sup>* =(*Vxe*,*<sup>i</sup>*

1

*<sup>k</sup>* +1 <sup>=</sup> *<sup>Δ</sup><sup>t</sup>*

is calculated from *f <sup>e</sup>*,*<sup>i</sup>*

by writing: *Δ<sup>X</sup> <sup>e</sup>*,*<sup>i</sup>*

is the two dimensional vector *Xe*,*<sup>i</sup>* =(*x*, *pxe*,*<sup>i</sup>*

<sup>2</sup> *<sup>V</sup>e*,*<sup>i</sup>*

a grid point).

Put

*Xe*,*<sup>i</sup>*

Then *f <sup>e</sup>*,*<sup>i</sup> n*+1 (*tn*)) is the point where the characteristic is originating at *tn* (not necessarily

, ,

*Vx p* <sup>D</sup> D = -D -D (10)

*Vx p* <sup>D</sup> D = -D -D (11)

<sup>D</sup> D = <sup>D</sup> **<sup>X</sup> V X <sup>X</sup>** (12)

=(*Δxe*,*<sup>i</sup>*

*<sup>k</sup>* ,*tn*+1/2), where we start the iteration with *Δ<sup>X</sup> <sup>e</sup>*,*<sup>i</sup>*

, *Δpxe*,*<sup>i</sup>*

) is the two dimensional

http://dx.doi.org/10.5772/57476

257

<sup>0</sup> =0 for *k*=0.

is implicit and is solved iteratively

,

and are solved by iteration. This iteration is

D = D = (9)

,

*xe i j n xe ij xe i n*

*x xt* - *p pt* -

*p x*

, , , ,, ( , ). <sup>2</sup> *<sup>x</sup> p xe i xe i xe i j xe i xe ij p*

, , , , , ( , ) <sup>2</sup> *xe i xe i x p xe i <sup>p</sup> p j xe i xe ij p*

and *Δpxe*,*<sup>i</sup>*

, , , , 1/2 ( - , ) <sup>2</sup> *ei ei ei ei n <sup>t</sup> <sup>t</sup>* <sup>+</sup>

), and *Δ<sup>X</sup> <sup>e</sup>*,*<sup>i</sup>*

*e i j xe ij e i n xe i n e i j xe i xe i p f x p f xt p t f x p* <sup>+</sup> <sup>=</sup> = -D D (13)

). Eq.(12) for *Δ<sup>X</sup> <sup>e</sup>*,*<sup>i</sup>*

Usually two or three iterations are sufficient to get a good convergence. The shifted values in Eqs.(10,11) are calculated by a two-dimensional interpolation using a tensor product of cubic *B*-splines [30]. We now write that the distribution function is constant along the characteristics.

*<sup>n</sup>* from the relation :

, , , , , ,, ( , ) ( ( ), ( )) ( 2 , -2 ). *x p <sup>x</sup> xe i*

Again the shifted values in Eq.(13) are calculated with a two-dimensional interpolation using a tensor product of cubic *B*-splines. Details have been presented in [30]. These methods

*xe i p*

*t*

*t*

effected as follows. We rewrite Eqs.(10,11) in the vectorial form:

*n*+1/2

*nn n*

(*Xe*,*<sup>i</sup>*

, *V pxe*,*<sup>i</sup> n*+1/2


( ) ( ) ; . <sup>2</sup> <sup>2</sup>

Stimulated Raman Scattering with a Relativistic Vlasov-Maxwell Code: Cascades of Nonstationary Nonlinear…

In our normalized units we have the following expressions for the normal current densities:

$$
\vec{J}\_{\perp} = \vec{J}\_{\perp e} + \vec{J}\_{\perp i} \quad ; \qquad \vec{J}\_{\perp e, i} = -\frac{\hat{a}\_{\perp}}{\mu\_{e, i}} \int\_{-\omega}^{\omega} \frac{f\_{e, i}}{\mathcal{V}\_{e, i}} d\, p\_{\chi e, i}.
$$

The longitudinal electric field is calculated from Ampère's equation: ∂*Ex* / ∂*t* = − *Jx* where

$$J\_x = \frac{1}{\mu\_i} \int\_{-\infty}^{+\infty} \frac{p\_{xi}}{\mathcal{I}\_i} f\_i \, dp\_{xi} - \frac{1}{\mu\_e} \int\_{-\infty}^{+\infty} \frac{p\_{xe}}{\mathcal{I}\_e} f\_e \, dp\_{xe} \tag{5}$$

Test runs were made in which Poisson's equation was used instead of Ampère's equation to obtain the longitudinal electric field, with identical results.

The Eulerian Vlasov code we use to solve Eqs.(1-5) was recently presented and applied in [11-16] for instance. We outline the main steps for the numerical solution of Eq.(1), using an Eulerian scheme. Given *f <sup>e</sup>*,*<sup>i</sup> <sup>n</sup>* at mesh points at time *t* =*nΔt* (we stress here that the subscript *i* denotes the ion distribution function), we calculate the new value *f <sup>e</sup>*,*<sup>i</sup> n*+1 at the grid points *jx*, and *jp* corresponding to the mesh points (*x <sup>j</sup> x* , *pxe*,*<sup>i</sup> <sup>j</sup> p* ) by writing that the distribution function is constant along the characteristics. The characteristics equations for Eq.(1) are given by:

$$\begin{split} \frac{d\mathbf{x}}{dt} &= m\_{e,i} \frac{p\_{xe,i}}{\mathcal{Y}\_{e,i}} = V\_{xe,i}(\mathbf{x}, p\_{xe,i}) \\ \frac{d p\_{xe,i}}{dt} &= \mp E\_x - \frac{m\_{e,i}}{2\mathcal{Y}\_{e,i}} \frac{\partial a\_{\perp}^2}{\partial \mathbf{x}} = V\_{p\_{xe,i}}(\mathbf{x}, p\_{xe,i}). \end{split} \tag{6}$$

We assume that at the time *tn*+1 ≡*tn* + *Δt*, *x* is at the grid point *jx*, and *pxe*,*<sup>i</sup>* is at the grid point *jp*. The following leapfrog scheme can be written for the solution of (6):

$$\frac{\mathbf{x}\_{j\_x} - \mathbf{x}(t\_n)}{\Delta t} = V\_{\mathbf{x}x, i}(\mathbf{x}^{n+1/2}, p\_{\mathbf{x}x, i}^{n+1/2}) = V\_{\mathbf{x}x, i}(\frac{\mathbf{x}\_{j\_x} + \mathbf{x}(t\_n)}{2}, \frac{p\_{\mathbf{x}x, i\bar{j}\_p} + p\_{\mathbf{x}x, i}(t\_n)}{2}) \tag{7}$$

$$\frac{p\_{xe,\vec{\eta}\_p} - p\_{xe,i}(t\_n)}{\Delta t} = V\_{p\_{xe,i}}(\mathbf{x}^{n+1/2}, p\_{xe,i}^{n+1/2}) = V\_{p\_{xe,i}}(\frac{\mathbf{x}\_{\vec{\eta}\_x} + \mathbf{x}(t\_n)}{2}, \frac{p\_{xe,\vec{\eta}\_p} + p\_{xe,i}(t\_n)}{2}) \tag{8}$$

where (*x*(*tn*), *pxe*,*<sup>i</sup>* (*tn*)) is the point where the characteristic is originating at *tn* (not necessarily a grid point).

Put

() . *<sup>y</sup> E J t x* ¶ ¶ <sup>±</sup> ± =-

> *d pxe*,*<sup>i</sup>* .

*J* → <sup>⊥</sup> = *J* → <sup>⊥</sup>*<sup>e</sup>* + *J* → <sup>⊥</sup>*<sup>i</sup>* ; *J* →

256 Computational and Numerical Simulations

Eulerian scheme. Given *f <sup>e</sup>*,*<sup>i</sup>*

and *jp* corresponding to the mesh points (*x <sup>j</sup>*

<sup>⊥</sup>*e*,*<sup>i</sup>* <sup>=</sup> <sup>−</sup> *<sup>a</sup>* ⇀ ⊥ *μe*,*<sup>i</sup> <sup>∫</sup>* −*∞*

+*∞ f e*,*i γe*,*<sup>i</sup>*

mg

denotes the ion distribution function), we calculate the new value *f <sup>e</sup>*,*<sup>i</sup>*

,

*dx <sup>p</sup> m V xp dt*

*xe i*

, <sup>2</sup> , ,

*e i xe i e i*

We assume that at the time *tn*+1 ≡*tn* + *Δt*, *x* is at the grid point *jx*, and *pxe*,*<sup>i</sup>*

, ,,

*xe i xe i xe i*

*Vx p V <sup>t</sup>*

*Vx p V <sup>t</sup>*

<sup>=</sup> <sup>=</sup> <sup>D</sup>

<sup>=</sup> <sup>=</sup> <sup>D</sup>

The following leapfrog scheme can be written for the solution of (6):

g

= =

, , ,

*e i xe i xe i*

,

*e i*

g

¶ =- = ¶ <sup>m</sup>

*dp <sup>m</sup> <sup>a</sup> <sup>E</sup> V xp dt <sup>x</sup>*

obtain the longitudinal electric field, with identical results.

In our normalized units we have the following expressions for the normal current densities:

The longitudinal electric field is calculated from Ampère's equation: ∂*Ex* / ∂*t* = − *Jx* where

 mg

Test runs were made in which Poisson's equation was used instead of Ampère's equation to

The Eulerian Vlasov code we use to solve Eqs.(1-5) was recently presented and applied in [11-16] for instance. We outline the main steps for the numerical solution of Eq.(1), using an

> *x* , *pxe*,*<sup>i</sup> <sup>j</sup> p*

is constant along the characteristics. The characteristics equations for Eq.(1) are given by:

(, )

*x p xe i*

^

( ) ( ) ( )

*<sup>p</sup> <sup>x</sup> <sup>x</sup> j n j n xe ij xe i n n n*

*xe ij xe i n j n xe ij xe i n n n*

*p pt x xt p pt*

+ + - <sup>+</sup> <sup>+</sup>

*x xt x xt p pt*

+ + - <sup>+</sup> <sup>+</sup>

, , , , , , 1/2 1/2 ,

*p xe i p*

*p p x xe i xe i*

,

( , ). <sup>2</sup> *xe i*

, , 1/2 1/2

( ) ( ) ( )

(,) ( , ) <sup>2</sup> <sup>2</sup>

(,) ( , ) <sup>2</sup> <sup>2</sup>

1 1 *xi xe x i xi e xe i i e e <sup>p</sup> <sup>p</sup> <sup>J</sup> f dp f dp*

+¥ +¥


¶ ¶ (4)

<sup>=</sup> - ò ò (5)

*<sup>n</sup>* at mesh points at time *t* =*nΔt* (we stress here that the subscript *i*

,

*n*+1

) by writing that the distribution function

at the grid points *jx*,

is at the grid point *jp*.

(6)

(7)

(8)

$$\Delta\_{\mathbf{x}e,i} = \frac{\mathbf{x}\_{j\_x} - \mathbf{x}(t\_n)}{\mathbf{2}} \quad ; \quad \Delta\_{p\_{\mathbf{x}e,i}} = \frac{p\_{\mathbf{x}e,i j\_p} - p\_{\mathbf{x}e,i}(t\_n)}{\mathbf{2}}.\tag{9}$$

Equations (7) and (8) can be rewritten as:

$$
\Delta\_{xe,i} = \frac{\Delta t}{2} V\_{xe,i} (\mathbf{x}\_{j\_x} - \Delta\_{xe,i}, p\_{xe,ij\_p} - \Delta\_{p\_{xe,i}}).\tag{10}
$$

$$
\Delta\_{p\_{\rm av,i}} = \frac{\Delta t}{2} V\_{p\_{\rm av,i}} (\mathbf{x}\_{j\_x} - \Delta\_{\rm xe,i}, p\_{\rm xe,j\_p} - \Delta\_{p\_{\rm xe,i}}) \tag{11}
$$

Which are implicit equations for *Δxe*,*<sup>i</sup>* and *Δpxe*,*<sup>i</sup>* and are solved by iteration. This iteration is effected as follows. We rewrite Eqs.(10,11) in the vectorial form:

$$
\Delta\_{\mathbf{X}\varepsilon,i} = \frac{\Delta t}{2} \mathbf{V}\_{e,i} \left( \mathbf{X}\_{e,i} \mathbf{\bar{-}} \Delta\_{\mathbf{X}e,i\prime} t\_{n+1/2} \right) \tag{12}
$$

*Xe*,*<sup>i</sup>* is the two dimensional vector *Xe*,*<sup>i</sup>* =(*x*, *pxe*,*<sup>i</sup>* ), and *Δ<sup>X</sup> <sup>e</sup>*,*<sup>i</sup>* =(*Δxe*,*<sup>i</sup>* , *Δpxe*,*<sup>i</sup>* ) is the two dimensional vector in Eq.(12) and *Ve*,*<sup>i</sup>* =(*Vxe*,*<sup>i</sup> n*+1/2 , *V pxe*,*<sup>i</sup> n*+1/2 ). Eq.(12) for *Δ<sup>X</sup> <sup>e</sup>*,*<sup>i</sup>* is implicit and is solved iteratively by writing: *Δ<sup>X</sup> <sup>e</sup>*,*<sup>i</sup> <sup>k</sup>* +1 <sup>=</sup> *<sup>Δ</sup><sup>t</sup>* <sup>2</sup> *<sup>V</sup>e*,*<sup>i</sup>* (*Xe*,*<sup>i</sup>* -*Δ<sup>X</sup> <sup>e</sup>*,*<sup>i</sup> <sup>k</sup>* ,*tn*+1/2), where we start the iteration with *Δ<sup>X</sup> <sup>e</sup>*,*<sup>i</sup>* <sup>0</sup> =0 for *k*=0. Usually two or three iterations are sufficient to get a good convergence. The shifted values in Eqs.(10,11) are calculated by a two-dimensional interpolation using a tensor product of cubic *B*-splines [30]. We now write that the distribution function is constant along the characteristics. Then *f <sup>e</sup>*,*<sup>i</sup> n*+1 is calculated from *f <sup>e</sup>*,*<sup>i</sup> <sup>n</sup>* from the relation :

$$f\_{e,i}^{n+1}(\mathbf{x}\_{\mathbf{j}\_x}, p\_{xe,\vec{\mu}\_{\mathbf{j}\_y}}) = f\_{e,i}^n(\mathbf{x}(t\_n), p\_{xe,i}(t\_n)) = f\_{e,i}^n(\mathbf{x}\_{\mathbf{j}\_x} - 2\Delta\_{xe,i}p\_{xe,i} - 2\Delta\_{p\_{xe,i}}).\tag{13}$$

Again the shifted values in Eq.(13) are calculated with a two-dimensional interpolation using a tensor product of cubic *B*-splines. Details have been presented in [30]. These methods compared favourably with other Eulerian methods for the numerical solution of the Vlasov equation [31].

The numerical scheme to advance Eq.(1) from time *tn* to *tn*+1 necessitates the knowledge of the electromagnetic field *E* <sup>±</sup> at time *tn*+1/2. This is done using a centered scheme where we integrate Eq.(4) exactly along the vacuum characteristics with *Δx* =*Δt*, to calculate *E* <sup>±</sup>*n*+1/2 as follows:

$$E^{\pm}(\mathbf{x} \pm \Delta t, t\_{n+1/2}) = E^{\pm}(\mathbf{x}, t\_{n-1/2}) - \Delta t \mathbf{J}\_y(\mathbf{x} \pm \Delta t \,/\ \mathbf{2}, t\_n) \tag{14}$$

*a*0 <sup>2</sup> <sup>=</sup> *<sup>I</sup>λ*<sup>0</sup>

*ω*0 <sup>2</sup> <sup>=</sup>*ωpe* <sup>2</sup> <sup>+</sup> *<sup>k</sup>*<sup>0</sup> 2 *c* 2

normalize time to *ω*<sup>0</sup>

<sup>2</sup> / 1.368*x*1018, *I* is the laser intensity in W/cm2

*<sup>E</sup>*<sup>0</sup> <sup>=</sup>*ω*0*a*0 in our units and *<sup>E</sup>* <sup>−</sup> =0 (no small seed is applied as *<sup>E</sup>* <sup>−</sup>

(normalized to *ωpe*). Hydrogen ions are used with *Mi*

−1

with the SRFS are roots of the equation [32]:

w

w ww*sB eB* m

parameters *<sup>μ</sup>* <sup>=</sup>*mec* <sup>2</sup> / *<sup>κ</sup>Te* <sup>=</sup>*<sup>c</sup>* <sup>2</sup> / *<sup>υ</sup>te*

by *ω*0, so that in this case we would have *E*<sup>0</sup> =*a*0).

, or in normalized units *ω*<sup>0</sup>

The frequencies are normalized to the plasma frequency *ωpe*, and the pump wave frequency *ω*0 of the injected laser beam is such that *ω*<sup>0</sup> / *ωpe* =1 / *n* / *ncr*, which corresponds to *ω*<sup>0</sup> =3.481

Stimulated Raman Scattering with a Relativistic Vlasov-Maxwell Code: Cascades of Nonstationary Nonlinear…

polarized wave is injected in the domain at the left boundary at *x*=0 with *<sup>E</sup>* <sup>+</sup> =2*E*0cos(*ω*0*t*),

normalized time and length should be multiplied by *ω*0, and the electric field should be divided

The frequency and wavenumber (*ω*0, *k*0) of the pump wave are related by the relation

<sup>2</sup> =1 <sup>+</sup> *<sup>k</sup>*<sup>0</sup> 2

coupling of a pump light wave to a daughter light wave and an electron plasma wave, the values of the electron plasma wavenumber *keB* associated with the SRBS, and *keF* associated

<sup>ë</sup> <sup>û</sup> <sup>ë</sup> <sup>û</sup> (15)

with *K* =*keλDe* and *Ω* =*ω*0 (normalized to *ωpe*). For the present problem we have the following

*keBλDe* =0.3377 for the plasma mode associated with the SRBS, and *keFλDe* =0.0666 for the plasma mode associated with the SRFS. As discussed in Bers *et al*., 2009, for these parameters the SRBS plasma wave is heavily damped, and the damping of the SRFS plasma wave is negligible. The heavily damped regime with *k λDe* >0.29 is called the kinetic regime (Kline *et al.*, 2005). In our normalized units the Debye length *λDe* =*υte* / *c* =0.04424 *Te* (normalized to *c* / *ωpe* in our units), so *λDe* =0.06256 for *Te* = 2 keV. We finally get *keB* =0.3377 / *λDe* =5.398 for the SRBS plasma wave, and *keF* =0.0666 / *λDe* =1.0645 for the SRFS plasma wave. The corresponding frequencies for the

2

m

( ) 1/2 <sup>4</sup> 2 1/2 2 <sup>é</sup> 15 / 4 6 ( 3 3) 2 1 (5 / 2 ) 2 1 (5 / 2 )( 1) 0 W - + + W- - W - + <sup>ù</sup>*K K*

<sup>2</sup> =1 / (0.04424 *Te*)

SRBS plasma wave and the SRFS plasma wave are solutions of the equation [32]:

 lw

(*ωsF* , *ksF* ) and the backward scattered electromagnetic wave (*ωsB*, *ksB*):

<sup>2</sup> 22 2 44 4 1 3 / 15 / 5 / (2 ) *De De*

 lw

Equation (16) has the following roots: *ωeB* =1.178 for the SRBS and *ωeF* =1.0066 for the SRFS. The selection rules give the following results for the forward scattered electromagnetic wave

<sup>0</sup> <sup>0</sup> 3.481 1.178 2 ;.303 3.481 1.0066 2.4744

w ww*sF eF* =- = - = =- = - = (17)

 m» + *k* + *k* - (16)

mé

and length to *c* / *ω*0, then in the results we present in this paper the

, and *λ*0 the laser wavelength in microns.

*/Me=*1836. A forward propagating linearly

, from which *k*<sup>0</sup> =3.3343. For SRS, or the

m

=255.8 for *Te* = 2 keV. The resulting roots are

ù *K* + W- - W- =

). (Note that if we choose to

http://dx.doi.org/10.5772/57476

259

with *Jy*(*x* ± *Δt* / 2,*tn*) = *Jy*(*x* ± *Δx*, *tn*) + *Jy*(*x*, *tn*) 2

From Eq.(2) we also have *a* → ⊥ *<sup>n</sup>*+1 =*a* → ⊥ *<sup>n</sup>* −*ΔtE* → ⊥ *n*+1/2 , from which we calculate *a* → ⊥ *<sup>n</sup>*+1/2 =(*a* → ⊥ *<sup>n</sup>*+1 + *a* → ⊥ *<sup>n</sup>* ) / 2. To calculate *Ex n*+1/2 , we use Ampère's equation: ∂*Ex* <sup>∂</sup>*<sup>t</sup>* <sup>=</sup> <sup>−</sup> *Jx*, from which *Ex <sup>n</sup>*+1/2 <sup>=</sup>*Ex <sup>n</sup>*−1/2 <sup>−</sup>*ΔtJx n*.

### **3. The relevant parameters**

We use a fine resolution grid in phase-space, with *N* = 60000 grid points in space, and 512 grid points in momentum space for the electrons and 256 grid points in momentum space for the ions (extrema of the electron momentum are ±0.5, and ±15 for the ion momentum). The initial distribution functions for the electrons and the ions are Maxwellian. The maximum of the density is normalized to *n=*0.0825*ncr*, where *ncr* is the critical density. The electron temperature is *Te* =2 keV. The ions have a temperature *Ti* = 0.5 keV. Ions are allowed to move, especially to adjust the sheath structure at the boundaries on both sides, but we noted at the end of the simulations beginning traces of stimulated Brillouin backscattering, which remained at a very weak level, and therefore ion dynamics can be ignored in the results we are presenting [This is not logical. Ion dynamics affects SRS saturation by creating IAW mediated LDI and similar instabilities. Given that EPW and KEEN waves live in these simulations, it is incorrect to assume a priori that IAWs could play no role. The only exception to this would be that you ran for a very short time in which case the results are not demonstrative of real situations which occur over 100s of ps at the very least at these intensities]. The initial flat profile of the uniform plasma with the density *ne=ni =*1 (normalized to *n*) extends over a length *L <sup>p</sup>* =1122.56*c* / *ωpe*. On either side of the slab the densities are smoothly brought down to zero through a parabolic profile of length *L edge* =7.81*c* / *ωpe*. An extra vacuum region of length *L vac* =16.81*c* / *ωpe* exists on each side of the slab, for a total length of the system of *L* =1171.89*c* / *ωpe*. In our normalized units *Δx* =*Δt*.

A characteristic parameter of laser beams is the normalized vector potential or quiver mo‐ mentum |*a* → <sup>⊥</sup> | = |*eA* → <sup>⊥</sup> / *Mec* | =*a*0, where *A* → <sup>⊥</sup> is the vector potential of the wave. We chose for the amplitude of the vector potential *a*<sup>0</sup> =0.025. For the linearly polarized wave *a*0 <sup>2</sup> <sup>=</sup> *<sup>I</sup>λ*<sup>0</sup> <sup>2</sup> / 1.368*x*1018, *I* is the laser intensity in W/cm2 , and *λ*0 the laser wavelength in microns. The frequencies are normalized to the plasma frequency *ωpe*, and the pump wave frequency *ω*0 of the injected laser beam is such that *ω*<sup>0</sup> / *ωpe* =1 / *n* / *ncr*, which corresponds to *ω*<sup>0</sup> =3.481 (normalized to *ωpe*). Hydrogen ions are used with *Mi /Me=*1836. A forward propagating linearly polarized wave is injected in the domain at the left boundary at *x*=0 with *<sup>E</sup>* <sup>+</sup> =2*E*0cos(*ω*0*t*), *<sup>E</sup>*<sup>0</sup> <sup>=</sup>*ω*0*a*0 in our units and *<sup>E</sup>* <sup>−</sup> =0 (no small seed is applied as *<sup>E</sup>* <sup>−</sup> ). (Note that if we choose to normalize time to *ω*<sup>0</sup> −1 and length to *c* / *ω*0, then in the results we present in this paper the normalized time and length should be multiplied by *ω*0, and the electric field should be divided by *ω*0, so that in this case we would have *E*<sup>0</sup> =*a*0).

compared favourably with other Eulerian methods for the numerical solution of the Vlasov

The numerical scheme to advance Eq.(1) from time *tn* to *tn*+1 necessitates the knowledge of the

∂*Ex*

We use a fine resolution grid in phase-space, with *N* = 60000 grid points in space, and 512 grid points in momentum space for the electrons and 256 grid points in momentum space for the ions (extrema of the electron momentum are ±0.5, and ±15 for the ion momentum). The initial distribution functions for the electrons and the ions are Maxwellian. The maximum of the density is normalized to *n=*0.0825*ncr*, where *ncr* is the critical density. The electron temperature

adjust the sheath structure at the boundaries on both sides, but we noted at the end of the simulations beginning traces of stimulated Brillouin backscattering, which remained at a very weak level, and therefore ion dynamics can be ignored in the results we are presenting [This is not logical. Ion dynamics affects SRS saturation by creating IAW mediated LDI and similar instabilities. Given that EPW and KEEN waves live in these simulations, it is incorrect to assume a priori that IAWs could play no role. The only exception to this would be that you ran for a very short time in which case the results are not demonstrative of real situations which occur over 100s of ps at the very least at these intensities]. The initial flat profile of the uniform

either side of the slab the densities are smoothly brought down to zero through a parabolic profile of length *L edge* =7.81*c* / *ωpe*. An extra vacuum region of length *L vac* =16.81*c* / *ωpe* exists on each side of the slab, for a total length of the system of *L* =1171.89*c* / *ωpe*. In our normalized

A characteristic parameter of laser beams is the normalized vector potential or quiver mo‐

the amplitude of the vector potential *a*<sup>0</sup> =0.025. For the linearly polarized wave

→

<sup>⊥</sup> / *Mec* | =*a*0, where *A*

Eq.(4) exactly along the vacuum characteristics with *Δx* =*Δt*, to calculate *E* <sup>±</sup>*n*+1/2

*Jy*(*x* ± *Δx*, *tn*) + *Jy*(*x*, *tn*) 2

> → ⊥ *<sup>n</sup>*+1 =*a* → ⊥ *<sup>n</sup>* −*ΔtE* → ⊥ *n*+1/2

, we use Ampère's equation:

1/2 1/2 ( , ) ( , ) ( / 2, ) *<sup>n</sup> <sup>n</sup> <sup>y</sup> <sup>n</sup> E x t t E x t tJ x t t* <sup>±</sup> <sup>±</sup>

at time *tn*+1/2. This is done using a centered scheme where we integrate

<sup>+</sup> - ±D = -D ±D (14)

, from which we calculate *a*

<sup>∂</sup>*<sup>t</sup>* <sup>=</sup> <sup>−</sup> *Jx*, from which *Ex*

→ ⊥ *<sup>n</sup>*+1/2 =(*a* → ⊥ *<sup>n</sup>*+1 + *a* → ⊥ *<sup>n</sup>* ) / 2. To

= 0.5 keV. Ions are allowed to move, especially to

<sup>⊥</sup> is the vector potential of the wave. We chose for

*=*1 (normalized to *n*) extends over a length *L <sup>p</sup>* =1122.56*c* / *ωpe*. On

*<sup>n</sup>*+1/2 <sup>=</sup>*Ex*

*<sup>n</sup>*−1/2 <sup>−</sup>*ΔtJx*

*n*.

as follows:

equation [31].

electromagnetic field *E* <sup>±</sup>

258 Computational and Numerical Simulations

with *Jy*(*x* ± *Δt* / 2,*tn*) =

calculate *Ex*

From Eq.(2) we also have *a*

*n*+1/2

**3. The relevant parameters**

plasma with the density *ne=ni*

units *Δx* =*Δt*.

mentum |*a*

→

<sup>⊥</sup> | = |*eA* →

is *Te* =2 keV. The ions have a temperature *Ti*

The frequency and wavenumber (*ω*0, *k*0) of the pump wave are related by the relation *ω*0 <sup>2</sup> <sup>=</sup>*ωpe* <sup>2</sup> <sup>+</sup> *<sup>k</sup>*<sup>0</sup> 2 *c* 2 , or in normalized units *ω*<sup>0</sup> <sup>2</sup> =1 <sup>+</sup> *<sup>k</sup>*<sup>0</sup> 2 , from which *k*<sup>0</sup> =3.3343. For SRS, or the coupling of a pump light wave to a daughter light wave and an electron plasma wave, the values of the electron plasma wavenumber *keB* associated with the SRBS, and *keF* associated with the SRFS are roots of the equation [32]:

$$\mathbb{E}\left[\left(15\Omega\,/4-6\right)\right]\mathbf{K}^4 + \left(\mu + 3\Omega - 3\right)\mathbf{K}^2 - 2\mu^{1/2}\left[\Omega^2 - 1 + \left(5/2\,\mu\right)\right]^{1/2}\mathbf{K} + 2\Omega - 1 - \left(5/2\,\mu\right)\left(\Omega - 1\right) = 0\tag{15}$$

with *K* =*keλDe* and *Ω* =*ω*0 (normalized to *ωpe*). For the present problem we have the following parameters *<sup>μ</sup>* <sup>=</sup>*mec* <sup>2</sup> / *<sup>κ</sup>Te* <sup>=</sup>*<sup>c</sup>* <sup>2</sup> / *<sup>υ</sup>te* <sup>2</sup> =1 / (0.04424 *Te*) 2 =255.8 for *Te* = 2 keV. The resulting roots are *keBλDe* =0.3377 for the plasma mode associated with the SRBS, and *keFλDe* =0.0666 for the plasma mode associated with the SRFS. As discussed in Bers *et al*., 2009, for these parameters the SRBS plasma wave is heavily damped, and the damping of the SRFS plasma wave is negligible. The heavily damped regime with *k λDe* >0.29 is called the kinetic regime (Kline *et al.*, 2005). In our normalized units the Debye length *λDe* =*υte* / *c* =0.04424 *Te* (normalized to *c* / *ωpe* in our units), so *λDe* =0.06256 for *Te* = 2 keV. We finally get *keB* =0.3377 / *λDe* =5.398 for the SRBS plasma wave, and *keF* =0.0666 / *λDe* =1.0645 for the SRFS plasma wave. The corresponding frequencies for the SRBS plasma wave and the SRFS plasma wave are solutions of the equation [32]:

$$
\rho o^2 \approx 1 + 3k^2 \lambda\_{\rm De}^2 / o^2 + 15k^4 \lambda\_{\rm De}^4 \,/\, o^4 - 5 \,/\, \text{(2\,\mu)}\tag{16}
$$

Equation (16) has the following roots: *ωeB* =1.178 for the SRBS and *ωeF* =1.0066 for the SRFS. The selection rules give the following results for the forward scattered electromagnetic wave (*ωsF* , *ksF* ) and the backward scattered electromagnetic wave (*ωsB*, *ksB*):

$$
\rho\_{\text{sB}} = \alpha\_0 - \alpha\_{\text{eB}} = 3.481 - 1.178 = 2.303; \text{ } \alpha\_{\text{sF}} = \alpha\_0 - \alpha\_{\text{eF}} = 3.481 - 1.0066 = 2.4744 \tag{17}
$$

$$\begin{aligned} k\_{s\text{B}} &= k\_{e\text{B}} - k\_0 = 5.398 - 3.3343 = 2.0637; \\ k\_{s\text{F}} &= k\_0 - k\_{e\text{F}} = 3.3343 - 1.0645 = 2.2698 \end{aligned} \tag{18}$$

(0,0, ) *<sup>y</sup> <sup>B</sup> <sup>E</sup> t x* ¶ ¶ = - ¶ ¶

<sup>→</sup> =(0,*py*,0), with *py* <sup>=</sup> <sup>−</sup> *<sup>p</sup>*0sin(*ψ*), and *p*<sup>0</sup> <sup>=</sup>*E*<sup>0</sup> / *<sup>ω</sup>*. The longitudinal Lorentz force is

Stimulated Raman Scattering with a Relativistic Vlasov-Maxwell Code: Cascades of Nonstationary Nonlinear…

We note in Figs.(1c) a mode with a wavenumber 5.616. These results are confirmed in Figure (3), where we present the spatial Fourier modes at the time *t*=527 in the same domain *x* ∈(250,410), and where the mode with a wavenumber 5.616 appears with its growing harmonics at 11.232 and 16.81. This mode at *kKEEN* =5.616 is different from the value of *keB* =5.398 for the SRBS plasma wave, and will be further discussed and identified as a KEEN wave, responsible for the small vortices we see in Figure (2a) and Figure (4,top left). We show on a logarithmic scale in Figure (1b) the distribution function around *pxe* =0.183, spatially averaged over a length *λKEEN* =2*π* / *kKEEN* =1.118 (which is the width of the small vortices we see in Figure (2a) and Figure (4,top left)) around *x*=280. At this stage at time *t*=351, it shows on a logarithmic scale the straight line of a Maxwellian. However, in Figure (2b) at time *t*=468, and in Figure (4,bottom right) at time *t*=761, the distribution function shows a slightly distorted (but not fully flattened) distribution function, which lower the damping rate and which can facilitate the excitation of the mode around *pxe* =0.183. The entire distribution function, spatially averaged over a length *λKEEN* =1.118 around *x*=280, is shown at time *t*=761 in Figure (4,bottom left). Figure (4,top right) shows a plot of the longitudinal electric field in *x* ∈(280,300), showing a wavelength *λeF* =2*π* / *keF* =5.925 with a small modulation with

Figure (5) presents the spatial Fourier spectrum at time *t*=761 in *x* ∈(250,410). We note in Figure (5a) for the longitudinal wave, the mode with a wavenumber 5.616 has now developed important harmonics at 11.232 and 16.81. We also identify the dominant mode with *keF* =1.06 of the SRFS plasma wave. Since we have a linear polarization, we have in the longitudinal

At this stage, we look in Figure (5b) to the wavenumber spectrum of the forward electromag‐

at *k*<sup>0</sup> =3.338 (3.334 in our theoretical results). We can also identify the contribution of the SRFS plasma wave at *ksF* =2.277 (2.2698 in our theoretical results in Eq.(18)). We have also a small peak at *kAS* =4.398, which corresponds to the anti-Stokes coupling *kAS* =*k*<sup>0</sup> + *ke*−*AS* = 3.334 + 1.064 = 4.398 in our theoretical results (*ke*−*AS* =1.064 is the plasma wavenumber for the anti-Stokes coupling, already present in the wide dominant peak *keF* =1.06 in Figure (5a)). The frequency *ωAS* =4.51 of the anti-Stokes wave is calculated from the relation

in the domain *x* ∈(250,410). We identify the dominant pump wave, appearing

<sup>2</sup> , in close agreement with the value calculated from the relation *ωAS* <sup>=</sup>

perturbation a mode present with 2*k*<sup>0</sup> =6.676 (at twice the wavenumber of the pump).

sin(2*ψ*). This drives a longitudinal response at the 2nd harmonic of the laser wave.

→

<sup>⊥</sup> = −∂*a* →

<sup>⊥</sup> / ∂*t* and *p*

→ <sup>⊥</sup> =*a* →

http://dx.doi.org/10.5772/57476

(19)

261

<sup>⊥</sup>, we get

r

<sup>→</sup> =(0,0,*Bz*) with *Bz* <sup>=</sup> *<sup>B</sup>*0cos(*ψ*), and *B*<sup>0</sup> <sup>=</sup>*E*0*<sup>k</sup>* / *<sup>ω</sup>*. From *<sup>E</sup>*

Then *B*

*pyBz* <sup>=</sup> <sup>−</sup> <sup>1</sup>

*λKEEN* =1.118.

netic wave *E* <sup>+</sup>

<sup>2</sup> =1 <sup>+</sup> *kAS*

*ωAS*

<sup>2</sup> *<sup>k</sup> <sup>p</sup>*<sup>0</sup> 2

*p*

The results in Eqs.(17-18) obey the dispersion relation for the electromagnetic wave: 1 + *ksF* <sup>2</sup> =6.152=*ωsF* 2 (from which we get *ωsF* =2.480), and 1 <sup>+</sup> *ksB* <sup>2</sup> =5.2588=*ωsB* 2 (from which we get 2.293). These results are very close to what is calculated in Eq.(17).
