## 2 Models and Numerical Methods  
### Discharge Model and Numerical Method  
The three-equation drift-diffusion fluid model contains the mass-conservation equations, an energy-conservation equation and a Poisson's equation for the voltage. The electron's mass and energy conservation equations are  
\[\frac{\partial}{\partial t}(n_{e})+\
abla\cdot\vec{\Gamma}_{e}=R_{e}, \tag{1}\]  
\[\frac{\partial}{\partial t}(n_{e})+\
abla\cdot\vec{\Gamma}_{e}+\vec{E}\cdot \vec{\Gamma}_{e}=R_{e}, \tag{2}\]  
where \(e\) and \(\varepsilon\) are the electron and electron-energy, respectively. \(n_{e}\) and \(n_{\varepsilon}\) are the number density and the energy density of the electrons. \(R_{e}\) and \(R_{e}\) are the reaction induced source terms. \(\vec{\Gamma}_{e}\) and \(\vec{\Gamma}_{e}\) are the fluxes with drift-diffusion approximation, which are  
\[\vec{\Gamma}_{e}=-(\mu_{e}\vec{\varepsilon})n_{e}-\
abla(D_{e}n_{e})+n_{e} \vec{v}_{c}, \tag{3}\]  
\[\vec{\Gamma}_{e}=-(\mu_{e}\vec{\varepsilon})n_{\varepsilon}-\
abla(D_{e}n_{e} )+n_{e}\vec{v}_{c}. \tag{4}\]  
The drift and diffusion coefficients are calculated by using a Boltzmann equation solver Bolsig+ [32] with the electron energy probability function (EEPF) \(f_{0}\)[33],  
\[\mu_{e}=-\frac{\gamma}{3N}\int_{0}^{\infty}\frac{\varepsilon}{\sigma_{m}}\left( \frac{\partial f_{0}}{\partial\varepsilon}\right)d\varepsilon, \tag{5}\]\[D_{\varepsilon}=\frac{\gamma}{3N}\int_{0}^{\infty}\frac{\varepsilon}{\sigma_{m}}f_{ 0}d\varepsilon, \tag{6}\]  
\[\mu_{\varepsilon}=-\frac{\gamma}{3N\varepsilon}\int_{0}^{\infty}\frac{ \varepsilon^{2}}{\sigma_{m}}\left(\frac{\partial f_{0}}{\partial\varepsilon} \right)d\varepsilon, \tag{7}\]  
\[D_{\varepsilon}=\frac{\gamma}{3N\varepsilon}\int_{0}^{\infty}\frac{ \varepsilon^{2}}{\sigma_{m}}f_{0}d\varepsilon, \tag{8}\]  
where \(N\) is the total gas number density and \(m_{\varepsilon}\) is the mass of an electron. \(\sigma_{m}\) is the effective momentum transfer cross section. \(\vec{v}_{c}\) is the external gas velocity. \(R_{e}\) and \(R_{\varepsilon}\) are \(\sum\left(\left(\alpha-\eta\right)\left|\widetilde{\Gamma}_{\varepsilon} \right|+k_{r}\Pi n_{r}\right)\) and \(\sum\left(\left(\alpha-\eta\right)\left|\widetilde{\Gamma}_{\varepsilon} \right|\Delta\varepsilon+k_{r}\Delta r\Pi n_{r}\right)\). \(\alpha\) and \(\eta\) are the ionization and the attachment coefficients. \(k_{r}\) is the reaction rate and \(n_{r}\) is the reactants' number density. \(\Delta\varepsilon\) and \(\Delta r\) are the energy loss of electronic ionization reactions and other electron impact reactions. The heavy species' mass-conservation equations are  
\[\rho\frac{\partial}{\partial t}(\omega_{k})=\
abla\cdot\vec{\Gamma}_{k}+R_{k}, \tag{9}\]  
where \(\omega_{k}\) is the mass fraction of \(k\). \(R_{k}\) is the source term due to the chemical reactions. \(\vec{\Gamma}_{k}\) are the fluxes  
\[\vec{\Gamma}_{k}=\rho\omega_{k}\left(D_{k}\frac{\
abla\omega_{k}}{\omega_{k}} +D_{k}\frac{\
abla M}{M}-z_{k}\mu_{k}\vec{E}\right)+\rho\left(\omega_{k}\vec{ v}_{c}\right), \tag{10}\]  
where \(\mu_{k}\) and \(D_{k}\) indicate the drift and diffusion coefficients, respectively. \(z_{k}\) and \(M\) are the charge number and the average molar mass. The drift coefficients of ions are determined by the electric field [34]. The diffusion coefficients of ions are calculated from the generalized Einstein relation, and diffusion coefficients of the neutrals are calculated by the classical gas kinetic theory [35; 36]. \(\rho\) is the gas density, which is determined by the equation of state for a perfect gas,  
\[\frac{\partial}{\partial t}(n_{e})+\
abla\cdot\vec{\Gamma}_{e}+\vec{E}\cdot \vec{\Gamma}_{e}=R_{e}, \tag{2}\]0  
\(R\) is the universal gas constant. \(p\) and \(T\) are the pressure and temperature, respectively.