## 2 Numerical Methods  
The Boltzmann BGK equation is used for UGKS to simulate plasma. Other collision term can also be adapted to this method via modelling particle relaxation time. The BGK equation can be written as  
\[\frac{\partial f}{\partial t}+u\frac{\partial f}{\partial x}+a\frac{\partial f }{\partial u}=\frac{g-f}{\tau}, \tag{1}\]where \(f\) is the particle distribution function of the space \(\mathbf{x}\), time \(t\), particle velocity \(\mathbf{u}=(u,v)\), collision time \(\tau\), and acceleration \(\mathbf{a}\) due to electromagnetic field. \(g\) is the equilibrium distribution function,  
\[g=\rho\left(\frac{\lambda}{\pi}\right)^{\frac{E+D}{2}}e^{-\lambda\left[(\mathbf{u} -\mathbf{U})^{2}+\xi^{2}\right]}, \tag{2}\]  
where \(K\) denotes internal freedom degree and \(D\) denotes dimension, \(\xi\) is the internal freedom of gas, \(\mathbf{U}\) is the macro velocity and \(\rho\) is the density of plasma. In this paper, only two dimensional cases are considered, so the relation between the macroscopic variables \(\mathbf{W}\) with the microscopic variables \(\mathbf{\psi}\) is  
\[\mathbf{W}=\left( \begin{array}{c}\rho\\ \rho U\\ \rho V\\ \rho E\end{array} \right)=\int\psi fdudvd\xi, \tag{3}\]  
where \(\rho E=\frac{1}{2}\rho\left(U^{2}+V^{2}+\frac{K+2}{2\lambda}\right)\), and \(\mathbf{\psi}=(1,u,v,1/2(u^{2}+v^{2}+\xi^{2}))^{T}\). The pressure and density can be related with \(\lambda\) as \(p=\rho/2\lambda\). In discrete particle velocity space, Equation (3) can be written as  
\[\mathbf{W}=\sum\sum\mathbf{\psi}f\Delta u\Delta v\Delta\xi. \tag{4}\]  
where \(\Delta u\), \(\Delta v\) and \(\Delta\xi\) are microscopic velocity intervals. In this discrete integration, the limit is decided by three times of initial Maxwell variance \([-3\sigma,3\sigma]\). In the unified gas kinetic scheme, at the cell interface the solution \(f_{cf}\) for Boltzmann BGK equation is constructed from an integral solution using the method of characteristics  
\[f_{j,k}\left(\mathbf{x}_{cf},\mathbf{u},t\right)=\frac{1}{\tau}\int_{0}^{t}g\left(\bm {x}^{\prime},\mathbf{u}^{\prime},t\right)e^{-(t-t^{\prime})/\tau}dt^{\prime}+e^{- t/\tau}\tilde{f}_{j,k}\left(\mathbf{x}_{cf}-\mathbf{u}t,\mathbf{u}-\mathbf{a}t,0\right). \tag{5}\]  
Here \(\mathbf{x}^{\prime}=\mathbf{x}_{cf}-\mathbf{u}t+\frac{1}{2}\mathbf{a}t^{2}\) and \(\mathbf{u}^{\prime}=\mathbf{u}-\mathbf{a}t\) are particle trajectory and velocity, respectively. \(\tilde{f}\) is the initial distribution function. Subscript \(cf\) denotes cell interface, \(j\) and \(k\) represent the index of particle velocities \(u\), \(v\). Inside each control volume of physical space and discrete particle velocity space, the gas distribution function at the beginning of time step is computed using a linear reconstruction [17]  
\[\tilde{f}_{j,k}\left(\mathbf{x}_{cf}-\mathbf{u}t,\mathbf{u}-\mathbf{a}t,0\right)=\left\{ \begin{array}{ll}f_{cf,jk}^{l}-\left(\frac{\partial f}{\partial\mathbf{x}} \right)_{cf,jk}^{l}\cdot\mathbf{u}t-\left(\frac{\partial f}{\partial\mathbf{u}} \right)_{cf,jk}^{l}\mathbf{a}t,&\mathbf{u}\cdot\mathbf{S}_{cf}\geq 0,\\ f_{cf,jk}^{\prime}-\left(\frac{\partial f}{\partial\mathbf{x}}\right)_{cf,jk}^{ \prime}\cdot\mathbf{u}t-\left(\frac{\partial f}{\partial\mathbf{u}}\right)_{cf,jk}^{ \prime}\mathbf{a}t,&\mathbf{u}\cdot\mathbf{S}_{cf}<0,\end{array} \right. \tag{6}\]  
where \(l\) and \(r\) represent left and right of cell interface. In a two dimensional configuration, \(i\) represents cell index of an unstructured mesh. Under unstructured mesh, spatial derivatives of distribution functions can be obtained via least square method. The approximation of the partial derivative of \(f\) with respect to the particle velocity \(\mathbf{u}\) can be written as  
\[\left(\frac{\partial f}{\partial u}\right)_{i}=\frac{f_{i,j+1,k}-f _{i,j-1,k}}{2\Delta u}, \tag{7}\] \[\left(\frac{\partial f}{\partial v}\right)_{i}=\frac{f_{i,j,k+1}- f_{i,j,k-1}}{2\Delta v}.\]  
For a second order scheme, it is an appropriate approach.