Simulation of Velocity Evolution of a Cold Collision-less Non- Magnetised Plasma by Particle-in-Cell Method

This work presents simple numerical simulation algorithm to analyse the velocity evolution of high density non-magnetized glow discharge (cold) collision-less plasma using Particle-in-Cell (PIC) method. In the place of millions of physical electrons and background ions, fewer particles called super particles are used for simulation to capture the plasma properties such as particle velocity, particle energy and electrical field of the plasma system. The plasma system which is of interest in this work is weakly coupled plasma having quasi-neutrality nature. Simulation results showed symmetric velocity distribution about zero with slight left skewness, indicating static system. The order of directional velocity of individual particle seems to agree with the input electron temperature of the considered plasma system. The particle and field energy evolution were observed having fluctuations about zero which indicates that the system is equilibrating. This work marks the preliminary work to study the transport of plasma species in plasma column of gliding arc discharge.


Introduction
Computer aided simulation is always a useful tool to predict results of expensive or tedious experiments, provided the formulation of simulation is valid by scientific theory. Due to complex phenomena involved in plasma, computer aided simulations have been very useful tool either to deepen the understanding of the underlying plasma physics or to explore new chemical reactions. Plasma simulation can be classified into two areas namely (i) simulations trying to understand plasma reaction kinetic and (ii) simulations based on plasma fluid description [1]. Fluid simulation solves the magnetohydrodynamic (MHD) equations assuming necessary transport properties. Kinetic description simulates plasma species interactions through electromagnetic (EM) field by solving plasma kinetic equations such as Vlasov or by Particle-In-Cell method (PIC) [1]. PIC simulation is used to predict laser-plasma interactions, electron acceleration and ion heating, ion-temperature-gradient and instabilities in fusion and dusty plasmas. In PIC, the motion of collective charged particles under externally applied field is computed. The system is represented using a 2D grid phase space by a small number of finitesize particles. Here, the direct coulomb law approach is replaced by obtaining field from charge density at the grid points [1]. Hence, the particles interact through the potential, only at distances beyond the overlap distance and corrects the effect of fewer particles at closer distances by the reduced interaction potential. With this approach, even few particles (electrons and ions) in a collision-less plasma simulation can represent millions of physical electrons and ions of an actual system. Hence, these particles in PIC simulation are called as super particles [2]. In this work, a weakly coupled plasma system is assumed and using PIC method, it is attempted to estimate the velocity distribution and equilibration with super particles rather than simulating the exact plasma density to verify the super particle approach.

Model Formulation and Simulation Methodology
A weakly coupled plasma system as defined in literature [2] i.e. cold glow discharge system is considered in this work. The system was represented by ions and electrons which are considered as particles in simulation. Figure 1 shows the grid terminology used in this work for numerical simulation where the particles were loaded and simulated. Figure 1. Grid terminology used for particle-in-cell simulation Before setting up the algorithm, numerical stability and assumption of the characteristic plasma system was verified using the conditions given in equation 1 to 7 as recommended in literature [1][2][3]. Table 1 shows the values of the parameters that represent the considered plasma system and the ones which are used to verify the conditions for numerical stability which are expressed as follows: where L is length of entire square grid space, is the Debye length, ℎ is size of single cell (Δx, Δy), is the plasma frequency, is the velocity of light and is the timestep used in the simulation. The following conditions in equation 5 to 7 were also satisfied to confirm that the parameters in this simulation represents weakly coupled plasma with quasi-neutrality.
These conditions are: (1) plasma density and electron temperature of system is relevant to the type of plasma to be simulated; (2) particle number density considered for the simulation (ND) >>1, which can be expressed as: This particle number density was reduced in the simulation (super particles) for 2D using the following expression: where rate of change of is zero; (3) ratio of mean interparticle potential energy to the mean plasma kinetic (thermal) energy was kept << 1. This ratio was defined as a plasma parameter (Λ) which was expressed as, The particles in simulation were equally divided among electrons and ions i.e. (ne) = (ni), with ions assumed as background charge at rest due to its low charge to mass (q/m) ratio compared to that of electrons [4]. These particles were having zero initial velocity (cold) at loading and later the initial velocity was calculated during its transition to first time step using leapfrog algorithm [5]. Figure 2 shows the flow of chart of the simulation and the relevant formulation equations used are shown in Supplementary-material. The coding was done in Matlab R2016b and simulated. The simulation details are provided in Table 1. The simulation was run twice for the same time step to check consistency of the velocity evolution and the distribution curve as the particle initial loading was random. The distribution fit of velocity of the particles was plotted using Matlab's 'Distribution Fitting App'   Figure 3a and 3b shows the plasma particle velocity distribution, in x and y directions respectively, obtained from the simulation for a high density non-magnetized glow discharge (cold) collision-less plasma system. As can be seen from figure 3a and 3b, the velocity distribution was symmetric about zero, but with a slight left skewness and the magnitude of the velocity was of the order of 105 m/s. This symmetric velocity distribution around zero is a typical signature of static cold plasma systems. For static equilibrium of particles, the net force on particles is zero due to cancellation of forces which have equal magnitude but acting in opposite directions. The symmetric velocity distribution of cold collision-less plasma system simulated in this work hence represents that the system is in static equilibrium.

Results & Discussion
As can be seen in figure 3, for the cold collision less plasma system studied, the super particles attained a velocity of order of 105 m/s, which is in agreement with the input electron temperature of plasma system [5]. This result demonstrates that the velocity attained by the particles in a cold collision-less glow discharge plasma system represented the input conditions and can be estimated by simulating very few particles instead of the entire 1016 m-3 of ions or electrons. Hence, super particles can be efficiently used to estimate plasma parameters using PIC simulation.  Figure 4a and 4b shows the evolution of particle energy and electric field energy over the timesteps for the simulated cold collision-less glow discharge plasma.
As can be seen from figure 4, there is no directionality in the energy distribution in both 4a and 4b. That is, both particle energy and electric field energy fluctuates around zero, which is a demonstration that the system is in equilibrium [6]. This result signifies that for the boundary condition simulated in this work, where the electrical potential is kept at zero, the charged particles within this boundary does not have net driving flow and energy. Typically, by altering the boundary condition potential, the motion of individual particles can be traced, and the transport properties of the species considered in the simulation can be estimated. Hence, PIC is used in many plasma applications for predicting the motion or behaviour of particles of interest that can respond to the electric field.

Conclusion
In this work, by using PIC simulation method, a glow discharge plasma system was simulated using few super particles (1800) instead of actual density (4×1010 for 2D) of the considered system. It was verified that the evolution of velocity space was in the same order of a glow discharge system. The energy fluctuations about zero indicated that system is equilibrating. This shows that the plasma parameters can be estimated for the considered system by simulating very few particles in simulation called "super particles". Though the concept is established before, this work marks the first step to test indigenous codes developed from initial stages for particle-in-cell simulations for estimating transport properties of the species in a plasma volume.