A unified gas kinetic scheme for transport and collision effects in plasma

In this study, the Vlasov-Poisson equation with or without collision term for plasma is solved by the unified gas kinetic scheme (UGKS). The Vlasov equation is a differential equation describing time evolution of the distribution function of plasma consisting of charged particles with long-range interaction. The distribution function is discretized in discrete particle velocity space. After the Vlasov equation is integrated in finite volumes of physical space, the numerical flux across a cell interface and source term for particle acceleration are computed to update the distribution function at next time step. The flux is decided by Riemann problem and variation of distribution function in discrete particle velocity space is evaluated with central difference method. A electron-ion collision model is introduced in the Vlasov equation. This finite volume method for the UGKS couples the free transport and long-range interaction between particles. The electric field induced by charged particles is controlled by the Poisson's equation. In this paper, the Poisson's equation is solved using the Green's function for two dimensional plasma system subjected to the symmetry or periodic boundary conditions. Two numerical tests of the linear Landau damping and the Gaussian beam are carried out to validate the proposed method. The linear electron plasma wave damping is simulated based on electron-ion collision operator. Compared with previous methods, it is shown that the current method is able to obtain accurate results of the Vlasov-Poisson equation with a time step much larger than the particle collision time. Highly non-equilibrium and rarefied plasma flows, such as electron flows driven by electromagnetic field, can be simulated easily.


Abstract
In this study, the Vlasov-Poisson equation with or without collision term for plasma is solved by the unified gas kinetic scheme (UGKS). The Vlasov equation is a differential equation describing time evolution of the distribution function of plasma consisting of charged particles with long-range interaction. The distribution function is discretized in discrete particle velocity space. After the Vlasov equation is integrated in finite volumes of physical space, the numerical flux across a cell interface and source term for particle acceleration are computed to update the distribution function at next time step. The flux is decided by Riemann problem and variation of distribution function in discrete particle velocity space is evaluated with central difference method. A electron-ion collision model is introduced in the Vlasov equation. This finite volume method for the UGKS couples the free transport and long-range interaction between particles.
The electric field induced by charged particles is controlled by the Poisson's equation. In this paper, the Poisson's equation is solved using the Green's function for two dimensional plasma system subjected to the symmetry or periodic boundary conditions. Two numerical tests of the linear Landau damping and the Gaussian beam are carried out to validate the proposed method. The linear electron plasma wave damping is simulated based on electron-ion collision op-

Introduction
The Vlasov equation describes time evolution of the distribution function of plasma which consists of charged particles with long-range interaction [? ].
Vlasov showed the difficulties in description of plasma with long-range Coulomb interaction when kinetic theory based on standard transport-collision is applied: 1) Theory of pair collisions disagrees with the discovery by Rayleigh, I. Langmuir and L. Tonks that vibrations exist in electron plasma; 2) Theory of pair collisions without Coulomb interaction will lead to divergence of kinetic term; 3) Theory of pair collisions cannot explain results of experiments by H. Merrill and H.
Webb that electrons scatter anomalously in gaseous plasma [? ]. In gas, binary interaction is taken as the rule. However, in plasma, waves, or organized motion of plasma, are very important because the particles can interact at long ranges through the electric and magnetic forces. Vlasov suggested a source term which contains the Coulomb interaction added in collisionless Boltzmann equation.
The Vlasov equation is a partial differential equation (PDE) for distribution function of particles, which represents probability that particle stays at a specific position and velocity. The acceleration term in equation describes redistribution of particle on particle phase space due to electromagnet force.
The microscopic method is a way to simulate particle motion directly. The trajectory of individual particle is traced and the distribution of particle is quantified with the number of particles per unit volume and per velocity interval.
The macroscopic variables can be obtained by moment results of microscopic variables. P. Degond et al(2010) proposed a particle-in-cell (PIC) method based on Vlasov equation [? ]. Individual particles in a Lagrangian frame are traced in physical space and particle velocity space. Under the Eulerian frame of reference, the mass and momentum conservation equations are obtained with velocity moments of the Vlasov equation [? ]. In this method, the particle details, such as the location, velocity, are exactly captured by the equation of motion. The PIC method also has a advantage on multi-component problems in plasma simulations. The biggest issue of PIC is the huge computational cost when the number of particles increases. For this reason, its application is limited in problems which consist of a small number of particles. The use of penalty function preserves the positivity of the electron distribution functions. The convergence is also greatly improved [? ]. To determine the coefficients in polynomial about phase space, the particle distribution in physical space and velocity space is needed and equation of particle motion is solved at each time step. A rather high computational cost still becomes a barrier for plasma simulation.
The spectral method is a class of methods which can be used to solve the The distribution function is projected on B-splines basis. The particle trajectory is used to decide the coefficient in the sum of basis functions. This method has a good robustness and allows a very large time step. However, because of lack of an effective reconstruction, the conservation is slightly lost, which causes a bit loss of accuracy [? ]. The spectral methods have lower computational cost than FEM. But in problems with complex geometries, it becomes less accurate. The coefficient of polynomial is decided by the particle characteristics. The equation of particle motion in electric field is discretized in time and solved in phase space. The distribution function can be obtained by the summation of the series after solving a linear system [? ]. The time step is restricted by the scale of particle motion and the computational cost does not show an advantage over PIC methods. These conservative methods partly improve the efficiency of PIC method and the conservative property. But they still fail to describe real particle motion in phase space with a low computational cost.
In recent years, the kinetic theory of gases is being studied to analyze the particle dynamics in plasma. At the microscopic level, particles are moving and interacting with each other randomly, but a collection of particles can be de-scribed with deterministic macroscopic properties. Kinetic theory analyzes the particle dynamics based on statistics, which links between the microscopic and macroscopic evolutions together. This method avoids huge computational cost since it does not need to trace individual particle. Based on the kinetic model, the complex physical process of nanosecond-pulse electrical discharges including the variations of multiple physical parameters are accurately captured. The defect of this method is the way of reconstruction at the cell interface. The characteristic curves for particle transport are lost in flux evaluation. An artificial approximation has replaced the physical situation at the cell interface. So a wrong dissipation is brought in because of the inaccurate reconstruction for particle state at the cell interface.
In current paper, a unified gas kinetic scheme (UGKS) is applied to sim- the Maxwellian distribution function can be easily used to integrate the Vlasov equation in a finite volume. But in plasma flow, we cannot make an assumption that particles stay in an equilibrium state. The real gas distribution function in the highly non-equilibrium region is never able to be described by a Maxwellian velocity distribution. In this paper, the particle velocity space is discretized so that distribution can be exactly described in all flow regimes [? ]. The flux reconstruction at the cell interface is treated as the Riemann problem along characteristic curves for particle transport [? ]. The free transport and interaction between particles are coupled so that the dissipation in the transport process is controlled by the source of the long-range interaction rather than numerical time step. With this technique, the artificial dissipation as the one in Ref.
[? ] will not be brought in numerical process.
For computational cost, the UGKS is a much more efficient method than PIC method. The PIC method includes the four typical procedures: 1) Equation of particle motion is integrated; 2) Properties of particles are interpolated to field mesh; 3) Field is reconstructed on mesh points; 4) Flow field is interpolated from mesh to particle locations. The information is communicated between particles and macroscopic field at every time step. In real system being studied, the number of particles is often extremely large, resulting in poor communication efficiency. The simulation becomes inefficient and even impossible at all. Different from the way of description for flows in PIC method, the UGKS describes flow motion using particle distribution function. It gives the number density of particles in the six-dimensional phase space. This method is a statistical method based on particle motion. Instead of following individual particle, the distribution function is used to show probability of particles which locate in a spatial point and a certain velocity interval. The macroscopic variables can be obtained by the integral of the microscopic variables in particle velocity space. Internal motions of particle, such as rotation and vibration, are taken into account through extra internal degrees of freedom [? ]. The UGKS avoids huge computational cost produced by solving the equation of particle motion and tracking particle as stated in PIC method. This high fidelity and efficiency in numerical simulation make the UGKS a better way to solve plasma problems. In solving the Vlasov equation, the UGKS can overcome the negative effect of truncation error caused by the approximation of the distribution function in FEM. In FEM the phase space is divided into a collection of sub-domains.
A set of equations are used to represent discretized PDE in every subdomain.
All sets of element equations are finally combined for solution of the whole discrete phase space. In construction of subdomains, the distribution function is approximated by polynomial. FEM tries to minimize the truncation error from this approximation. So, integration of inner product of residual and weight functions is set to zero. In this procedure, error is caused by trial functions for distribution function in every subcell of phase space. In the UGKS, the Vlasov equation is integrated and the flux is evaluated at cell interface of phase space. This flux is obtained by the solution for the local Riemann problem [? ]. At the cell interface, the direction of propagation of particle information is simulated.
And then the flux is reconstructed by using difference biased in the direction determined by the sign of the characteristic speeds [? ]. The particle velocities propagated to cell interface are the characteristics in the UGKS, which gives a correct representation for physical motion. The only source of error comes from the dissipation which depends on the resolution of a mesh cell.
The UGKS also has a more concrete physical foundation than spectral method. Spectral method writes distribution function in the Vlasov equation as a sum of basis functions. The solution must satisfy the control equation as well as possible so coefficients before basis should be carefully chosen. In the case that the solution is smooth, the spectral method has a high order accuracy [? ]. However, a highly non-equilibrium state will challenge this method because a large error appears in a strong discontinuity at cell interface of phase space.
The approximation using polynomial can not achieve a fidelity when the distribution function deviates from the Maxwellian distribution in non-equilibrium state. In physics, the distribution function can have arbitrary geometric features in plasma so polynomial often fails to capture actual distribution in phase space. In equilibrium state, the Vlasov equation can be integrated under continuous particle velocity space in the flux evaluation based on FVM. Because in this regime, analytical formulas for gas distribution function is well-defined.
However, a plasma has properties unlike those of fluid flow. Although particles are unbound, their transport is not entirely free. When a charged particle is moving, time-varying electric field is produced around it. In plasma, the movement of a charged particle affects and is affected by the general field created by the movement of other charged particles. This process causes collective behavior with many degrees of variation [? ? ]. Collision interactions are often weak in hot plasmas and external forcing can drive the plasma far from local equilibrium and lead to a significant population of unusually fast particles. For these reason, we can not have a confidence that we can still approximate distribution function of particles in plasma under continuous velocity space. In the UGKS, the particle velocity space is discretized into a number of points which represent every dynamical state. The Vlasov equation is split into PDEs for distribution function at each velocity point. In fact, we solve the Vlasov equation for distribution function at each point of discrete phase space. So, a high fidelity is obtained when distribution function is solved at each time step, which is the advantage over spectral method. Theoretically, the physical content in the UGKS is far too richer than methods based on numerical approximation using polynomials.
Unlike neutral gas, plasma consists of a significant number of charge carriers, which makes it electrically conductive. Because of this feature, plasma responds strongly to electromagnet field. Plasma does not have certain shape or certain volume if not enclosed in a container, but its motion and state are

Numerical methods
The plasma Vlasov equation can be written as where f is the particle distribution function of the space x, time t, particle velocity u, and acceleration a due to electromagnetic field. The relation between the macroscopic variables W with the microscopic variables ψ is where ρE = 1 2 ρ U 2 + V 2 + K+2 2λ , and ψ = (1, u, v, 1/2(u 2 + v 2 + ξ 2 )) T . The pressure and density can be related with λ as p = ρ/2λ. u = (u, v) is the particle velocity and ξ is the internal freedom of gas.
In discrete particle velocity space, Eq. (2) can be written as The particle velocity u is discretized with intervals ∆u and ∆v. In finite volume, we integrate the Vlasov equation in physical space for each discrete velocity [? ]. The result of integration can be written as where i is the cell index of a unstructured mesh, j and k present the index of particle velocity. The distribution function at the n + 1 time step can be obtained through Eq. (4).
In current work, the charged particles interact with each other by longrange electromagnet force. The advection term can be evaluated as numerical flux along characteristic curves. In the Vlasov model, The distribution function around the cell interface can be reconstructed as [? ] The electromagnet force term is obtained by central difference method. For a dimension, difference between f and particle velocity u and v can be written At the boundary of particle velocity space for two dimensional plasma, we where N u and N v are the total numbers of discrete particle velocities in each dimension.
Next, we will introduce the process of solving the electromagnet Poisson's equation. Solving the Poisson's equation amounts to finding the electric potential Φ for a given charge distribution. The mathematical details behind the Poisson's equation in electrostatics are described by Gauss's law for electricity [? ]. The Poisson's equation applied in our work is written as where ρ c is the charge density and ε is the permittivity of the medium. Eq. (8) is solved using the Green's function [? ]. First, we find the particular solution.
where N c represents the number of cells, i represents the cell being studied, j represents cells in discrete computational domain and j = i. Q j represents charge in cell j. G (x i , x j ) is the Green's function for the Poisson's equation, which can be written as In the Green's function, x 1 and x 2 represent radius vectors for two locations of charges. To give a unique Green's function, symmetry or periodic boundary conditions should be implemented [? ]. According to the Faraday's law, the electric field from charged particles is a conservative vector field. The electric potential Φ can be defined [? ], such that Eq. (11) is used for electrostatic part of electric field. The electrodynamic part of electric field is induced by magnet field due to motion of charged particles.
If magnet field is taken into consideration, the electric field induced by charged particles can be written as where A is the magnetic vector potential [? ]. The relation between magnet field B and A can be written as The component of particle acceleration due to the electrostatic part of electric field obtained by particular solution via the Green's function and electrodynamic part of electric field induced by magnet field is computed as, where q i is the charge of the particle in cell i, m is the mass of the particle.
Besides the electric field induced by the charged particles, the external electric field E appl also contributes to variation of distribution function. So the acceleration term in Eq. (1) includes two parts.
where q m is charge to mass ratio. In this paper, the non-dimensional parameter satisfies relation q m = 1. After substituting Eq. (5), Eq. (6) and Eq. (15) into Eq. (4), we can update the distribution function at n + 1 step. The above procedure can be repeated in the next time level.

Linear Landau damping
The Landau damping is the effect of damping of longitudinal space charge waves in plasma [? ]. It prevents an instability from developing and creates a region of stability in the parameter space [? ]. Energy exchange between electromagnetic wave in phase space and charged particles in plasma is the cause of Landau damping. In evolution, particles interact strongly with the wave [? ].
In current work, we simulate the Landau damping based on the UGKS to study the way in which electrons interact and exchange energy with electromagnetic field. This case is used to test the accuracy of the UGKS for solving the Vlasov equation in describing time-variant electromagnetism parameters such as electric field and electric potential energy. First, an electromagnetic wave should be added to flow field. This wave is caused by an initial disturbance of distribution function [? ], which can be written as, where g is a Maxwellian distribution, In current work, 2D initial condition is set to where λ = 0.5, U = V = 0 , α = 0.05 and the wave numbers k x = k y = 0.5 in current problem.
In the linear Landau damping case, a periodic boundary condition is applied.
The four dimensional phase space contains 64 points per dimension. The particle velocity is truncated at 6.0. The boundary condition for 2D linear Landau damping is presented in Fig. 1. The length for computational domain in each dimension is set to be 4π. The initial electronic density distribution is given in Fig. 2. The constant time step ∆t is chosen as 0.01. The evolution of density contour is presented in Fig. 3. The energy of electric field is being carried away by electrons moving at the phase space. To show this process, the evolution of symmetric electric field and electric energy will be given. First, we computed the average of symmetric electric field. The logarithm scale of it is used for damping process of electric field, which is presented in Fig. 4. Now, the electric potential energy is computed to show the process of energy transport and conversion during interaction between the electric field and electrons. The electric potential energy is a potential energy that results from coulomb forces between charges. This kind of energy is associated with the configuration of a defined system which contains a certain number of charged particles. In physics, the electric energy of a system is the energy required for assembling charges from an infinite distance by bringing them close together. If an object has electric potential energy U E , two key elements are crucial: 1) this object keeps its own charge; 2) it stays at a relative position to other electrically charged objects.
The evolution of potential energy in systems with time-variant electric fields is given in Fig. 5.
In current work, we obtain the same evolution process of the electric field and electric energy as Ref. [? ]. Many works have been done to study the physical picture of Landau damping. Linearized theory is the simplest and rather complete way [? ], but it is not physically reasonable because interaction between charged particles and electromagnet field is coupled with particle transport in the Vlasov equation. Linearized method is not appropriate to be used in this nonlinear system. However, nonlinearity has been a longstanding problem. To model nonlinear level in the Vlasov equation, a class of exponentially damped solutions of the Vlasov-Poisson equation is proposed [? ]. This method fails to explain the mechanism of energy transport although it approximates the damping curves in a mathematical way. The UGKS shows a great advantage on correctly representation for the Landau damping with a much lower computational cost than particle based methods.

Nonlinear Landau damping
In linear Landau damping case discussed above, the non-equilibrium phenomenon is still not very significant because the initial perturbation of density is very small. In this example, we increase amplitude of the initial perturbation of density in 1D space. The perturbation rate is set to be α = 0.5. The wave number and periodic length remain the same k = 0.5 and L = 4π respectively.
The initial particle distribution function reads: In UGKS, the discrete particle velocity space is applied. The 1D velocity space is truncated at [−6, 6] with 128 intervals. When the velocity space is being discretized, we used cosine to decide the positions of discrete velocity points so that a high resolution can be obtained near the peak of distribution function.

This equation is written as
where v max is the maximum of particle velocity.
The evolution of the kinetic entropy is taken as a benchmark solution. Traditional FEMs always overestimate this value.The particle based methods may be an alternative, but computation is overwhelming. The UGKS is able to obtain the same accuracy as particle based method but calls for much lower The results of UGKS fit well with the ones of the particle based method.
With a much lower computational cost, the UGKS has the same advantage in accuracy as particle based method. The UGKS has a good conservation property and good accuracy in reconstruction so that spurious dissipation is avoided. The nonlinear Landau damping case is a highly non-equilibrium phenomenon. To show this, particle distribution functions at different time steps at the center of space are presented in Fig. 8. Fig. 8 shows the distribution function in velocity space. To give a better description of particle distribution in phase space, the distribution function in the (x, u) space is presented in Fig. 9.
In nonlinear Landau damping case, the nonlinear effects are too important so previous traditional Galerkin methods always fail to simulate this highly non-equilibrium case accurately. In UGKS, the particle behavior is accurately obtained by solving the Vlasov equation with discrete particle velocity space.

Gaussian beam
In optics, the Gaussian beam is a monochromatic electromagnetic radiation which has a Gaussian intensity profile. As a result, it has a transverse magnet electromagnet field given by the Gaussian function. Such a beam can be expanded and focused through lens thus becomes different Gaussian beam at different time step. According to the theory of quantum mechanics, electrons can both interact with other particles and be diffracted. It has properties of particle and wave [? ]. For Gaussian beam case, a beam of electrons is a good model since both interaction between particles due to electromagnet force and intensity of light should be simulated. In current paper, we choose a electron beam as model to be deformed by the electromagnet field. And the external electric field plays a role of lens. In evolution of beam, the external electric field has an opposite direction to electric field induced by electron itself. In the Vlasov equation, the acceleration term is from the electric field which includes two parts: Electric field E self induced by electron given by the Poisson's equation and external electric field E appl . The external electric field is decided by self-consistent field method (SCF) using initial electric field induced by electrons. In current work, the external electric field is linear with respect to spatial coordinates. The relation between E self and E appl can be written as where ω represents the difference between the initial self-consistent electric field due to particle charges and the linear external electric field. Its value is decided by a self-consistent domain.
All values in Eq. (23) are obtained by the root mean square (RMS) thermal velocity v th . The RMS results for space and particle velocity are computed by For Gaussian beam case, the initial self-consistent electric field is linear in space. So we can obtain the two unknown variables ω and a using RMS thermal where r a is the radius of beam. Now we introduce a tune depression factor χ = ω0 ω such that we can adjust external electric field with initial self-consistent electric field. The linear external electric field whose direction is opposite to initial self-consistent electric field can be written as In evolution of a beam, the external electric field E appl is used to focalize the beam. The dimensionless form of Eq. (1) can be rewritten as The negative sign " − " is because of the negative charge carried by electron.
The 2D initial condition is written as In kinetic theory, λ can be obtained by [? ] where m is the particle mass, k B is the Boltzmann constant, T is the temperature. The mesh for beam case has a round boundary, which is shown in Fig. 10.
4032 mesh cells are used to discretize the computational domain.
Absorbing boundary condition (ABC) is applied in Gaussian beam case. The buffer zone is added between the computational domain and boundary with a range 0.5 2 < x 2 + y 2 < 1.0 2 . The flux damps along warp of a circular zone in buffer zone.
where δ is damping rate, r o and r i are the two location on a warp line respectively. The damping rate can be computed as r i and r o represent the two limit of buffer zone. In our work, we have r i = 0.5, The location for F ro and F ri is described in Fig. 11.
In evolution, the electric field induced by charged particles and particle free transport will expand the beam. Applied external electric field focalizes the beam. So there is a process in which expansion and contraction take place alternately. We describe evolution process of beam using electron density contour at different time step. The particle velocity is truncated at 6.0 with 64 points.
We simulated the evolution of a Gaussian beam and compared our results with the work of Ref. [? ], which proves the accuracy that the UGKS obtains in plasma simulation. Now, let's compare the computational cost with particle based methods introduced in Ref. [? ]. The test model is put in a 64 × 64 × 64 × 64 phase space.
We compared computation time of the UGKS and particle based scheme under the same phase space in Table 1. The comparison results show a high efficiency achieved by the UGKS in solving the Vlasov equation. Although electron has wave-particle dualism, in this paper, we only study Gaussian beam under an electrodynamics frame. Based on the Vlasov-Poisson system, the motion of electrons is described in phase space with particle distribution function.
Electrons are moving in space because of existence of electromagnet force. So the peak value of intensity of the beam varies with time in Fig. 12. To obtain a correct evolution process of beam, the term of transport and acceleration in the Vlasov equation should be accurately coupled. Although the wave nature for electrons is out of the range in this paper, the Schrödinger equation can be used to construct a numerical method for quantum results of Gaussian beam case in the future.

Linear Electron Plasma Wave Damping
In this section, the effect of collisions on linear damping of spatially one- where the collision operator −C ei f is given by where ν ei,th denotes the thermal electron-ion collision frequency and u th is velocity of particle thermal motion. Tensor P is defined by In the EPW case, the wavenumber is set to be 0.3. In the low perturbation amplitude regime, the initial electron distribution is where α = 0.0001 and k = π 250 . The 1D computational domain is discretized with 64 grid points. In microscopic velocity space, 128 discrete points are used.
To compared with the result of Ref. [? ], the thermal electron-ion collision frequency ν ei,th is set to be 0.05. The amplitude of the perturbation's electric field E is plotted on a log scale as a function of time in Fig. 13.
Where ω pe is the circular frequency of initial perturbation. The result obtained from present methods is consistent with Ref. [? ]. The decaying process fit well with exponential decaying line. Because of the initial perturbation of particle distribution, the electric field is also damped in oscillatory process. In traditional numerical method based on continuous assumption, its very hard to capture these phenomena for the lack of direct modeling of particle motion. For particle based method, the huge number of particles always produces unacceptable computational cost. In UGKS, particle motion is described in modeling of flux through cell interface under discrete microscopic velocity space. For every discrete distribution function, computation cost only comes from process of FVM. But accuracy is the same as particle based method since all particles are considered based on statistical concept using distribution function. With suitable collision model, different state of plasma can be accurately obtained.

Conclusions
In this paper, we solve the Vlasov equation based on the UGKS. This approach gives a full representation for particle moving in phase space with longrange electromagnet interaction. An integral solution of the kinetic model under discrete physical and particle velocity space is applied. In non-equilibrium state, this method can describe distribution function with discrete phase space. Compared to Lagrangian method, the UGKS has a much lower computational cost but achieves the same accuracy in plasma simulation. Many methods for solution of the Vlasov equation, such as FEM and spectral method, still have to solve equation of particle motion to decide the coefficients in basis functions for particle velocity dimensions. As a result, the time step is restricted according to dissipative length scale determined by the physical transport. Time step in current work can be decided by CFL strategy without spurious dissipation. In order to simplify the solution of the Poisson's equation by using the Green's function, we restrict the boundary conditions to the periodic or symmetric boundary conditions. The method presented in this paper can be used in simulation for flows of charged particles in plasma. In multi-scale problems of plasma flows from free transport regime to collisional effect, UGKS can give pretty good description of physical evolution. Compared to particle based method, computational efficiency has also been greatly improved. From the results via the UGKS, it is proved to be a better way than previous particle methods to study plasma problems because of its good accuracy and low computational cost.