Particle-in-cell Simulations of Raman Forward Scattering Instability in Low-density Plasmas: A Computational Study

Electromagnetic-wave propagation through equilibrium p lasmas can lead to medium instabilities whose formation and growth conditions are of importance to study. Raman forward scattering (RFS) instability in low-density plasmas is investigated through particle-in-cell (PIC) simulation, which is employed due to plasmas’ high complexity, costly experiment process, and nonlinear behavior in physical terms. A 2D electromagnetic PIC code is used to model RFS instability. Doing so, a plasma medium containing ions and electrons is PIC simulated and ions are assumed to form a constant background because of high inertia . The results indicate that as a plane electromagnetic wave enters the plasma medium, large-amplitude longitudinal waves, whose growth rate is well consistent with the theoretical results, begin to grow, suggesting the appearance of RFS instability.


Introduction
Various parameters such as instabilities are significant to describe equilibriu m plasmas. A system is in instable equilibriu m whenever there is a gro wing perturbation affecting the system's total configuration or macroscopic quantities. In a simp le classificat ion, the instabilities are divided into two main categories of position-space instabilities (macroscopic) and velocity-space instabilities (microscopic), wh ich include electrostatic and electro magnetic instabilities. A plas ma current creates a magnetic field shrinking plasma and increasing current density, wh ich causes electromagnetic instability. In this instance, 0 E ∇ × ≠   and the set of Maxwell equations must be comp letely solved [1].
Here, we d ep loy co mp uter simu lation to inv estigate p las ma in s tab ilit ies . Th is met ho d is class ified as a co mp utation al s cien ce in wh ich th e start po int is to s cientifically s tud y th e math ematical mo d el o f th e phenomenon. The equations of the mathematical model must become algebraically discrete in order to be understandable in terms of numerical solution. When the discrete equations are exp ressed as a series o f co mputer co mmands, they describe the simu lat ion model as a co mputer s imu lat ion program. A mass of data is produced even by the simp lest simu lation calculations. Hence, most research efforts in this field have been focused to obtain appropriate simu lation methods for physical systems that can be investigated using available co mputer resources. The main difference between a simu lated plasma and a real plasma is in the number of charges, fields, and spatial/temporal criteria. A real plas ma is composed of a mu ltitude of positive electrons and ions. In a simu lation process, a charged particle is essentially a homogeneous ensemble of numerous charged particles existing in a real plas ma. Therefore, its charge and mass is several times larger than real particles'. However, the simu lated particles' charge/mass ratio is the same as that of real part icles; consequently, we can substitute a lower number of charges in simu lation and simu late physical phenomena by a limited nu mber of particles. One of the main advantages of simu lation is the low volu me of arrays and calculations. The application of the word "particle" in computer memo ry represents the concept of charged particles occupying the real space.
Note that, many theoretical investigations regarding Raman instability have been well carried out and this paper aims to reconcile the PIC simu lation results with the theoretical results.

RFS Instability and the Used Equations
In Raman instability, an input large-amplitude light wave is transformed into a scattered light wave and a plas ma electron wave. The phase velocity in direct forward Instability in Low-density Plasmas: A Computational Study scattering is close to the speed of light and there are a tiny number of particles in the plasma background, absorbing energy required to be captured by these waves and accelerated consequently. If the input wave is propagated through a low-density plasma in a relat ively long distance, a large-amp litude plas ma electron wave will be created that can produce high-energy electrons [2].
An input wave ( A wave must penetrate into the plasma in order to create Raman instability. Since the min imu m frequency of a light wave in p lasma is pe ω (plasma electron frequency), it is obvious that 0 2 pe ω ω ≥ , based on the matching conditions.
It means / 4 cr n n ≤ , where n denotes plasma density and cr n represents critical p lasma density [5]. Considering the plasma dispersion relation and frequency, we obtain

PIC Simulation Method
PIC simulat ion is emp loyed to trace the trajectory of charged particles, which are calculated on constant-lattice points, in a self-consistent electromagnetic field. The simu lation outline is based on the concept that fields are self-consistent. If the distribution of particles at each time step and the values of the fields at the previous time step are known, the values of the field at the new t ime step can be determined using the equations governing the field (ampere-faraday). Such a field affects the motion of particles and rearranges them. This effect can be determined using the Lorentz force solution. Accordingly, these electromagnetic fields and the distribution of part icles affect each other. It is the self-consistency of electro magnetic fields. The general algorith m of the written 2D electro magnetic PIC code is defined as follo ws: 1) First, the program of the in itial conditions is written by determining the part icles' position and velocity. Then, various parameters such as charge density and current density can be obtained in each time interval and on each lattice point through interpolation. By obtaining charge density and current density, the electromagnetic fields on the lattice points can be calculated by the field equations.
2) The position of the phase points in the phase space is determined through the Newton-Lorentz equations. However, to do such calculations, we need to have the values of the fields on the phase points, which are obtained by extrapolating the values of the fields on the lattice points [6]. The values obtained for the position and velocity of the phase points are used as the initial conditions of the next time step, and all the program steps will fo llo w the same procedure. All the above stages are written in a C++ code.
There are so me criteria to ensure that the code is properly implemented in terms of the number of particles, temporal interval, and spatial interval. The number of part icles should be much h igher than that of lattice points, wh ich ensures the fact that all the latt ice squares constantly averagely contain several particles during the simu lation process. The spatial interval should not exceed the Debye length ( D λ ) in order for the distribution of potential to be computable by charge separation. The temporal interval should be lesser than the inverse plasma electron frequency ( 1 pe ω − ) in order to produce plasma oscillations [7].
The first step to calculate a system of equations is to make them dimensionless to avoid the errors due to very large or small coefficients. The nondimensionalization of the 2D code with respect to the input wave's properties is defined as: where and respectively denote the frequency and wavenumber of the input electro magnetic plane wave, and are the mass and electric charge of particles respectively, is the electric field, is the magnetic field, is the speed of light in vacuu m, is the current density, is the velocity of particles, is the position of particles, and represents time. The interaction between an electro magnetic plane wave and a plasma lobe can be simulated through the Maxwell equations, which are defined in the Gaussian system as: Solving these equations in three dimensions leads to a quite comp lex ensemble of arrays and calculat ions. Hence, we try to simplify the calculat ions through reducing the problem's dimensions. Here, we study the evolution of an electro magnetic plane wave with linear polarization propagated along the x-axis in an interaction with plasma. Accordingly, the vector potential components become (0, , 0) , the electro magnetic-field co mponents are defined as While a plane wave with such properties is passing, if the electrons do not have any initial velocity along the z-axis, the motion of the particles will remain in the xy -plane and there will not be any evolution along the z-axis (Appendix A).

Initial Conditions and the Simulation Model
The init ial conditions and dimensionless data we used in the 2D electro magnetic code are as follows: The length of the simulation bo x along the x-axis and y-axis is chosen 50, the spatial distance between the lattice points along the x-axis and y-axis is 0.1, there is one particle in each cell, the total number of particles is 5

10
× , the number of lattice points along the x-axis and y-axis is 501, the amplitude of the incident electric wave is 0.05, the wavelength of the incident wave is 6, the incident-wave wavenumber is 100/3, the frequency of the dimensional incident wave is the Maxwell equations, we investigate longitudinal-wave growth or RFS instability. Note that the theoretical investigations in this field have been previously performed and this paper aims to compare those results with the results obtained from the PIC simu lation.

Results and Discussion
The simu lation bo x should be initially empty of plas ma to compare the transmission of the plane wave through plasma and vacuum. It can be shown that the propagation of an electro magnetic plane wave through the vacuum is in compliance with the theoretical results and there is not any longitudinal electric wave. These numerical results suggest this code's proper performance [8].
Now, regarding RFS instability, considering the init ial conditions and the simulat ion model, the electro magnetic plane wave first enters the vacuum boundary, and then enters the low-density plasma slab. This wave returns the vacuum region after passing this slab.
Here we consider a thin p lasma slab and investigate the propagation of electro magnetic waves through it. Given an input plane wave (   Figure 1 shows how the longitudinal electric field grows due to the propagation of the input wave through plasma (the number of iteration is 5000). As evident, opposed to the vacuum state, the longitudinal electric field starts to grow and reaches the maximu m of 4

10 − ×
, wh ich is considerable in co mparison to the input-wave amplitude. The maximu m longitudinal-wave amplitude is grown approximately one hundredth of the input-wave amp litude.
Given that /4 cr n n ≤ in this code, the plasma is expected to be transparent and the input wave easily passes through it. The correctness of this case has been already proved.
Note that the growth of the longitudinal-wave amp litude is limited by turbulence resulting fro m holes in plasma. This issue imp lies the plasma nonlinear effects [9][10][11]. In this state, the particles captured by the wave continuously exchange energy with the wave. As can be seen in Figure 1, the longitudinal wave deteriorates as it progresses.
However, since the density of the plasma particles is still much lesser than the critical density, the plasma is transparent and the propagation of the input wave through plasma is the same as what Figure 2 shows. The growth rate of instability is the most important parameter that should be measured. Figure 3 indicates the logarith mic diagram of the longitudinal electric field versus time. As can be observed, although it fluctuates, the field amp litude is substantial and shows a gentle trend for 25 t > .
The slope of the diagram at this reg ion can be considered as the field gro wth rate. Linear least-square fitting is the best technique to calculate the slope of this diagram. The slope of this diagram or the instability growth rate is obtained   Regarding the phase space, Figure 4 shows the maximu m velocity along the x-axis and y-axis in each temporal interval for a 50-length simulation bo x and 5000 iterat ion steps. As evident, the velocity rises along both x-axis and y-axis, and declines for the next temporal intervals.

Conclusions
In this paper, RFS instability was simulated by a 2D electro magnetic PIC code written by C++. First, because we were sure about this code's proper performance, the propagation of electromagnetic waves through vacuum and low-density plasma was investigated. We found that if the med iu m is a low-density plasma, there will not be a difference in the propagation of electro magnetic waves through plasma and vacuum. It means that a low-density plasma acts as a transparent mediu m against electromagnetic waves. Also, it was observed that large-amp litude longitudinal waves are produced and grown when a plane wave passes through a low-density plasma. This case indicates instability in the propagation of electromagnetic waves through low-density plasmas. The slight difference between the calculated growth rate and the theoretical growth rate suggests this code's proper performance.

A. Si mulation of fiel d equations
The Maxwell equations govern fields and electromagnetic waves, whereby we can study the evolution of electro magnetic fields in time and position. Given the spatial-temporal dependence of electro magnetic fields (Faraday's law and Ampè re-Maxwell law), the discretizat ion of these dimensionless equations ) is defined as the following procedure.
Considering the above equations (the temporal derivation of E  and B  at the left-handed side and their existence in the right-handed side), the temporal derivation has been performed using the leap-frog scheme in integer temporal intervals and temporal and spatial semi-integer intervals [11].
It is the time centering method and the accuracy of the temporal derivatives is of the second order. In the 2D space, fields are divided into transverse electric (TE) and transverse magnetic (TM ) fields. All the spatial variab les, including denotes the x-co mponent of the electric field at n t ∆ ,  Assuming that the th charged particle is placed on and , using the particles' current density and charge density, a spatial lattice in which and is considered to obtain the fields. This lattice's size, which is obtained from , should be small enough to avoid numerical errors. To calculate each particle's current and charge density, we consider a surface weight function. These weights are given by: where is one cell's charge density. The Boris method to calculate the current density is defined as: is obtained, consequently. The four nearest points surrounding each particle are used to calculate that particle's current density [12].

B. Simul ation of single-particle 2D motion
The following Newton-Lorentz equation is employed to study the motion of particles.
which can be solved by three simultaneous scalar equations. However, the Boris method is simp ler and o f the second order in wh ich the Newton-Lorentz equation is discretized as: where m denotes the particle rest mass (electron or ion where v′  2) The magnetic field affects v −  , and v +  is obtained through rotation.
3) The electric field affects v +  , and (1/2) n v +  is obtained. Accordingly, velocity at the mo ment (1/ 2) n + can be obtained by knowing the fields at n and velocity at (1/ 2) n − .