Ion Acceleration and Plasma Jet Formation, Amplification of a Whistler Pulse along a Magnetic Field, Interaction of a Relativistic Electron Beam with a Plasma: Eulerian Vlasov Codes for Laboratory and Space Plasmas Simulation

Magdi Shoucri^*

Institut de Recherche d’Hydro-Québec (IREQ), Varennes, Canada.

DOI: 10.4236/jhepgc.2025.114096                

Abstract

The application of Eulerian codes for the numerical solution of the kinetic equations of plasmas has become a powerful numerical tool. In the present work, we use an Eulerian Vlasov code for the numerical solution of the one-dimensional (1D) relativistic Vlasov-Maxwell set of equations to illustrate the application of these methods to problems pertinent to laboratory and space plasmas. In Section 2, we present a study on the acceleration of ions and plasma jet formation during the interaction of a high intensity circularly polarized laser beam normally incident on an overdense plasma target. Both electrons and ions are treated with a relativistic Vlasov equation. We follow the evolution of this system when the ponderomotive pressure of the wave pushes the electrons at the plasma surface and leads to the formation of a neutral plasma jet ejected from the back of the target. In Section 3, we use a 1D relativistic Vlasov code in a magnetized plasms to study the amplification of an ultra-short seed pulse mediated via a Brillouin backscattering of energy from a long pump pulse. In both Section 2 and Section 3, details on the numerical codes will be presented. In Section 4, we discuss the generation of electromagnetic radiation in a relativistic electron beam-plasma interaction known as the Type III solar radio bursts. In Section 5, we discuss the problem of the auroral kilometric radiation, the colourful dance of light in the sky in the Northern and Southern hemispheres of the Earth, when weakly relativistic electron beams are propagating along the magnetic field in the north and south poles. The absence of noise in the Eulerian Vlasov code allows a detailed representation of the phase-space structures of the distribution functions, especially in the low density regions of the phase-space.

Keywords

Ion Acceleration, Whistler Amplification, Relativistic Vlasov-Maxwell

1. Introduction

Transport of charged particles and collisions are main aspects of fundamental kinetic processes which we need to study in electronic devices and plasmas. Considerable amount of effort is devoted for the application of numerical techniques to the study and comprehension of these kinetic processes, which are generally highly nonlinear problems.

There are two important numerical approaches which are generally used in these studies. The first approach includes essentially Lagrangian methods which have become known as the particle-in-cell (PIC) methods. In this case, the plasma is approximated by a finite number of macro-particles which are advanced in time in the self-consistent fields computed on a background mesh. This approach has been widely developed, with successful applications to several, if not all, aspects of plasma phenomena and devices. An important drawback with this approach is that there is a numerical noise inherently associated with PIC methods. This problem becomes important when the physics of interest is in the low density regions of the phase-space where there is only a relatively small fraction of the macro-particles, as for instance in the high energy tail of the distribution function. An alternative approach is to look for the direct calculation of the particle distribution function $f\left(r,\upsilon ,t\right)$ , where $r$ and $\upsilon$ are respectively the position and velocity vectors at time t. The methods which discretize the full phase-space on a multi-dimensional grid and perform the direct solution of the distribution function on this grid, by solving numerically the partial differential equation for $f\left(r,\upsilon ,t\right)$ in the multi-dimensional phase-space, are called Eulerian methods. In the last few decades there have been important developments in the application of these methods for the numerical solutions of the kinetic equations of plasmas. Of special interest are the Eulerian codes based on the method of characteristics for hyperbolic differential equations, because they have a low noise level and they allow an accurate study of the phase-space dynamics, especially in the low density region of the phase-space. This approach has become increasingly important, especially for the numerical solution of the Vlasov equation which is used to study the kinetic processes when the collisions between the particles are negligible, as in high temperature and low-density plasmas. A review of some of the progress work done for the application of these methods is presented and reviewed in [1]. A great advantage of the Eulerian codes is their very low noise level, which allows accurate representation of the fine scale structures associated with the particles distribution function, especially in the low density regions of the phase-space. The areas of the kinetic processes in plasmas are broad, and of special interest are the applications to the study of collisionless transport processes using the Vlasov equation for a wide range of phenomena displaying collective behaviour, and in different many-particles systems in plasmas, charged-particles beams, laser-plasma interactions, inviscid fluid turbulence, and in a variety of other domains such as galaxies and other self-gravitating systems, and more recently in semiconductors and thin films.

We present in this paper few examples of the application of the Eulerian methods to the kinetic equations of plasmas for the numerical solution of the appropriate relativistic Vlasov-plasma equations. Two problems are discussed respectively in Sections 2 and 3. For each of these problems, the details of the numerical methods developed will be presented. The first example we present in Section 2 is for the problem of the ions acceleration and plasma jet formation in the interaction of an intense laser beam normally obtained from these simulations which could have relevance in different fields, in a number of physical situations in laboratory plasma where plasma jets are observed [2]-[4], and more recently in ion-driven fast ignition [5]. And in astrophysics, it is believed to provide the main acceleration mechanism for cosmic rays at the front of collisionless shock waves generated during supernova explosions by the interaction of high intensity electromagnetic waves with plasmas [6]. Detailed results of the phase-space evolution for the case we study will be presented. Reference [7] contains additional results for related problems. Section 3 presents the applications for the numerical solution of a relativistic Vlasov-Maxwell equations, for a study of the problem of stimulated Brillouin backscattering in the interaction of a laser beam with a plasma. This problem is important for the application of the idea of plasma based laser amplifiers for the direct amplification of ultra-short laser pulses. It has been studied in unmagnetized plasmas [8], and in a plasma embedded in an external magnetic field [9] [10], to study the transfer of energy from a long pump electromagnetic wave, to a counter-propagating ultra-short seed pulse, mediated by a resonant ion wave via stimulated Raman or Brillouin backscattering. We study in Section 3 a case which includes the presence of an applied external magnetic field for the Brillouin backscattering, which also allows the amplification to whistler pulses, a problem of importance for the amplification of whistler waves along the magnetic field of the Earth. In addition to the problem of the amplification of ultra-short seed pulses, there is an obvious academic interest in the study of problems of the excitation of electromagnetic waves observed in many situations in laboratory plasmas and space plasmas in the magnetosphere and other geophysical situations in the environments of planets, where important variations in the presence and strength of magnetic fields are observed. Section 4 will study the problem of Type III solar radio bursts excited when relativistic bunches of electron beams are ejected from the sun [11]. And in Section 5, we discuss the problem of the auroral kilometric radiation [12] [13], the colourful dance of lights in the sky in the Northern and Southern hemispheres of the Earth, when weakly relativistic electron beams are precipitating along the magnetic field in the north and south poles.

2. Ions Acceleration and Plasma Jet Formation in the Interaction of an Intense Laser Beam Normally Incident on an Overdense Plasma: A Vlasov Code Simulation

The field of the interaction of a high intensity (>10²² W/cm²) and extremely short femtoseconds (10⁻¹⁵ s) laser beam [14] [15] incident on foil targets has attracted considerable attention in recent years. Under these conditions, a solid target transforms very rapidly, in a few cycles of the light wave, into a plasma. During the short period of the laser pulse, the laser interacts directly with a plasma with density of the order of the density of the solid. This interaction takes place in the relativistic domain because of the high intensity of the laser field. Collisions become negligible since the electrons quiver energy is relativistic, and the interaction at the plasma surface will be dominated by collective effects. Interesting interactions between the wave and the plasma take place, with promising applications in a variety of areas in physics and medicine have been demonstrated [16]-[19]. Because of the production of high energy accelerated ions to multi-MeV energies, especially in the case of a circular polarization of the laser wave, the phenomenon of energization and expansion of particles at the front of a high intensity wave and the subsequent formation and evolution of collisionless shock waves have been of particular interest since they play a fundamental role in laboratory laser-plasma interactions, and in a number of physical situations in laboratory and space plasmas, where plasma jets of high speed are observed [2]-[4], and more recently in ion-driven fast ignition [5]. It is believed to provide the main acceleration mechanism for cosmic rays at the front of collisionless shock waves generated during supernova explosions [6]. Also of great interest is the generation of coherent high order harmonics when the laser wave is linearly polarized, which has important potential applications as radiation sources with attosecond (10⁻¹⁸ s) duration and with unprecedented intensities, offering a combination of short wavelengths and very high time resolution [20].

In the present Section, we study the interaction of a high intensity laser beam normally incident on an overdense plasma, using an Eulerian Vlasov code for the numerical solution of the 1D relativistic Vlasov-Maxwell equations for both electrons and ions [21]-[29]. If the intensity of the wave is sufficiently high to make the oscillation of the electrons relativistic, interesting interactions between the wave and the surface of the plasma take place. Of particular interest for the problem of the acceleration of the ions is the case of circular polarization of the laser wave, when the free space wavelength of the wave ${\lambda }_0$ is greater than the scale length of the jump in the plasma density at the plasma edge $L_{edge}$ (${\lambda }_0>L_{edge}$ ), and the plasma density is such that $n/n_{cr}\gg 1$ , where $n_{cr}=1.1\times {10}^{21}{\lambda }_0^{-2}{\text{cm}}^{-3}$ (${\lambda }_0$ is the laser wavelength expressed in microns). In these regimes, the radiation pressure causes a charge separation and an electric field to appear at the plasma edge, and produces a plasma surface with a steep density profile. In the case of circular polarization of the laser wave, there is a non-oscillatory radiation pressure applied at the target surface by the incident laser wave. This radiation pushes the electrons at the plasma surface through the ponderomotive pressure, producing a sharp density gradient at the target plasma surface. There is a build-up of the electron density at this surface which creates a large space-charge, giving rise to a strong longitudinal electric field. The combined effects of the edge electric field and the radiation pressure results in an important acceleration of the ions in the forward direction. In the ongoing acceleration process, if the laser beam intensity is a Gaussian-shaped pulse, the accelerated ions reach a free-streaming expansion phase where the electrons neutralize the charge of the expanding ions, and a neutral plasma jet is formed. When the Gaussian shaped pulse decays, a dense compact bunch of a neutral plasma jet is formed and propagates at a constant speed in the forward direction.

In the present Section, we present numerical simulations to study the case of a 1D model where a circularly polarized laser beam pulse with a Gaussian-shaped intensity, is normally incident on the surface of a target. We will follow the formation of a free-streaming dense compact bunch of a neutral plasma jet.

2.1. The 1D Relativistic Vlasov-Maxwell Model

Time t is normalized to the inverse electron plasma frequency ${\omega }_{pe}^{-1}$ , length is normalized to $l_0=c{\omega }_{pe}^{-1}$ , velocity and momentum are normalized respectively to the velocity of light c and to $m_{e}c$ . The general form of the Vlasov equation is written for the present problem first in a 4D phase-space for the electron distribution function $F_e\left(x,p_{xe},p_{ye},p_{ze},t\right)$ and the ion distribution function $F_i\left(x,p_{xi},p_{yi},p_{zi},t\right)$ (one spatial dimension) as follows:

$\frac{\partial F_{e,i}}{\partial t}+M_{e,i}\frac{p_{xe,i}}{{\gamma }_{e,i}}\frac{\partial F_{e,i}}{\partial x}\mp \left(E+\frac{p_{e,i}\times B}{{\gamma }_{e,i}}\right)\cdot \frac{\partial F_{e,i}}{\partial p_{e,i}}=0$ (1)

With

${\gamma }_{e,i}={\left(1+M_{e,i}^2\left(p_{xe,i}^2+p_{ye,i}^2+p_{ze,i}^2\right)\right)}^{1/2}$ (2)

(the upper sign in Equation (1) is for electrons and the lower sign for ions, and subscripts e or i denote electrons or ions respectively). In our normalized units $M_e=1$ and $M_i=\frac{m_e}{m_i}$ .

We write the Hamiltonian of a particle in the electromagnetic field of the wave:

$H_{e,i}=\frac{1}{M_{e,i}}\left({\gamma }_{e,i}-1\right)\mp \phi \left(x,t\right)$ . (3)

where $\phi \left(x,t\right)$ is the scalar potential. Equation (1) can be reduced to a two-dimensional phase-space Vlasov equation as follows. The canonical momentum $P_{ce,i}$ is connected to the particle momentum $p_{e,i}$ by the relation $P_{ce,i}=p_{e,i}\mp a$ . $a=eA/m_{e}c$ is the normalized vector potential. From Equation (3), we can write:

$H_{e,i}=\frac{1}{M_{e,i}}\left({\left(1+M_{e,i}^2{\left(P_{ce,i}\pm a\right)}^2\right)}^{1/2}-1\right)\mp \phi \left(x,t\right)$ . (4)

Choosing the Coulomb gauge ($diva=0$ ), we have for our one dimensional problem $\frac{\partial a_x}{\partial x}=0$ , hence $a_x=0$ . The vector potential $a=a_{\perp }\left(x,t\right)$ , and we also have the following relation along the longitudinal direction:

$\frac{\text{d}{P}_{cxe,i}}{\text{d}t}=-\frac{\partial H_{e,i}}{\partial x}$ (5)

And since there is no transverse dependence

$\frac{\text{d}{P}_{c\perp e,i}}{\text{d}t}=-{\nabla }_{\perp }{H}_{e,i}=0$ . (6)

This last equation means $P_{c\perp e,i}=\text{const}$ . We can choose this constant to be zero without loss of generality, which means that initially all particles at a given (x, t) have the same perpendicular momentum $p_{\perp e,i}=\pm a_{\perp }\left(x,t\right)$ . The Hamiltonian is now written:

$H_{e,i}=\frac{1}{M_{e,i}}\left({\left(1+M_{e,i}^2{p}_{xe,i}^2+M_{e,i}^2{a}_{\perp }^2\left(x,t\right)\right)}^{1/2}-1\right)\mp \phi \left(x,t\right)$ . (7)

The 4D distribution function $F_{e,i}\left(x,p_{xe,i},p_{\perp e,i},t\right)$ can now be reduced to a 2D distribution function $f_{e,i}\left(x,p_{xe,i},t\right)$ :

$F_{e,i}\left(x,p_{xe,i},p_{\perp e,i},t\right)=f_{e,i}\left(x,p_{xe,i},t\right)\delta \left(p_{\perp e,i}\mp a_{\perp }\right)$ . (8)

$f_{e,i}\left(x,p_{xe,i},t\right)$ verify the relation:

$\frac{\text{d}{f}_{e,i}}{\text{d}t}=\frac{\partial f_{e,i}}{\partial t}+\frac{\partial H_{e,i}}{\partial p_{xe,i}}\frac{\partial f_{e,i}}{\partial x}-\frac{\partial H_{e,i}}{\partial x}\frac{\partial f_{e,i}}{\partial p_{xe,i}}=0$ . (9)

Which gives the following Vlasov equations for the electrons and the ions::

$\frac{\partial f_{e,i}}{\partial t}+M_{e,i}\frac{p_{xe,i}}{{\gamma }_{e,i}}\frac{\partial f_{e,i}}{\partial x}+\left(\mp E_x-\frac{M_{e,i}}{2{\gamma }_{e,i}}\frac{\partial a_{\perp }^2}{\partial x}\right)\frac{\partial f_{e,i}}{\partial p_{xe,i}}=0$ . (10)

where ${\gamma }_{e,i}={\left(1+{\left(M_{e,i}{p}_{xe,i}\right)}^2+{\left(M_{e,i}{a}_{\perp }\right)}^2\right)}^{1/2}$ .

$E_x=-\frac{\partial \phi }{\partial x}$ and $E_{\perp }=-\frac{\partial a_{\perp }}{\partial t}$ (11)

and Poisson’s equation is given by:

$\frac{{\partial }^2\phi }{\partial x^2}={\displaystyle \int f_e\left(x,p_{xe}\right)\text{d}{p}_{xe}}-{\displaystyle \int f_i\left(x,p_{xi}\right)\text{d}{p}_{xi}}$ (12)

An alternative way, used in the present work, is to calculate $E_x^{}$ from Ampère’s equation: $\frac{\partial E_x}{\partial t}=-J_x$ ,

The transverse electromagnetic fields $E_y,B_z$ and $E_z,B_y$ for the circularly polarized wave obey Maxwell’s equations. With $E^{\pm }=E_y\pm B_z$ and $F^{\pm }=E_z\pm B_y$ , we write these equations in the following form:

$\left(\frac{\partial }{\partial t}\pm \frac{\partial }{\partial x}\right)E^{\pm }=-J_y;\left(\frac{\partial }{\partial t}\mp \frac{\partial }{\partial x}\right)F^{\pm }=-J_z$ (13)

Which are integrated along their vacuum characteristic x = t [21]-[29]. In our normalized units we have the following expressions for the normal current densities:

$J_{\perp }=J_{\perp e}+J_{\perp i}\text{; }{J}_{\perp e,i}=-a_{\perp }{M}_{e,i}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{f_{e,i}}{{\gamma }_{e,i}}\text{d}{p}_{xe,i}}$ . (14)

2.2. The Numerical Scheme

We use a method of characteristics for the numerical solution of Equation (10), as discussed in [21]-[31], for the numerical solution of hyperbolic type partial differential equations These methods are Eulerian methods which use a computational mesh to discretize the equations on a fixed grid, and have been successfully applied to different important problems in plasma physics (see for instance the articles [21]-[31]).

The numerical scheme to advance Equation (10) from time $t_n$ to $t_{n+1}$ necessitates the knowledge of the electromagnetic field $E^{\pm }$ and $F^{\pm }$ at time $t_{n+1/2}$ . This is done using a centered scheme where we integrate Equation (13) exactly along the characteristics with $\Delta x=\Delta t$ , to calculate$E^{\pm n+1/2}$ and $F^{\pm n+1/2}$ as follows:

$E^{\pm }\left(x\pm \Delta t,t_{n+1/2}\right)=E^{\pm }\left(x,t_{n-1/2}\right)-\Delta t{J}_y\left(x\pm \Delta t/2,t_n\right)$ (15)

with $J_y\left(x\pm \Delta t/2,t_n\right)=\frac{J_y\left(x\pm \Delta x,t_n\right)+J_y\left(x,t_n\right)}{2}$

A similar equation can be written for $F^{\pm n+1/2}$ with J_z. From Equation (11), we also have $a_{\perp }^{n+1}=a_{\perp }^n-\Delta t{E}_{\perp }^{n+1/2}$ , from which we calculate $a_{\perp }^{n+1/2}=\left(a_{\perp }^{n+1}+a_{\perp }^n\right)/2$ . We also have:

$J_x^n=M_i{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{p_{xi}}{{\gamma }_i}{f}_i^n\text{d}{p}_{xi}}-M_e{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{p_{xe}}{{\gamma }_e}{f}_e^n\text{d}{p}_{xe}}$ (16)

We calculate $E_x^{n+1/2}$ from Ampère’s equation:$\frac{\partial E_x}{\partial t}=-J_x$ , from which $E_x^{n+1/2}=E_x^{n-1/2}-\Delta t{J}_x^n$ . Now given $f_{e,i}^n$ at mesh points (we stress here that the subscript i denotes the ion distribution function), we calculate the new value $f_{e,i}^{n+1}$ at mesh points from the relation [21]-[31]:

$f_{e,i}^{n+1}\left(x,p_{xe,i}\right)=f_{e,i}^n\left(x-2{\Delta }_{xe,i},p_{xe,i}-2{\Delta }_{p_{xe,i}}\right)$ . (17)

${\Delta }_{xe,i}$ and ${\Delta }_{p_{xe,i}}$ are calculated from the solution of the characteristics equations for Equation (10), which are given by:

$\begin{array}{l}\frac{\text{d}x}{\text{d}t}=M_{e,i}\frac{p_{xe,i}}{{\gamma }_{e,i}}\\ \frac{\text{d}{p}_{xe,i}}{\text{d}t}=\mp E_x-\frac{M_{e,i}}{2{\gamma }_{e,i}}\frac{\partial a_{\perp }^2}{\partial x}\end{array}$ (18)

This calculation is effected as follows. We rewrite Equation (17) in the vectorial form:

$f_{e,i}^{n+1}\left(X_{e,i}\right)=f_{e,i}^n\left(X_{e,i}-2{\Delta }_X{}_{e,i}\right)$ . (19)

$X_{e,i}$ is the two dimensional vector $X_{e,i}=\left(x,p_{xe,i}\right)$ , and ${\Delta }_X{}_{e,i}=\left({\Delta }_{xe,i},{\Delta }_{p_{xe,i}}\right)$ is the two dimensional vector calculated from the implicit relation [21]-[31]:

${\Delta }_X{}_{e,i}=\frac{\Delta t}{2}{V}_{e,i}\left(X_{e,i}-{\Delta }_X{}_{e,i},t_{n+1/2}\right)$ . (20)

$V_{e,i}=\left(M_{e,i}\frac{p_{xe,i}}{{\gamma }_{e,i}},\mp E_x^{n+1/2}-\frac{M_{e,i}}{2{\gamma }_{e,i}}\frac{\partial {\left(a_{\perp }^{\left(n+1/2\right)}\right)}^2}{\partial x}\right)$ . Equation (20) for ${\Delta }_X{}_{e,i}$ is implicit and is solved iteratively: ${\Delta }_{Xe,i}^{k+1}=\frac{\Delta t}{2}{V}_{e,i}\left(X_{e,i}-{\Delta }_{Xe,i}^k,t_{n+1/2}\right)$ , where we start the iteration with ${\Delta }_{Xe,i}^0=0$ for iteration number k = 0. Two or three iterations are sufficient to get a good convergence. Then $f_{e,i}^{n+1}$ is calculated from $f_{e,i}^n$ in Equation (17) by calculating the shifted value using a two-dimensional interpolation in the two dimensional phase-space $\left(x,p_{xe,i}\right)$ . Similarly in Equation (20), the shifted value was calculated at every iteration using a two-dimensional interpolation. In all cases, these two-dimensional interpolations are executed using a tensor product of cubic B-splines [21]-[23]. As mentioned, this code applies a numerical scheme based on a two-dimensional advection technique, of second order accuracy in time-step, where the value of the distribution function is advanced in time by interpolating in two dimensions along the characteristics using a tensor product of cubic B-splines. Cubic splines have been tested in several applications and shown to give accurate results [21]-[23]. See also [28] for the integration of the Vlasov equation in one and two dimensions.

2.3. A Circularly Polarized Laser Beam Normally Incident on an Overdense Plasma (n/n_cr = 25)

The forward propagating circularly polarized laser beam enters the system at the left boundary (x = 0), where the forward propagating fields are $E^{+}=2E_0{P}_r\left(t\right)\cos \left(\tau \right)$ and $F^{-}=-2E_0{P}_r\left(t\right)\sin \left(\tau \right)$ , as defined in Equation (13). Time is in units of ${\omega }^{-1}$ . In these expressions, $\tau =t-1.5t_p$ and $t_p=24$ is the pulse duration at full width at half maximum of the beam intensity. In our normalized units $E_0=a_0$ . For a Gaussian time dependence, the shape factor $P_r\left(t\right)$ is given by:

$P_r\left(t\right)=\exp \left(-2\ln \left(2\right){\left(\tau /t_p\right)}^2\right)$ (21)

so that the forward propagating laser beam behaves as a Gaussian pulse, which reaches its peak at $t=1.5t_p=36$ .(The units of time and space can be easily translated in units of ${\omega }_{pe}^{-1}$ and $c/{\omega }_{pe}$ by multiplying them by a factor of 5, since ${\omega }_{pe}/\omega =5$ , which corresponds to $n/n_{cr}=25$ , since the laser frequency $\omega /{\omega }_p=1/\sqrt{n/n_{cr}}=0.2$ ).

A characteristic parameter of high power laser beams is the normalized vector potential or quiver momentum $\left|a_{\perp }\right|=\left|e{A}_{\perp }/M_{e}c\right|=a_0$ , where $A_{\perp }$ is the vector potential of the wave. We chose for the amplitude of the vector potential. $a_0=25/\sqrt{2}$ For the circularly polarized wave $2a_0^2=I{\lambda }_0^2/1.368{\text{x10}}^{\text{18}}$ , I being the intensity in W/cm², and ${\lambda }_0$ the laser wavelength in microns. We use N = 5500 grid points in space, and in momentum space we use 1600 grid points for the electrons and 4800 grid points for the ions (extrema of the electron momentum are ±4, and ±360 for the ion momentum). The Lorentz factor for the transverse oscillation of an electron in the field of the wave is${\gamma }_0=\sqrt{1+a_0^2}=17.7$ . Deuterium ions are used with$M_i/M_e=2\times \text{1836}$ , where M_i is the deuterium mass. The initial distribution functions for electrons and ions are Maxwellians with temperatures T_e = 1 keV for the electrons and T_i = 0.1 keV for the deuterons. The length of the system is L = 22 in units of $c/{\omega }_{pe}$ . From the parameters presented we have a grid spacing $\Delta x=0.004$ , and we use a time step $\Delta t=\Delta x$ . We have a vacuum region of length L_vac = 9.0 on each side of the plasma slab. The steep ramp in density at the plasma edge on each side of the flat top density of the slab target has a length of $L_{edge}=0.3$ . The length of the central plasma slab, with flat top density of 1 (normalized to $25n_{cr}$ ), is $L_p=2.884$ , for a total length of 22. In our units, the skin depth $c/{\omega }_{pe}=\frac{c}{\omega }\frac{\omega }{{\omega }_{pe}}=0.2\frac{c}{\omega }$ . So the length of the large plasma slab is about $L_p\approx 15$ skin depth $c/{\omega }_{pe}$ . In free space $\omega =k$ for the electromagnetic wave, and as ${\lambda }_0=2\pi$ , the condition that ${\lambda }_0\gg L_{edge}$ is very well satisfied.

We first present in Figure 1 some of the results obtained for the density profiles against distance. The electron and ion density profiles at times t = 24, 40, 43, 50, 58, and 77 (time is in units of ${\omega }^{-1}$ ) are given in Figure 1 (full curve for the electrons, broken curves for the ions). The electric field (dashed-dotted curve) is also plotted divided by a factor of 5, to have it on the same figure. We concentrate these plots in the region where the plasma slab is located, since otherwise on both sides we have vacuum. In Figure 1 at t = 24 electrons are pushed and accelerated at the wave-front plasma-edge interface under the effect of the radiation pressure and form a steep density spike. Ions (dashed curve) have moved very little, and their profile gives essentially the initial plasma slab profile. An electric field is formed at the edge due to the charge separation appearing there, and its penetration in the plasma at the sharp gradient at the surface of the target is of the order of about $0.2c/\omega$ , which is the value of the skin depth in our units as we previously mentioned. In Figure 1 at t = 40 the steepening of the electrons is increased, due to the increase in the radiation pressure caused by the increase of the incident wave amplitude. The ion density profile is also increasing at the target surface.

Figure 1. Electron (full curves) and ion (dashed curves) density profiles. The electric field (dashed-dotted curves, divided by a factor of 5). Plots are at: t = 24, t = 40, t = 43, t = 50, t = 58, t = 77.

And then, very rapidly, there is a huge jump in the ion density at the target edge, as shown in Figure 1 at t = 43, reaching a peak of about 27. We note from Equation (21) that the laser intensity reaches a peak at 3t_p/2 = 36 at the boundary x = 0. This peak has to travel a distance of about x = 9 to reach the target surface (see for instance the position of the edge in Figure 1 at t = 24). So the maximum laser intensity will reach the target surface at about t ≈ 45, i.e. about the time where we observe a strong acceleration of the ions at the target surface in Figure 1 at t = 43. In Figure 1 at t = 50 we see the electron peak separating from the ion peak, which decreases compared to the peak in Figure 1 at t = 43. In Figure 1 at t = 58 the electron peak seems to come back closer to the ion peak, and a shock-like structure of an expanding neutral plasma jet propagating to the right is formed. In Figure 1 at t = 77 we see to the right two electron peaks coinciding with the two ion peaks, and the structure of a neutral plasma jet propagating in the forward direction. The solitary-like structure of the neutral plasma propagating to the right is now very close to the position of the back edge position of the initial target at about x = 13. In Figure 1 at t = 77 the electron and ion peaks to the left at the target surface are decreasing, due to the decay of the intensity of the Gaussian incident laser pulse, which is now exercising negligible pressure at the target surface at the left (see Figure 5 at t = 77). At this stage this wave-front plasma-edge interface is not propagating to the right anymore, since the radiation pressure has vanished. We have in Figure 1 at t = 77 a neutral plasma jet population at the right, with two peaks, which is now emerging from the back of what was the target and travelling in the forward direction to the right. And at the left, a standing wave-front plasma-edge interface.

Figure 2 shows the phase-space contour plots of the ion distribution function at t = 24, 43 (during the rapid acceleration of the ions at the target surface presented in Figure 1) at t = 43), 50, and 77. The velocity of the free streaming jet at the right can be calculated from Figure 2 at t = 77, for instance, where we see the edge of the free streaming deuterium ions in the forward direction having a normalized momentum around 350. Hence $M_i\upsilon /M_{e}c\approx 350$ , or $\upsilon /c\approx 350/\left(2\times \text{1836}\right)=0.095$ . This corresponds to the energy of the deuterium ions of $M_i{\upsilon }^2/2\approx M_i{c}^2{\left(\upsilon /c\right)}^2/2\approx 938\times 0.009=8.46\text{MeV}$ . The appearance of loops at the top in Figs. (2) at t = 50, t = 77 appear because slower ions (lower part of the loop), which were accelerated by the shock during the rise of the laser pulse, are caught up by faster ones (upper part of the loop), which were accelerated by the pulse when it was at its maximum. The rightmost parts of these loops correspond to the shock appearing at the right of the expanding neutral plasma edge, as we see for instance in Figure 1 at t = 77. These loops are close to what has been presented in [29] using PIC codes. See also the results in [30] [31].

Figure 3 presents the phase-space contour plots of the electron distribution function, at the same time as the corresponding phase-space contour plots of the ion distribution function presented in Figure 2. The steep edge at the target surface is due to the effect of the ponderomotive pressure pushing the electrons. In particular, in the initial phase, Figure 3 at t = 24 shows on the upper side of the contour small electron jets, a picture close to what has been reported in [32], and in [33] for the study of wave breaking and electron acceleration using a Vlasov code (a more detailed presentation of this phase of wave breaking is presented in Figure 6 below).

Figure 2. Phase-space plots of the ion distribution function. Plots at: t = 24, t = 43, t = 50, t = 77.

Note in Figure 3 at t = 43 and t = 50, beyond the solitary structure at the sharp front edge of the target, the long tail remaining from the initial electron distribution function, which is slowly disappearing. At t = 77 there is no ponderomotive pressure of the wave at the target surface (see Figure 5(c) below), the plasma at this surface location at x = 11.5 in Figure 1 at t = 77 is not moving anymore, and the electron and ion densities at this location are decaying slowly, with electrons slightly expanding in the backward direction as seen in Figure 1 at t = 77. But beyond the target surface, the neutral plasma jet is propagating to the right in the forward direction, at a constant free-streaming speed as discussed above. Note in the right of Figure 3 at t = 77 the two peaks of the neutral plasma jet we see at the right of Figure 1 at t = 77.

We follow in Figure 4 at t = 86 the two bump free-streaming plasma jet structure at the right in Figure 3 at t = 77. We concentrate the phase-space plots in Figure 4(a) on the ions, in Figure 4(b) on the electrons, and in Figure 4(c) we present the corresponding density and electric field plots.

We have indeed a stable solitary-like structure with a double peak, of a neutral

Figure 3. Phase-space plots of the electron distribution function. Plots at: t = 24, t = 43, t = 50, t = 77.

plasma free-streaming at constant speed to the right. The electron amd ion density curves (full curve and dashed curve respectively) are identical. An estimation of this speed can be obtained by following the right boundary of the jet. In Figure 1 at t = 77 the right boundary is located at x = 13.4. In Figure 4(c) the right boundary is at x = 14.3 at t = 86. This corresponds to a speed of 0.1, very close to the value of 0.095 estimated above from Figure 2 at t = 77.

Figure 5 presents in the left figures the forward propagating wave $E^{+}$ (full curves) and the backward reflected wave $E^{-}$ (dashed curves), and in the right figures presents the corresponding results for the forward propagating wave $F^{-}$ (full curves) and the backward propagating wave $F^{+}$ (dashed curves).

In Figure 5(a) the curves are calculated at t = 43 at the time of the rapid acceleration of the ions at the target surface observed in Figure 1 at t = 43, when the maximum of the wave at the target surface is very close to the peak of the pulse. From Equation (21), the laser intensity reaches a peak at 3t_p/2 = 36 at the boundary x = 0. This peak has to travel a distance of about x = 9 to reach the target surface (see Figure 1 at t = 24). So the maximum laser intensity will reach the target surface around t ≈ 45, i.e. about the time where we see the strong acceleration

Figure 4. Plots of the free-streaming plasma jet at t = 86: (a) the ion distribution function; (b) the electron distribution function; (c) plot of the ion density (dashed curve), of the electron density (full curve), and of the electric field (dashed-dotted curve, essentially zero in the free-streaming at constant velocity neutral plasma jet).

of the ions at the target surface in Figure 1 at t = 43. The electromagnetic wave damps in the plasma over a distance of the order of the skin depth $c/{\omega }_{pe}$ , equal to 0.2 in our normalized units. The strong increase of the ion and electron densities at the wave-front plasma-edge interface makes the plasma more opaque, with the steep plasma edge acting as a moving mirror for the incident light. Note that when $E^{+}$ is at a peak, the value of $F^{-}$ is zero and vice-versa, so the wave is always maintaining a pressure on the surface of the plasma., until the decay of the Gaussian pulse. The frequency of the backward reflected wave is slightly down-shifted by the moving reflecting plasma surface. We can check that in Figure 5, the reflected waves (dashed curves) in vacuum have a wavelength (and a period) slightly bigger than the corresponding ones for the incident wave (full curves). Figure 5(b) presents the incident and reflected waves at t = 60, at the time of the acceleration we observe in Figure 1 at t = 77, and in Figure 5(c) at t = 77 the Gaussian pulse is decaying and the ponderomotive pressure of the wave at the target surface is negligible. The neutral two-peaks structure at the right of Figure 1 at t = 77 is now free streaming as we previously mentioned, with a constant speed of 0.1.

Figure 5. Incident waves E⁺ and F⁻ (full curves) and reflected waves E⁻ and F⁺ (dashed curves) at (from top to bottom): a) t = 43, b) t = 60, c) t = 77.

We finally present in Figure 6 a sequence showing the wave breaking at the front edge of the electron phase-space, as presented in Figure 3 at t = 24. This is due to the longitudinal field generated at the front edge of the plasma due to the electron-ion separation. We see in the sequence of Figure 6 how peaks are building up at the target surface in the electron phase-space, then the accelerated electrons are ejected forming this finger-like structure we observe. The sequence is presented at t = 15.5, 17, and 20.

We have presented in Figure 6 the finger-like structures in the electron distribution function due to a wave-breaking mechanism, because these structures are well resolved by the fine (x, p_x) computational mesh of the Vlasov code, a picture close to what has been reported in [32], and in [33] for the study of wave breaking and electron acceleration using a Vlasov code. Our presentation of this phase of wave breaking we presented in Figure 6 contains more details because of the precision of the Vlasov code.

Figure 6. Phase-space plots of the electron distribution function at the target surface, showing the evolution of the wave braking process at the beginning, when the laser beam pulse is reaching the target surface.

The detailed analysis presented above in this section for the formation and ejection at constant speed of the neutral structure of a plasma jet with two peaks at t = 77, when the Gaussian electromagnetic pulse has damped, and applies no more pressure on the neutral population advancing at a constant speed as mentioned in the discussion above, testify to the power of the Eulerian codes for the numerical resolution of the details of this problem. Similar neutral plasma jets are obtained in several simulations where the ponderomotive pressure of the wave plays an important role in the formation and acceleration of the neutral plasma jet. See, for instance Figs. (1) and (5) in [7], and Fig. 1 in [34].

3. A Vlasov Code Simulation of the Amplification of Seed Pulses by Brillouin Backscattering in Magnetized Plasmas

Special attention has been given recently to the idea of plasma based laser amplifiers for their application to the direct amplification of ultra-short laser pulses. This problem has been studied with an Eulerian Vlasov code in an unmagnetized plasmas in [8], and in a magnetized plasma in [9] [10], to study the amplification of an ultra-short seed pulse via stimulated Raman or Brillouin backscattering of energy from a long pump pulse, assumed at constant amplitude, in a plasma embedded in an external magnetic field. In the present work, we apply an Eulerian Vlasov code to study the transfer of energy from a long pump electromagnetic wave, to a counter-propagating ultra-short seed pulse, mediated by a resonant ion wave via stimulated Brillouin backscattering (SBS), in a plasma embedded in an external magnetic field. Detailed analysis of the spectra developed during the amplification process is presented, together with the evolution showing the pump depletion, accompanied by the counter-propagating seed-pulse amplification, compression and increased steepness of the waveform. This problem is of importance for the amplification of seed pulses of whistler waves along the magnetic field of the Earth, and other similar planetary situations. The numerical code solves a one-dimensional relativistic Vlasov-Maxwell set of equations for a plasma in a magnetic field for both electrons and ions. The absence of noise in the Eulerian Vlasov code allows us to follow the evolution of the system with an accurate representation of the phase-space structures of the distribution functions.

The Relevant Equations of the Eulerian Vlasov Code and the Numerical Method

Time t is normalized to the inverse electron plasma frequency ${\omega }_{pe}^{-1}$ , velocity and momentum are normalized respectively to the velocity of light c and to M_ec, where M_e is the electron mass, length is normalized to $l_0=c{\omega }_{pe}^{-1}$ , and the electric field is normalized to $E_0=M_e{\omega }_{pe}c/e$ . The 1D relativistic Vlasov equations for the electrons and ions distribution functions $f_{e,i}\left(x,p_{xe,i},t\right)$ verify the relation [35]:

$\frac{\partial f_{e,i}}{\partial t}+\frac{p_{xe,i}}{m_{e,i}{\gamma }_{e,i}}\frac{\partial f_{e,i}}{\partial x}\mp \left(E_x+{\left(u_{\perp e,i}\times B_{\perp }\right)}_x\right)\cdot \frac{\partial f_{e,i}}{\partial p_{xe,i}}=0$ (22)

where ${\gamma }_{e,i}\approx {\left(1+{\left(p_{xe,i}/m_{e,i}\right)}^2\right)}^{1/2}$ . The upper sign in Equation (22) is for the electron equation and the lower sign for the ion equation, and subscripts e and i denote electrons and ions, respectively. In our normalized units $m_e=1$ , and $m_i=M_i/M_e=1836$ is the ratio of hydrogen ion to electron masses. We use a lower case letter for the normalized masses in this section, to differentiate from the upper case letter for the normalized masses for the deuterium used in Section 2. In this 1D model, the transverse velocity $u_{\perp e,i}\left(x,t\right)$ is assumed non-relativistic and is calculated from a transverse “fluid” equation, but the kinetic features of the plasma are fully preserved along the external magnetic field in the longitudinal direction x. The transverse momentum $p_{\perp e,i}\left(x,t\right)=m_{e,i}{u}_{\perp e,i}$ is obtained through the introduction of the generalized canonical momentum defined by $P_{\perp e,i}\left(x,t\right)=p_{\perp e,i}\left(x,t\right)\mp a_{\perp }\left(x,t\right)$ ($A$ is the vector potential of the wave, and $a=eA/M_{e}c$ is the normalized vector potential), which obeys the ‘fluid’ equation [35]:

$\left(\frac{\partial }{\partial t}+u_{xe,i}\frac{\partial }{\partial x}\right)P_{\perp e,i}=P_{\perp e,i}\times {\Omega }_{ce,i}\pm a_{\perp }\times {\Omega }_{ce,i}$ (23)

With $u_{xe,i}={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{p_{xe,i}}{m_{e,i}{\gamma }_{e,i}}{f}_{e,i}\text{d}{p}_{xe,i}}$ .

${\Omega }_{ce,i}={\Omega }_{ce,i}{e}_x$ , where ${\Omega }_{ce,i}$ is the electron (ion) cyclotron frequency, and $e_x$ is the unit vector in the direction of the external magnetic field.

We adopt the Coulomb gauge $\nabla \cdot a=0$ , where the vector potential is in the perpendicular plane. The longitudinal electric field $E_x$ is calculated either from the gradient of a potential $E_x=-\partial \phi /\partial x$ or from Ampére’s law:

$\frac{\partial E_x}{\partial t}=-J_x,J_x={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{p_{xi}}{m_i{\gamma }_i}{f}_i\text{d}{p}_{xi}}-{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{p_{xe}}{m_e{\gamma }_e}{f}_e\text{d}{p}_{xe}}$ (24)

The transverse electric field $E_{\perp }$ is calculated from the relations:

$E_{\perp }=-\frac{\partial a_{\perp }}{\partial t}$ (25)

The transverse electromagnetic fields $E_y,B_z$ and $E_z,B_y$ for the circularly polarized wave obey Maxwell’s equations. Defining $E^{\pm }=E_y\pm B_z$ and $F^{\pm }=E_z\pm B_y$ we have the following equations:

$\left(\frac{\partial }{\partial t}\pm \frac{\partial }{\partial x}\right)E^{\pm }=-J_y;\left(\frac{\partial }{\partial t}\mp \frac{\partial }{\partial x}\right)F^{\pm }=-J_z$ (26)

which are integrated along their vacuum characteristics $x=\pm t$ . In our normalized units we have the following expressions for the normal current densities:

$J_{\perp }=J_{\perp e}+J_{\perp i}\text{; }{J}_{\perp e,i}=\mp n_{e,i}\left(x,t\right)u_{\perp e,i}\left(x,t\right);n_{e,i}\left(x,t\right)={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{f}_{e,i}\text{d}{p}_{xe,i}}$ (27)

The numerical scheme applies a direct solution method of the Vlasov equation as a partial differential equation in phase-space. In Equation (22) ${\gamma }_{e,i}\approx {\left(1+{\left(p_{xe,i}/m_{e,i}\right)}^2\right)}^{1/2}$ is a function $p_x$ only, which allows to apply a time-splitting scheme with the separation of $x$ and $p_x$ variables by a method of fractional steps, as discussed for Eulerian Vlasov codes in several publications [21]-[23] [35]-[38].

Results

The normalized vector potential or quiver momentum for the laser beam is $a_{0p}$ . A forward propagating circularly polarized laser beam of constant amplitude enters the system at the left boundary, where for the forward pump $E^{+}=E_y+B_z=2E_{0p}\cos \left({\omega }_{0p}t\right)$ and $F^{-}=E_z-B_y=2E_{0p}sin\left({\omega }_{0p}t\right)$ at x = 0, ${\omega }_{0p}=3.2$ (normalized to ${\omega }_{pe}$ ), and $a_{0p}=0.02$ . In our normalized units $E_{0p}={\omega }_{0p}{a}_{0p}$ . We use a system of length L = 600. The number of grid points in space is N = 120000, so $\Delta x=\Delta t=0.005$ . The extrema of the momentum for the electrons are ±0.5, with 1000 grid points in velocity space, and for the ions these extrema are ±1.5 × 10⁻³ with 500 grid points in velocity space.

The pump precursor reaches the right boundary x = 600 at t = 600, since in our normalized units x = t. The seed pulse is injected at the right boundary x = 600 towards the backward direction just before the arrival of the pump, in the time interval 540 < t < 600, with$E^{-}=E_y-B_z=2E_{0s}{P}_r\left(t\right)\cos \left({\omega }_{0s}\left(t-540\right)\right)$ , $F^{+}=E_z+B_y=-2E_{0s}{P}_r\left(t\right)sin\left({\omega }_{0s}\left(t-540\right)\right)$ . The SBS produces a very small frequency downshift in the scattered wave ${\omega }_{0p}={\omega }_{0s}+{\omega }_{ion}$ . We neglect this small downshift ${\omega }_{ion}\ll 1$ , and set the seed frequency ${\omega }_{0s}\approx {\omega }_{0p}=3.2$ , therefore $k_{0p}\approx k_{0s}$ . The SBS is backward propagating, so the selection rule $k_{0p}=-k_{0s}+k_{ion}$ will result in an ion wave $k_{ion}\approx 2k_{0p}$ . Also $a_{0p}=a_{0s}=0.02$ and $E_{0s}=E_{0p}$ . The initial distribution functions for electrons and ions are Maxweellian with temperatures T_e = 0.3 keV and T_i = 0.02 keV. respectively. The shape factor P_r(t) for the injected seed pulse has a Gaussian time dependence $P_r\left(t\right)=\exp {\left(-0.5\left(t-570\right)/t_p\right)}^2$ , for 540 < t < 600, with t_p = 10. The backward propagating seed pulse reaches its peak at t = 570, and has completely penetrated the plasma at t = 600, i.e. the time the pump precursor reaches the right boundary. The electron cyclotron frequency (normalized to the electron plasma frequency) is ${\omega }_{ce}=0.5$ . From the linear dispersion relation $k^2={\omega }^2-1/\left(1\mp {\omega }_{ce}/\omega \right)$ the corresponding wavenumber for the left hand circularly polarized (LCP) wave is 3.062, and for the right hand circularly polarized wave (RCP) the wavenumber is 3.01.

Figure 7 shows at t = 1200 successively: in the first frame the incident pump wave $E^{+}$ (full curve) and the growing backward seed pulse $E^{-}$ (dashed curve), in the second frame is the longitudinal electric field $E_x$ , the third and fourth frames show the electron and ion density.

Figure 7. From top to bottom, first row: the incident pump wave $E^{+}$ (full curve). The growing backward seed pulse $E^{-}$ (dashed curve), then the electric field $E_x$ ; second row, the electron density and the ion density respectively.

Figure 8 presents contour plots for the electron distribution function, and Figure 9 for the ion distribution function, at t = 1200. In Figure 10, the wavenumber spectra associated with the forward wave $E^{+}$ , the backward wave $E^{-}$ , the longitudinal electric field $E_x$ and the ion density, are calculated by Fourier transform of the corresponding signal in the domain x = (75,403) at t = 1200, at a time where the backward seed pulse has reached the left boundary. The spectrum of $E^{+}$ in Figure 10) shows the peak at $k_{0p}=3.049$ (the linear dispersion relation calculates a wavenumber of 3.062 for the LCP wave for a frequency of 3.2 as mentioned above). The same mode is appearing in the wavenumber spectrum of the backward wave $E^{-}$ .

Figure 8. Contour plots for the electron distribution function.

Figure 9. Contour plots for the ion distribution function.

The presence of these two modes will result in the appearance of a SBS mode at $k_{ion}\approx 2k_{0p}=6.089$ , very close to the peak at 6.12 in the wavenumber spectrum of the ion density in Figure 10. Note the harmonics of 6.12 at 12.25, 18.37 and 24.5. We note in the wavenumber spectrum of $E^{+}$ the presence of a forward scattered mode at 2.01, with the selection rule 3.049 = 2.01 + k, from which k = 1.039, very close to the peak at 1.05 in the spectrum of $E_x$ . The corresponding frequencies we see in Figure 11) by taking the temporal Fourier transform of the signals registered at x = 90 in the time t = (1070, 1234) give a frequency of 3.2 for the pump $E^{+}$ and 2.186 for the forward scattered mode (the corresponding wavenumber calculated from the linear dispersion relation is 1.991, very close to the value of 2.01 we see in the wavenumber spectrum of $E^{+}$ ). They obey the selection rule $3.2=2.186+\omega$ , from which a plasma mode at $\omega =1.014$ , very close to the peak at 1.035 we see in the frequency spectrum of $E_x$ in Figure 11. Also in the wavenumber spectrum of $E^{-}$ a backward scattered mode at $k_s=2.01$ is dominating the spectrum, such that $k_{0p}=-k_s+k_e$ , from which $k_e=3.049+2.01=5.059$ , very close to the value of 5.043 we see in the wavenumber spectrum of $E_x$ . Note the harmonics of 5.043 at 10.086, 15.13 and 20.13. The forward scattered mode at 2.01 and the backscattered mode at −2.01 will couple through Brillouin scattering to give another SBS mode at $k_{ion}\approx 2\times 2.01=4.02$ , very close to the value of 4.0 we see in the ion spectrum in Figure 10). So there are two SBS events, at 4.0 and 6.12. Other resonances are also

Figure 10. Wavenumber spectra in x = (75,403) at t = 1200.

Figure 11. Temporal Fourier transform in t = (1070, 1234) at x = 90.

present. We see in the spectrum of $E^{+}$ the anti-Stokes at 4.103 obtained by the coupling of the plasma mode at 1.05 and the pump $E^{+}$ at the wavenumber 3.049, such that 3.049 + 1.05 = 4.099, very close to the value of 4.103 in the wavenumber spectrum of $E^{+}$ , with the corresponding frequencies 3.2 + 1.035 = 4.235, very close to the value of 4.29 we see in the frequency spectrum of $E^{+}$ in Figure 11).

Again, accurate results have been presented in Section 3 for the problem of the amplification of seed pulses by Brillouin backscattering in magnetized plasmas. For more information, see [8]-[10]. The relevance to the amplification of whistler pulses along the Earth’s magnetic field, and along the magnetic field in other planetary and space plasmas situations, is an example of how space plasma can be used as a laboratory to test the plasma theory.

4. Generation of Electromagnetic Radiation in a Relativistic Beam-Plasma Interaction

Electromagnetic radiation around the plasma frequency generated in beam-plasma interaction are frequent in observations of solar radio bursts when relativistic bunches of electrons are ejected from the sun [39] [40], generally known as Type III solar radio bursts, as well as in many tokamak experiments in the presence of runaway electrons [41]. These emissions are usually thought to be generated by a strong beam-plasma interaction, and the waves are subsequently converted by wave mixing to modes which can be radiated. It has been pointed out however that in the case of Type III solar radio emissions, the electromagnetic emission at the harmonic at twice the plasma frequency can be frequently stronger or at the same level as the emission at the fundamental plasma frequency [42]. We use a relativistic Vlasov-Maxwell code similar to what has been presented in Sections 2 and 3, which has been presented in [36] [37], to study this problem in one spatial dimension, and the results have been presented in details in [11]. The relativistic electron beam interacts with a plasma, and a plasma wave is excited, which interacts through a wave-mixing process with an electromagnetiv wave at the plasma frequency to excite an electromagnetic wave at twice the plasma frequency. This is essentially a stimulated Raman scattering process in a bump on tail instability, in a relativistic beam-plasma interaction, During the evolution of this system, the total electromagnetic action is conserved. The scattered mode at twice the plasma frequency would grow to a level equal to the decaying low frequency pump. Increasing the beam density further increases the beam-plasma instability and the level of saturation of the scattered wave (see Fig. 1 in [11]). Full details are discussed in [11]. For such an experiment with a relativistic electron beam, space plasma provides the laboratory to verify the theory.

5. The Generation of Electromagnetic Radiation by a Magnetic Field Aligned Electron Beam; the Auroral Kilometric Radiation

We discuss the problem of the auroral kilometric radiation [12] [13] [43], the colourful dance of lights in the sky in the Northern and Southern hemispheres of the Earth, when weakly relativistic electron beams are precipitating along the magnetic field in the north and south poles. This radiation is emitted inside cold-plasma depleted regions, with the plasma frequency normalized to the emission frequency ${\omega }_p/\omega =0.1$ much less than 1, and the emission frequency is close to the extraordinary mode cutoff, that is very close to the electron cyclotron frequency. Under this condition, the two branches of the extraordinary mode have an eigenfrequency very close to the electron cyclotron frequency (see Fig. 2 in [12]), and that a weakly relativistic beam (velocity $\upsilon /c$ of the order of 0.2) propagating along the magnetic field lines will couple these two branches at a frequency very close to the electron cyclotron frequency. It is suggested that this mechanism can be the source for the auroral kilometric radiation. Full details are presented in [12]. Analysis using the effect of a beam thermal spread, a horseshoe distribution function close to the experimental observation, on this instability has been investigated in [13] [43]. The theory together with the observations of the events are providing good information on these events, It is to be conjectured that the same process would apply to other planetary radio sources such as the decametric Jovian radiation [44], and the kilometric emission from Saturn [45], as well as to other astrophysical situations such as microwave spikes observed in many solar flares [46], and electromagnetic radiation associated with pulsars [47].

6. Conclusions

We have applied Eulerian codes for the numerical solution of the one-dimensional relativistic Vlasov-Maxwell set of equations, for problems pertinent to laboratory and space plasmas. We have presented in Section 2 the case of a circular polarization laser beam normally incident on an overdense plasma target. We have followed the effect of a non-oscillatory radiation pressure applied at the target surface of the plasma by the incident circularly polarized laser wave. This radiation pushes the electrons at the plasma surface through the ponderomotive pressure, producing a sharp density gradient at the target plasma surface. There is a build-up of the electron density at this surface, which creates a large space-charge, giving rise to a strong longitudinal electric field. This electric field at the wave-front plasma-edge interface accelerates the ions in the forward direction. In the ongoing acceleration process, the accelerated ions reach a free-streaming expansion phase where the electrons neutralize the charge of the expanding ions, and a neutral plasma jet is formed. When the Gaussian pulse of the wave finally damps, a dense compact bunch of neutral plasma jet is created and is free-streaming in the forward direction, as we observe in Figs. (1) and (3) at t = 77, and in Figure 4 at t = 86. Comparison presented in [29], of this Eulerian Vlasov code results, with the results of a particle-in-cell PIC code, to study a similar problem, has shown good agreement between the macroscopic quantities calculated with both codes. We have complemented the work we present here with many references to related problems. See also the interesting article in [48]. It is believed this mechanism provides the main acceleration mechanism for cosmic rays at the front of collisionless shock waves generated during supernova explosions by the interaction of high intensity electromagnetic waves with plasmas [6].

In Section 3, we have applied an Eulerian Vlasov code to study the transfer of energy from a long pump electromagnetic wave to a counter-propagating ultra-short seed pulse, mediated by a resonant ion wave via stimulated Brillouin backscattering (SBS), in a plasma embedded in an external magnetic field. Detailed analysis of the spectra developed during the amplification process is presented, together with the evolution showing the pump depletion, accompanied by the counter-propagating seed-pulse amplification, compression and increased steepness of the waveform. This problem is of importance for the amplification of seed pulses of whistler waves along the magnetic field of the Earth, and along the magnetic field in other planetary and space plasmas situations. Section 4 presented the application of a Vlasov code for the simulation of electromagnetic radiations around the plasma frequency generated in beam-plasma interaction. These are frequent in observations of solar radio bursts of Type III, when relativistic bunches of electrons are ejected from the sun [39] [40] [42]. The results we discuss here have been presented in [11]. The relativistic electron beam interacts with a plasma, and a plasma wave is excited, which interacts through a wave-mixing process with an electromagnetic wave at the plasma frequency to excite an electromagnetic wave at twice the plasma frequency. This is essentially a stimulated Raman scattering process in a bump-on-tail instability, in a relativistic beam-plasma interaction. Another example is space plasma, providing a laboratory to test the theory.

We have discussed in Section 5 the problem of the auroral kilometric radiation [12] [13] [43], when weakly relativistic electron beams are precipitating along the magnetic field at the north and south poles. This radiation is emitted inside cold-plasma depleted regions, with the plasma frequency normalized to the emission frequency ${\omega }_p/\omega =0.1$ much less than 1, and the emission frequency is close to the extraordinary mode cutoff, which is very close to the electron cyclotron frequency. Under these conditions, the two branches of the extraordinary mode have an eigenfrequency very close to the electron cyclotron frequency (see Fig. 2 in [12]), and a weakly relativistic beam (velocity $\upsilon /c$ of the order of 0.2) propagating along the magnetic field lines can couple these two branches. It is suggested that this mechanism is the source of the auroral kilometric radiation. Analysis using the effect of a beam thermal spread on this instability has been investigated in [13] [43]. The theory, together with the observations of the events, provides good information on these events. It is to be conjectured that the same process would apply to other planetary radio sources, such as the decametric Jovian radiation [44], and the kilometric emission from Saturn [45], as well as to other astrophysical situations, such as microwave spikes observed in many solar flares [46], and electromagnetic radiation associated with pulsars [47].

There has been important development in the direct numerical solution of the relativistic Vlasov equation as a partial differential equation in phase-space [21]-[28] [35]-[38]. Interest in Eulerian grid-based solvers associated with the method of characteristics for the numerical solution of the Vlasov equation arises from the absence of numerical noise associated with these numerical methods, which allows the accurate study of the low density regions of the phase-space where particles are accelerated. We have presented a few examples pertinent to laboratory and space plasmas. It is worth mentioning some of the early pioneering works that applied Eulerian codes for the numerical solution of the Vlasov equation on a fixed grid [49]-[59], who have shown how to split the multidimensional Vlasov equation problem into a succession of 1D advections. This greatly simplified the task of obtaining numerical solution of the Vlasov equation while maintaining second order accuracy in time.

NOTES

*Retired senior scientist, Hydro-Quebec Research Institute.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Shoucri, M. (2010) Eulerian Codes for the Numerical Solution of the Kinetic Equations of Plasmas. Nova Science Publishers.
[2]	Klimo, O., Psikal, J., Limpouch, J. and Tikhonchuk, V.T. (2008) Monoenergetic Ion Beams from Ultrathin Foils Irradiated by Ultrahigh-Contrast Circularly Polarized Laser Pulses. Physical Review Special Topics—Accelerators and Beams, 11, Article ID: 031301.[CrossRef]
[3]	Kar, S., Borghesi, M., Bulanov, S.V., Key, M.H., Liseykina, T.V., Macchi, A., et al. (2008) Plasma Jets Driven by Ultraintense-Laser Interaction with Thin Foils. Physical Review Letters, 100, Article ID: 225004.[CrossRef] [PubMed]
[4]	Romagnani, L., Bulanov, S.V., Borghesi, M., Audebert, P., Gauthier, J.C., Löwenbrück, K., et al. (2008) Observation of Collisionless Shocks in Laser-Plasma Experiments. Physical Review Letters, 101, Article ID: 025004.[CrossRef] [PubMed]
[5]	Fernández, J.C., Honrubia, J.J., Albright, B.J., Flippo, K.A., Gautier, D.C., Hegelich, B.M., et al. (2009) Progress and Prospects of Ion-Driven Fast Ignition. Nuclear Fusion, 49, Article ID: 065004.[CrossRef]
[6]	Berezinskii, V.S., Bulanov, V.S., Dogiel, V.A., Ginzburg, V.L. and Ptuskin, V.S. (1990) Astrophysics of Cosmic Rays. Elsevier.
[7]	Shoucri, M., Vidal, F. and Matte, J. (2016) Formation of Double Layers and Evolution of the Distribution Functions during Ion Acceleration Driven by a High-Intensity Short Laser Pulse Normally Incident on Thin Foils. Laser and Particle Beams, 34, 714-728.[CrossRef]
[8]	Shoucri, M., Matte, J. and Vidal, F. (2015) Relativistic Eulerian Vlasov Simulations of the Amplification of Seed Pulses by Brillouin Backscattering in Plasmas. Physics of Plasmas, 22, Article ID: 053101.[CrossRef]
[9]	Shoucri, M. (2016) Numerical Simulation of Raman and Brillouin Laser-Pulse Amplification in a Magnetized Plasma. Laser and Particle Beams, 34, 315-337.[CrossRef]
[10]	Shoucri, M. (2016) A Vlasov Code Simulation of the Amplification of Seed Pulses by Brillouin Backscaring in Magnetized Plasmas. Proceedings of the 43rd European Physical Society Conference on Plasma Physics, Leuven, 4-8 July 2016, 1.068.
[11]	Shoucri, M., Bertrand, P., Ghizzo, A., Lebas, J., Johnston, T., Feix, M., et al. (1991) Generation of Electromagnetic Radiation in a Relativistic Beam-Plasma Interaction. Physics Letters A, 156, 76-80.[CrossRef]
[12]	Shoucri, M. (1985) Generation of Electromagnetic Radiation by a Field-Aligned Electron Beam. The Physics of Fluids, 28, 666-668.[CrossRef]
[13]	Bingham, R. and Cairns, R.A. (2000) Generation of Auroral Kilometric Radiation by Electron Horseshoe Distributions. Physics of Plasmas, 7, 3089-3092.[CrossRef]
[14]	Mourou, G.A., Tajima, T. and Bulanov, S.V. (2006) Optics in the Relativistic Regime. Reviews of Modern Physics, 78, 309-371.[CrossRef]
[15]	Borghesi, M., Kar, S., Romagnani, L., Toncian, T., Antici, P., Audebert, P., et al. (2007) Impulsive Electric Fields Driven by High-Intensity Laser Matter Interactions. Laser and Particle Beams, 25, 161-167.[CrossRef]
[16]	Robson, L., Simpson, P.T., Clarke, R.J., Ledingham, K.W.D., Lindau, F., Lundh, O., et al. (2006) Scaling of Proton Acceleration Driven by Petawatt-Laser-Plasma Interactions. Nature Physics, 3, 58-62.[CrossRef]
[17]	Liseykina, T.V., Borghesi, M., Macchi, A. and Tuveri, S. (2008) Radiation pressure acceleration by ultraintense laser pulses. Plasma Physics and Controlled Fusion, 50, Article ID: 124033.[CrossRef]
[18]	Schlegel, T., Naumova, N., Tikhonchuk, V.T., Labaune, C., Sokolov, I.V. and Mourou, G. (2009) Relativistic Laser Piston Model: Ponderomotive Ion Acceleration in Dense Plasmas Using Ultraintense Laser Pulses. Physics of Plasmas, 16, Article ID: 083103.[CrossRef]
[19]	Robinson, A.P.L., Gibbon, P., Zepf, M., Kar, S., Evans, R.G. and Bellei, C. (2009) Relativistically Correct Hole-Boring and Ion Acceleration by Circularly Polarized Laser Pulses. Plasma Physics and Controlled Fusion, 51, Article ID: 024004.[CrossRef]
[20]	Shoucri, M. and Afeyan, B. (2009) Harmonics Generation in the Reflection of a Linearly Polarized Laser Beam Normally Incident on an Overdense Plasmas. Proceedings of the 36th European Physical Society Conference on Plasma Physics, Sofia, 29 June-3 July 2009, P1.022.
[21]	Shoucri, M. (2008) Numerical Solution of Hyperbolic Differential Equations. Nova Science Publishers.
[22]	Shoucri, M. (2008) Eulerian Codes for the Numerical Solution of the Vlasov Equation. Communications in Nonlinear Science and Numerical Simulation, 13, 174-182.[CrossRef]
[23]	Shoucri, M. (2008) The Application of the Method of Characteristics for the Numerical Solution of Hyperbolic Differential Equations. In: Colombo, S.P., et al., Eds., Numerical Simulation Research Progress, Nova Science Publishers, 1-98.
[24]	Ruhl, H. and Mulser, P. (1995) Relativistic Vlasov Simulation of Intense Fs Laser Pulse-Matter Interaction. Physics Letters A, 205, 388-392.[CrossRef]
[25]	Califano, F. and Mangeney, A. (2010) A Vlasov Approach to Collisionless Space and Laboratory Plasmas. In: Shoucri, M., Ed., Eulerian Codes for the Numerical Solution of the Kinetic Equations of Plasmas, Nova Science Publishers, 23-63.
[26]	Mangeney, A., Califano, F., Cavazzoni, C. and Travnicek, P. (2002) A Numerical Scheme for the Integration of the Vlasov-Maxwell System of Equations. Journal of Computational Physics, 179, 495-538.[CrossRef]
[27]	Macchi, A., Cattani, F., Liseykina, T.V. and Cornolti, F. (2005) Laser Acceleration of Ion Bunches at the Front Surface of Overdense Plasmas. Physical Review Letters, 94, Article ID: 165003.[CrossRef] [PubMed]
[28]	Shoucri, M., Gerhauser, H. and Finken, K. (2003) Integration of the Vlasov Equation along Characteristics in One and Two Dimensions. Computer Physics Communications, 154, 65-75.[CrossRef]
[29]	Shoucri, M., Lavocat-Dubuis, X., Matte, J. and Vidal, F. (2011) Numerical Study of Ion Acceleration and Plasma Jet Formation in the Interaction of an Intense Laser Beam Normally Incident on an Overdense Plasma. Laser and Particle Beams, 29, 315-332.[CrossRef]
[30]	Shoucri, M. (2010) Numerical Solution of the Relativistic Vlasov-Maxwell Equations for the Study of the Interaction of a High Intensity Laser Beam Normally Incident on an Overdense Plasma. In: Shoucri, M., Ed., Eulerian Codes for the Numerical Solution of the Kinetic Equations of Plasmas, Nova Science Publishers, 163-236.
[31]	Shoucri, M. (2012) Ion Acceleration and Plasma Jet Formation in the Interaction of an Intense Laser Beam Normally Incident on an Overdense Plasma: A Vlasov Code Simulation. Computational Science & Discovery, 5, Article ID: 014005.[CrossRef]
[32]	Bergmann, A. and Mulser, P. (1993) Breaking of Resonantly Excited Electron Plasma Waves. Physical Review E, 47, 3585-3589.[CrossRef] [PubMed]
[33]	Sircombe, N.J., Arber, T.D. and Dendy, R.O. (2004) Accelerated Electron Populations Formed by Langmuir Wave-Caviton Interactions. Physics of Plasmas, 12, Article ID: 012303.[CrossRef]
[34]	Shoucri, M., Vidal, F. and Matte, J. (2016) A Numerical Study of Ponderomotive Ion Acceleration in a Dense Plasma Driven by a Circularly Polarized High-Intensity Laser Beam Normally Incident on Thin Foils. Laser and Particle Beams, 34, 242-262.[CrossRef]
[35]	Bertrand, P., Ghizzo, A., Karttunen, S.J., Pättikangas, T.J.H., Salomaa, R.R.E. and Shoucri, M. (1992) Simulations of Wave-Particle Interactions in Stimulated Raman Forward Scattering in a Magnetized Plasma. Physics of Fluids B: Plasma Physics, 4, 3590-3607.[CrossRef]
[36]	Ghizzo, A., Bertrand, P., Shoucri, M.M., Johnston, T.W., Fualkow, E. and Feix, M.R. (1990) A Vlasov Code for the Numerical Simulation of Stimulated Raman Scattering. Journal of Computational Physics, 90, 431-457.[CrossRef]
[37]	Bertrand, P., Ghizzo, A., Johnston, T.W., Shoucri, M., Fijalkow, E. and Feix, M.R. (1990) A Nonperiodic Euler-Vlasov Code for the Numerical Simulation of Laser–plasma Beat Wave Acceleration and Raman Scattering. Physics of Fluids B: Plasma Physics, 2, 1028-1037.[CrossRef]
[38]	Shoucri, M. and Afeyan, B. (2016) Vlasov-Maxwell Simulations of Backward Raman Amplification of Seed Pulses in Plasmas. Laser and Particle Beams, 34, 576-600.[CrossRef]
[39]	Schrijver, H. (1973) The Time Dependence of Turbulent Resistivity in a Linear Plasma Column. Physica, 70, 339-357.[CrossRef]
[40]	Hutchinson, I.H., Molvig, K. and Yuen, S.Y. (1978) Mechanism for the ω Peradiation in Tokamaks. Physical Review Letters, 40, 1091-1094.[CrossRef]
[41]	Dendy, R.O., Lashmore-Davies, C.N. and Shoucri, M.M. (1985) A Triple-Wave Resonance Model for the Emission from Tokamaks of Narrow-Band Radiation at the Plasma Frequency. Nuclear Fusion, 25, 721-726.[CrossRef]
[42]	Gurnett, D.A. and Frank, L.A. (1975) Electron Plasma Oscillations Associated with Type III Radio Emissions and Solar Electrons. Solar Physics, 45, 477-493.[CrossRef]
[43]	Speirs, D.C., Bingham, R., Cairns, R.A., Vorgul, I., Kellett, B.J., Phelps, A.D.R., et al. (2014) Backward Wave Cyclotron-Maser Emission in the Auroral Magnetosphere. Physical Review Letters, 113, Article ID: 155002.[CrossRef] [PubMed]
[44]	Burke, B.F. and Franklin, K.L. (1955) Observations of a Variable Radio Source Associated with the Planet Jupiter. Journal of Geophysical Research, 60, 213-217.[CrossRef]
[45]	Kaiser, M.L., Desch, M.D., Warwick, J.W. and Pearce, J.B. (1980) Voyager Detection of Nonthermal Radio Emission from Saturn. Science, 209, 1238-1240.[CrossRef] [PubMed]
[46]	Slottje, C. (1978) Millisecond Microwave Spikes in a Solar Flare. Nature, 275, 520-521.[CrossRef]
[47]	Mikhailovskii, A.B. (1980) Oscillations of an Isotropic Relativistic Plasma. Plasma Physics, 22, 133-149.[CrossRef]
[48]	Schmidt, P., Boine-Frankenheim, O. and Mulser, P. (2015) Optimum Laser Parameters for 1D Radiation Pressure Acceleration. Laser and Particle Beams, 33, 387-396.[CrossRef]
[49]	Cheng, C.Z. and Knorr, G. (1976) The Integration of the Vlasov Equation in Configuration Space. Journal of Computational Physics, 22, 330-351.[CrossRef]
[50]	Cheng, C.Z. (2010) Splitting Methods for the Vlasov-Maxwell Equations in Plasma. In: Shoucri, M., Ed., Eulerian Codes for the Numerical Solution of the Kinetic Equations of Plasmas, Nova Science Publishers, 1-21.
[51]	Gagné, R.R.J. and Shoucri, M.M. (1977) A Splitting Scheme for the Numerical Solution of a One-Dimensional Vlasov Equation. Journal of Computational Physics, 24, 445-449.[CrossRef]
[52]	Shoucri, M.M. and Gagné, R.R.J. (1977) Numerical Solution of a Two-Dimensional Vlasov Equation. Journal of Computational Physics, 25, 94-103.[CrossRef]
[53]	Shoucri, M.M. and Gagné, R.R.J. (1978) Splitting Schemes for the Numerical Solution of a Two-Dimensional Vlasov Equation. Journal of Computational Physics, 27, 315-322.[CrossRef]
[54]	Shoucri, M.M. (1978) Computer Simulation of the Sideband Instability. The Physics of Fluids, 21, 1359-1365.[CrossRef]
[55]	Shoucri, M.M. (1979) Numerical Solution of the Two-Dimensional Vlasov Equation. IEEE Transactions on Plasma Science, 7, 69-72.[CrossRef]
[56]	Simon, A. and Short, R.W. (1988) Comments on ‘‘Motion of an Electron Bunch through a Plasma’’ [Phys. Fluids 2 9, 262 (1986)]. The Physics of Fluids, 31, 217-217.[CrossRef]
[57]	Ghizzo, A., Shoucri, M.M., Bertrand, P., Feix, M. and Fijalkow, E. (1988) Nonlinear Evolution of the Beam-Plasma Instabilities. Physics Letters A, 129, 453-458.[CrossRef]
[58]	Strozzi, D.J., Shoucri, M.M., Bers, A., Williams, E.A. and Langdon, A.B. (2006) Vlasov Simulations of Trapping and Inhomogeneity in Raman Scattering. Journal of Plasma Physics, 72, 1299-1302.[CrossRef]
[59]	Sonnendrücker, E., Roche, J., Bertrand, P. and Ghizzo, A. (1999) The Semi-Lagrangian Method for the Numerical Resolution of the Vlasov Equation. Journal of Computational Physics, 149, 201-220.[CrossRef]