## 2 Numerical Methods
An efficient fast spectral method can also be implemented in UGKS solver for high order accuracy [42]. At the boundary of particle velocity space for two dimensional plasma, distribution function is set to be zero. For integral term of equilibrium distribution function in Equation (5), continuous particle velocity space is applied to evaluate this Maxwell distribution function. Around cell interface \(cf\), it can be expanded as  
\[g\left(\mathbf{x}^{\prime},\mathbf{u}^{\prime},t\right)=\left\{ \begin{array}{l}g_{0} \left(\mathbf{x},\mathbf{u},0\right)\left[1-\left(t-t^{\prime}\right)\bar{\mathbf{a}}^{l} \cdot\mathbf{u}-\left(t-t^{\prime}\right)\bar{\mathbf{b}}^{l}\cdot\mathbf{a}+\bar{A}t^{ \prime}\right],\ \mathbf{u}\cdot\mathbf{S}_{cf}\geq 0,\\ g_{0}\left(\mathbf{x},\mathbf{u},0\right)\left[1-\left(t-t^{\prime}\right)\bar{\mathbf{a}}^ {r}\cdot\mathbf{u}-\left(t-t^{\prime}\right)\bar{\mathbf{b}}^{r}\cdot\mathbf{a}+\bar{A}t^{ \prime}\right],\ \mathbf{u}\cdot\mathbf{S}_{cf}<0,\end{array} \right. \tag{8}\]  
where \(g_{0}(\mathbf{x},\mathbf{u},0)\) is local Maxwellian located at cell interface which reads  
\[g_{0}\left(\mathbf{x},\mathbf{u},0\right)=\rho_{0}\left(\frac{\lambda_{0}}{\pi}\right) ^{\frac{K+D}{2}}e^{-\lambda_{0}\left[\left(\mathbf{u}-\mathbf{U}_{0}\right)^{2}+\xi^{ 2}\right]}, \tag{9}\]  
where \(\bar{\mathbf{a}}\) and \(\bar{A}\) are derivatives of Maxwellian in space and time. They can be obtained from Taylor expansion of Maxwellian. \(\bar{\mathbf{b}}\) is the derivative term of Maxwellian \(g\) with respect to the particle velocity \(\mathbf{u}\). In a two dimensional configuration, they have the following form,  
\[\bar{a}_{n} =\bar{a}_{n}^{1}+\bar{a}_{n}^{2}u+\bar{a}_{n}^{3}v+\bar{a}_{n}^{4} \frac{1}{2}\left(u^{2}+v^{2}+\xi^{2}\right), \tag{10}\] \[\bar{a}_{t} =\bar{a}_{t}^{1}+\bar{a}_{t}^{2}u+\bar{a}_{t}^{3}v+\bar{a}_{t}^{4} \frac{1}{2}\left(u^{2}+v^{2}+\xi^{2}\right),\] \[\bar{A}_{n} =\bar{A}_{n}^{1}+\bar{A}_{n}^{2}u+\bar{A}_{n}^{3}v+\bar{A}_{n}^{4} \frac{1}{2}\left(u^{2}+v^{2}+\xi^{2}\right),\]  
where \(n\) and \(t\) represent the normal and tangential direction component respectively. Macroscopic variables in Equation (9) are determined by the compatibility condition of BGK model. The conservation constraints are given as  
\[\mathbf{W}_{0}=\int\psi g_{0}\left(\mathbf{x},\mathbf{u},0\right)d\Xi=\int\mathbf{\psi}\left( g^{l}H\left[u_{n}\right]+g^{r}\left(1-H\left[u_{n}\right]\right)\right)d\Xi, \tag{11}\]  
where \(H\) is Heaviside function defined by  
\[H\left[x\right]=\left\{ \begin{array}{ll}0,&x<0,\\ 1,&x\geq 0.\end{array} \right. \tag{12}\]  
The normal derivatives term can be obtained by the relation between the spatial derivatives of the distribution function and the conservation variables,  
\[\frac{\mathbf{W}_{0}-\mathbf{W}^{l}}{\rho_{0}\Delta x^{l}} =\int\bar{a}_{n}^{l}\psi g_{0}\left(\mathbf{x},\mathbf{u},0\right)d\Xi, \tag{13}\] \[\frac{\mathbf{W}^{r}-\mathbf{W}_{0}}{\rho_{0}\Delta x^{r}} =\int\bar{a}_{n}^{r}\psi g_{0}\left(\mathbf{x},\mathbf{u},0\right)d\Xi.\]  
For the two dimensional problems, the the normal and tangential component of \(\bar{\mathbf{b}}\) in Equation (8) read  
\[g_{0}\left(\mathbf{x},\mathbf{u},0\right)\bar{b}_{n} =\frac{\partial g\left(\mathbf{x},\mathbf{u},0\right)}{\partial u}=-2 \lambda\left(u-U\right)g_{0}\left(\mathbf{x},\mathbf{u},0\right), \tag{14}\] \[g_{0}\left(\mathbf{x},\mathbf{u},0\right)\bar{b}_{t} =\frac{\partial g\left(\mathbf{x},\mathbf{u},0\right)}{\partial v}=-2 \lambda\left(v-V\right)g_{0}\left(\mathbf{x},\mathbf{u},0\right).\]  
Now, we have determined all parameters of the initial distribution function and the Maxwell distribution function in Equation (5). After substituting these parameters into the formula of gas distribution function, we have \[