## 2 Numerical Methods
\begin{split} f_{j,k}\left(\mathbf{x}_{cf},\mathbf{u},t\right)=& \left(1-e^{-\frac{t}{t}}\right)g\left(\mathbf{x}^{\prime},\mathbf{u}^{ \prime},t\right)+\left(e^{-\frac{t}{\tau}}\left(t+\tau\right)-\tau\right) \left(\overline{a}_{n}^{l}\mathbf{u}\cdot\mathbf{n}+\overline{b}_{n}^{l}a_{n}\right)H \left[\mathbf{u}\cdot\mathbf{n}\right]g\left(\mathbf{x}^{\prime},\mathbf{u}^{\prime},t\right) \\ &+\left(e^{-\frac{t}{\tau}}\left(t+\tau\right)-\tau\right)\left( \overline{a}_{n}^{r}\mathbf{u}\cdot\mathbf{n}+\overline{b}_{n}^{r}a_{n}\right)\left(1 -H\left[\mathbf{u}\cdot\mathbf{n}\right]\right)g\left(\mathbf{x}^{\prime},\mathbf{u}^{\prime },t\right)\\ &+\left(e^{-\frac{t}{\tau}}\left(t+\tau\right)-\tau\right)\left( \overline{a}_{l}v+\overline{b}_{l}a_{l}\right)g\left(\mathbf{x}^{\prime},\mathbf{u}^ {\prime},t\right)\\ &+\left[t-\tau\left(1-e^{-\frac{t}{\tau}}\right)\right]\overline {A}g\left(\mathbf{x}^{\prime},\mathbf{u}^{\prime},t\right)\\ &+e^{-\frac{t}{\tau}}\left(\left(f_{0,jk}^{l}-\left(\frac{\partial f }{\partial x}\right)_{cf,jk}^{l}\cdot\mathbf{u}t-\left(\frac{\partial f}{\partial u }\right)_{cf,jk}^{l}\mathbf{u}t\right)H[\mathbf{u}\cdot\mathbf{n}]\right)\\ &+e^{-\frac{t}{\tau}}\left(\left(f_{0,jk}^{r}-\left(\frac{\partial f }{\partial x}\right)_{cf,jk}^{r}\cdot\mathbf{u}t-\left(\frac{\partial f}{\partial u }\right)_{cf,jk}^{r}\mathbf{u}t\right)\left(1-H[\mathbf{u}\cdot\mathbf{n}]\right)\right), \end{split}  \tag{15}\]  
where \(\mathbf{n}\) denotes the normal vector of interface. Finally, the macroscopic variables and discrete gas distribution function are updated respectively. For macroscopic conservative variables, we have  
\[\mathbf{W}^{n+1}=\mathbf{W}^{n}-\frac{1}{\Omega}\mathbf{F}\cdot\mathbf{S}+\mathbf{S}\mathbf{r}, \tag{16}\]  
where \(\Omega\) is volume of a cell and \(\mathbf{S}\) is the area of all the interfaces of the cell. The flux of macroscopic conservative variables reads  
\[\mathbf{F}=\int_{0}^{\Delta t}\int\mathbf{u}\cdot\mathbf{n}\mathbf{\psi}f_{cf}\left(\mathbf{x}, \mathbf{u},t\right)d\Xi dt. \tag{17}\]  
The source term \(\mathbf{S}\mathbf{r}\) is induced by inaccurate computation of reconstruction. The acceleration term does not appear in the evaluation of the particle trajectory when solving the gas distribution function at the cell interface. So a source term should be added in the projection stage [43], which reads  
\[\mathbf{S}\mathbf{r}=\left( \begin{array}{c}0\\ \rho a_{x}\Delta t\\ \rho a_{y}\Delta t\\ \widetilde{F}_{\rho}\cdot\overline{a}\end{array} \right). \tag{18}\]  
In above equations, integration under particle velocity space takes a discrete form. For discrete gas distribution function in each cell of physical space and particle velocity space, the update is realized by  
\[f_{i,j,k}^{n+1}=f_{i,j,k}^{n}-\frac{1}{\Omega}\int_{0}^{\Delta t}\mathbf{u}\cdot \mathbf{S}f_{cf}\left(\mathbf{x},\mathbf{u},t\right)dt+\frac{\Delta t}{2}\left(\frac{g_{ i,j,k}^{n+1}-f_{i,j,k}^{n+1}}{\tau^{n+1}}+\frac{g_{i,j,k}^{n}-f_{i,j,k}^{n}}{ \tau^{n}}\right). \tag{19}\]  
In all above formulas, the relaxation time \(\tau\) is related to different scales of particles in plasma, which reads [35]  
\[\tau_{\alpha}=\frac{1}{\sum\limits_{k}\frac{\theta_{\mathrm{xx}}\rho_{\mathrm{ x}}}{m_{\alpha}}}, \tag{20}\]  
where \(\rho_{\mathrm{x}}\) represents density of component \(\kappa\), \(m_{\kappa}\) represents particle mass of component \(\kappa\) and \(\theta_{\mathrm{xx}}\) represents collision frequency coefficient between component \(\alpha\) and \(\kappa\),  
\[\theta_{\mathrm{xx}}=\frac{4\sqrt{\pi}}{3}\left(\frac{2k_{B}T_{\alpha}}{m_{ \alpha}}+\frac{2k_{B}T_{\kappa}}{m_{\kappa}}\right)\left(\frac{d_{\alpha}+d_{ \kappa}}{2}\right)^{2}, \tag{21}\]  
where \(T\) represents temperature, \(k_{B}\) represents Boltzmann constant, and \(d\) represents diameter of particle. Next, we will introduce the process of solving the electromagnetic [PERSON]'s equation. Solving the [PERSON]'s equation amounts to finding the electric potential \(\Phi\) for a given charge distribution.