An asymptotic preserving well-balanced scheme for the isothermal fluid equations in low-temperature plasma applications

TL;DR

This paper introduces a novel asymptotic preserving well-balanced numerical scheme for low-temperature plasma fluid equations, accurately capturing plasma dynamics without resolving small scales like the Debye length.

Contribution

The scheme efficiently handles quasi-neutral regimes and sheath regions without implicit solvers, improving accuracy over standard explicit methods.

Findings

Accurately models plasma dynamics at low speeds and near electron thermal speed.

Preserves quasi-neutral limit without resolving Debye length.

Effectively captures plasma sheaths in low-temperature discharges.

Abstract

We present a novel numerical scheme for the efficient and accurate solution of the isothermal two-fluid (electron and ion) equations coupled to Poisson's equation for low-temperature plasmas. The model considers electrons and ions as separate fluids, comprising the electron inertia and charge separation. The discretization of this system with standard explicit schemes is constrained by very restrictive time steps and cell sizes related to the resolution of the Debye length, electron plasma frequency, and electron sound waves. Both sheath and electron inertia are fundamental to fully explain the physics in low-pressure and low-temperature plasmas. However, most of the phenomena of interest for fluid models occur at speeds much slower than the electron thermal speed and are quasi-neutral, except in small charged regions. In this work, we present a scheme based on the Lagrange-projection operator splitting that preserves the asymptotic regime where the plasma is quasi-neutral with massless electrons. As a result, the quasi-neutral regime is treated without the need of an implicit solver nor the resolution of the Debye length and electron plasma frequency. Additionally, the scheme proves to accurately represent the dynamics of the electrons both at low speeds and when the electron speed is comparable to the thermal speed. In addition, a well-balanced treatment of the ion source terms is proposed in order to tackle problems where the ion temperature is very low compared to the electron temperature. The scheme significantly improves the accuracy both in the quasi-neutral limit and in the presence of plasma sheaths when the Debye length is resolved. The scheme is assessed in two specifically designed test-cases: a quasi-neutral two-stream periodic perturbation with analytical solution and a low-temperature discharge that includes sheaths.

Figures

Tables

Table 1. Table 1: Characteristic values of an Argon RF discharge at 1 1 1 Pa [ 1 ] .

Dimensional quantities				Dimensionless quantities
Neutral density	$n_{𝔫}$	$1.25\times {10}^{20}$	${\text{m}}^{-3}$	Electron-to-ion mass ratio	$\epsilon =m_{𝔢}/m_{𝔦}$	$1.36\times {10}^{-5}$
Electron density	$n_{𝔢,𝔦}$	${10}^{16}$	${\text{m}}^{-3}$	Ion-to-electron temperature ration	$\kappa =T_{𝔦}/T_{𝔢}$	$0.025$
Neutral and ion temperature	$T_{𝔫,𝔦}$	$0.05$	eV	Normalized Debye length	$\lambda ={\lambda }_D/l$	$3.5\times {10}^{-3}$
Electron temperature	$T_{𝔢}$	$2$	eV	Ionization level	$n_{𝔢,𝔦}/n_{𝔫}$	$8\times {10}^{-5}$
Distance between plates	$l$	$3$	cm	Electron Knudsen number	${\text{Kn}}_{𝔢}$	$1.7$
Ion-neutral collisional cross section	${\sigma }_{𝔦𝔫}$	${10}^{-18}$	${\text{m}}^2$	Ion Knudsen number	${\text{Kn}}_{𝔦}$	$0.17$
Electron-neutral collisional cross section	${\sigma }_{𝔢𝔫}$	${10}^{-19}$	${\text{m}}^2$	Normalized ionization rate	${\overline{\nu }}^{iz}$	$0.0139$
Ionization constant	$K_{ion}$	$8.16\times {10}^{-18}$	${\text{m}}^3{\text{s}}^{-1}$	Normalized electron collision rate	${\overline{\nu }}_{𝔢}$	153.8
Ionization potential	${\epsilon }_{ion}$	$17.44$	eV	Normalized ion collision rate	${\overline{\nu }}_i$	$0.94$
Electron plasma period	${\omega }_{p𝔢}^{-1}$	$1.77\times {10}^{-10}$	s	Normalized plasma period	${\overline{\omega }}_{p𝔢}^{-1}={\omega }_{p𝔢}^{-1}{u}_B/l$	$1.29\times {10}^{-5}$

Table 2. Table 2: Solution of the two-stream periodic perturbation in a collisionless low-temperature plasma: Comparison of the computational performance with the standard and the AP Lagrange-projection scheme.

Numerical scheme and resolution	Cells	Iters. for one period	CPU time per iter. [s]	CPU time for one period [s]
Standard $\Delta x=\lambda$, $\Delta t=0.1{\omega }_{p𝔢}^{-1}$	$10000$	${10}^7$	$41.7175369$	$4.2\times {10}^8$
Standard $\Delta x=100\lambda$, $\Delta t=0.1{\omega }_{p𝔢}^{-1}$	$100$	${10}^7$	$1.1\times {10}^{-3}$	$1.1\times {10}^4$
AP $\Delta x=100\lambda$, $\Delta t=2500{\omega }_{p𝔢}^{-1}$	$100$	$400$	$1.1\times {10}^{-3}$	$0.44$

Equations

$${\partial }_t{n}_{\text{EuFrak}e}+{\partial }_x(n_{\text{EuFrak}e}{u}_{\text{EuFrak}e})$$

$${\partial }_t{n}_{\text{EuFrak}i}+{\partial }_x(n_{\text{EuFrak}i}{u}_{\text{EuFrak}i})$$

$$m_{\text{EuFrak}e}{\partial }_t(n_{\text{EuFrak}e}{u}_{\text{EuFrak}e})+{\partial }_x\left(m_{\text{EuFrak}e}{n}_{\text{EuFrak}e}{u}_{\text{EuFrak}e}^2+p_{\text{EuFrak}e}\right)$$

$$m_{\text{EuFrak}i}{\partial }_t(n_{\text{EuFrak}i}{u}_{\text{EuFrak}i})+{\partial }_x\left(m_{\text{EuFrak}i}{n}_{\text{EuFrak}i}{u}_{\text{EuFrak}i}^2+p_{\text{EuFrak}i}\right)$$

$${\partial }_{xx}^2\varphi$$

$${\partial }_{\bar{t}}{\bar{n}}_{\text{EuFrak}e}+{\partial }_{\bar{x}}({\bar{n}}_{\text{EuFrak}e}{\bar{u}}_{\text{EuFrak}e})$$

$${\partial }_{\bar{t}}{\bar{n}}_{\text{EuFrak}i}+{\partial }_{\bar{x}}({\bar{n}}_{\text{EuFrak}i}{\bar{u}}_{\text{EuFrak}i})$$

$${\partial }_{\bar{t}}({\bar{n}}_{\text{EuFrak}e}{\bar{u}}_{\text{EuFrak}e})+{\partial }_{\bar{x}}\left[{\bar{n}}_{\text{EuFrak}e}\left({\bar{u}}_{\text{EuFrak}e}^2+{\epsilon }^{-1}\right)\right]$$

$${\partial }_{\bar{t}}({\bar{n}}_{\text{EuFrak}i}{\bar{u}}_i)+{\partial }_{\bar{x}}\left[{\bar{n}}_{\text{EuFrak}i}\left({\bar{u}}_{\text{EuFrak}i}^2+\kappa \right)\right]$$

$${\partial }_{\bar{x}\bar{x}}^2\bar{\varphi }$$

$$\text{Electron-ion mass ratio: }\epsilon =m_{\text{EuFrak}e}/m_{\text{EuFrak}i},$$

$$\text{squared non-dimensional Debye length: }\chi =\frac{{ϵ}_0{\text{k}}_{\text{B}}{T}_{\text{EuFrak}e}}{n_{\text{EuFrak}e}{e}^2{L}_0^2},$$

$${\bar{n}}_{\text{EuFrak}e}$$

$${\partial }_{\bar{t}}{\bar{n}}_{\text{EuFrak}i}+{\partial }_{\bar{x}}({\bar{n}}_{\text{EuFrak}i}{\bar{u}}_{\text{EuFrak}i})$$

$${\partial }_{\bar{x}}({\bar{n}}_{\text{EuFrak}e}{\bar{u}}_{\text{EuFrak}e})$$

$${\partial }_{\bar{t}}({\bar{n}}_{\text{EuFrak}i}{\bar{u}}_{\text{EuFrak}i})+{\partial }_{\bar{x}}({\bar{n}}_{\text{EuFrak}i}{\bar{u}}_{\text{EuFrak}i}^2+\kappa {\bar{n}}_{\text{EuFrak}i})$$

$${\partial }_{\bar{x}}\bar{\varphi }$$

$$f(x,t)=f^{(0,\chi )}+\epsilon f^{(1,\chi )}+\mathcal{O}({\epsilon }^2).$$

$${\partial }_{\bar{t}}{\bar{n}}_{\text{EuFrak}e}^{(0,\chi )}+{\partial }_{\bar{x}}({\bar{n}}_{\text{EuFrak}e}^{(0,\chi )}{\bar{n}}_{\text{EuFrak}i}^{(0,\chi )})$$

$${\partial }_{\bar{t}}{\bar{n}}_{\text{EuFrak}i}^{(0,\chi )}+{\partial }_{\bar{x}}({\bar{n}}_{\text{EuFrak}i}^{(0,\chi )}{\bar{u}}_{\text{EuFrak}i}^{(0,\chi )})$$

$$\frac{1}{{\bar{n}}_{\text{EuFrak}e}^{(0,\chi )}}{\partial }_{\bar{x}}{\bar{n}}_{\text{EuFrak}e}^{(0,\chi )}$$

$${\partial }_{\bar{t}}\left({\bar{n}}_{\text{EuFrak}i}^{(0,\chi )}{\bar{u}}_{\text{EuFrak}i}^{(0,\chi )}\right)+{\partial }_{\bar{x}}\left[{\bar{n}}_{\text{EuFrak}i}^{(0,\chi )}\left({\bar{u}}_{\text{EuFrak}i}^{(0,\chi {)}^2}+\kappa \right)\right]$$

$${\partial }_{\bar{x}\bar{x}}^2{\varphi }^{(0,\chi )}$$

$$f(x,t){=}^{(\epsilon ,0)}f+\chi {}^{(\epsilon ,1)}f+\mathcal{O}({\chi }^2).$$

$${\partial }_{\bar{t}}^{(\epsilon ,0)}{\bar{n}}_{\text{EuFrak}e}+{\partial }_{\bar{x}}\left({}^{(\epsilon ,0)}{\bar{n}}_{\text{EuFrak}e}^{(\epsilon ,0)}{\bar{u}}_{\text{EuFrak}e}\right)$$

$${\partial }_{\bar{t}}^{(\epsilon ,0)}{\bar{n}}_{\text{EuFrak}i}+{\partial }_{\bar{x}}\left({}^{(\epsilon ,0)}{\bar{n}}_{\text{EuFrak}i}^{(\epsilon ,0)}{\bar{u}}_{\text{EuFrak}i}\right)$$

$${\partial }_{\bar{t}}\left({}^{(\epsilon ,0)}{\bar{n}}_{\text{EuFrak}e}^{(\epsilon ,0)}{\bar{u}}_{\text{EuFrak}e}\right)+{\partial }_{\bar{x}}\left[{}^{(\epsilon ,0)}{\bar{n}}_{\text{EuFrak}e}\left({}^{(\epsilon ,0)}{\bar{u}}_{\text{EuFrak}e}^2+{\epsilon }^{-1}\right)\right]$$

$${\partial }_{\bar{t}}\left({}^{(\epsilon ,0)}{\bar{n}}_{\text{EuFrak}i}^{(\epsilon ,0)}{\bar{u}}_{\text{EuFrak}i}\right)+{\partial }_{\bar{x}}\left[{}^{(\epsilon ,0)}{\bar{n}}_{\text{EuFrak}i}\left({}^{(\epsilon ,0)}{\bar{u}}_{\text{EuFrak}i}^2+\kappa \right)\right]$$

$${}^{(\epsilon ,0)}{\bar{n}}_{\text{EuFrak}e}$$

$${\partial }_{\bar{t}}{\bar{n}}_{\text{EuFrak}e}^{(0,0)}+{\partial }_{\bar{x}}({\bar{n}}_{\text{EuFrak}e}^{(0,0)}{\bar{n}}_{\text{EuFrak}i}^{(0,0)})$$

$${\partial }_{\bar{t}}{\bar{n}}_{\text{EuFrak}i}^{(0,0)}+{\partial }_{\bar{x}}({\bar{n}}_{\text{EuFrak}i}^{(0,0)}{\bar{u}}_{\text{EuFrak}i}^{(0,0)})$$

$$\frac{1}{{\bar{n}}_{\text{EuFrak}e}^{(0,0)}}{\partial }_{\bar{x}}{\bar{n}}_{\text{EuFrak}e}^{(0,0)}$$

$${\partial }_{\bar{t}}\left({\bar{n}}_{\text{EuFrak}i}^{(0,0)}{\bar{u}}_{\text{EuFrak}i}^{(0,0)}\right)+{\partial }_{\bar{x}}\left[{\bar{n}}_{\text{EuFrak}i}^{(0,0)}\left({\bar{u}}_{\text{EuFrak}i}^{(0,0{)}^2}+\kappa \right)\right]$$

$${\bar{n}}_{\text{EuFrak}e}^{(0,0)}$$

$${\partial }_{\bar{t}}^{(0,0)}{\bar{n}}_{\text{EuFrak}e}+{\partial }_{\bar{x}}\left({}^{(0,0)}{\bar{n}}_{\text{EuFrak}e}^{(0,0)}{\bar{u}}_{\text{EuFrak}e}\right)$$

$${\partial }_{\bar{t}}^{(0,0)}{\bar{n}}_{\text{EuFrak}i}+{\partial }_{\bar{x}}\left({}^{(0,0)}{\bar{n}}_{\text{EuFrak}i}^{(0,0)}{\bar{u}}_{\text{EuFrak}i}\right)$$

$$\frac{1}{{}^{(0,0)}{\bar{n}}_{\text{EuFrak}e}}{\partial }_{\bar{x}}^{(0,0)}{\bar{n}}_{\text{EuFrak}e}$$

$${\partial }_{\bar{t}}\left({}^{(0,0)}{\bar{n}}_{\text{EuFrak}i}^{(0,0)}{\bar{u}}_{\text{EuFrak}i}\right)+{\partial }_{\bar{x}}\left[{}^{(0,0)}{\bar{n}}_{\text{EuFrak}i}\left({}^{(0,0)}{\bar{u}}_{\text{EuFrak}i}^2+\kappa \right)\right]$$

$${}^{(0,0)}{\bar{n}}_{\text{EuFrak}e}$$

$$\text{Dimensional: }{n}_{\text{EuFrak}e}{u}_{\text{EuFrak}e}{|}_{wall}=n_{\text{EuFrak}e}\sqrt{\frac{{\text{k}}_{\text{B}}{T}_{\text{EuFrak}e}}{2\pi m_{\text{EuFrak}e}}}\text{and}\text{dimensionless: }{\bar{n}}_{\text{EuFrak}e}{\bar{u}}_{\text{EuFrak}e}{|}_{wall}=\frac{{\bar{n}}_{\text{EuFrak}e}}{\sqrt{2\pi \epsilon }}.$$

$$2n_{\text{EuFrak}e}{u}_{\text{EuFrak}e}{|}_{wall}={\int }_0^l{n}_{\text{EuFrak}e}{\nu }^{iz}dx.$$

$${\varphi }_p=-\frac{{\text{k}}_{\text{B}}{T}_{\text{EuFrak}e}}{2e}\text{and}{\varphi }_W=\frac{{\text{k}}_{\text{B}}{T}_{\text{EuFrak}e}}{e}\ln {\left(\frac{m_{\text{EuFrak}e}}{2\pi m_i}\right)}^{1/2},$$

$$\left(\begin{array}[]{c}\bar{n}_{\EuFrak{e}}\\ \bar{n}_{\EuFrak{i}}\\ \bar{n}_{\EuFrak{e}}\bar{u}_{\EuFrak{e}}\\ \bar{n}_{\EuFrak{i}}\bar{u}_{\EuFrak{i}}\\ \end{array}\right)_{j}^{{n+1}}=\left(\begin{array}[]{c}\bar{n}_{\EuFrak{e}}\\ \bar{n}_{\EuFrak{i}}\\ \bar{n}_{\EuFrak{e}}\bar{u}_{\EuFrak{e}}\\ \bar{n}_{\EuFrak{i}}\bar{u}_{\EuFrak{i}}\\ \end{array}\right)_{j}^{n}-\frac{\Delta t}{\Delta x}\left[\left(\begin{array}[]{c}\bar{n}_{\EuFrak{e}}\bar{u}_{\EuFrak{e}}\\ \bar{n}_{\EuFrak{i}}\bar{u}_{i}\\ \bar{n}_{\EuFrak{e}}(\bar{u}_{\EuFrak{e}}^{2}+\varepsilon^{-1})\\ \bar{n}_{\EuFrak{i}}(\bar{u}_{i}^{2}+\kappa)\\ \end{array}\right)_{j+1/2}^{n}-\left(\begin{array}[]{c}\bar{n}_{\EuFrak{e}}\bar{u}_{\EuFrak{e}}\\ \bar{n}_{\EuFrak{i}}\bar{u}_{i}\\ \bar{n}_{\EuFrak{e}}(\bar{u}_{\EuFrak{e}}^{2}+\varepsilon^{-1})\\ \bar{n}_{\EuFrak{i}}(\bar{u}_{i}^{2}+\kappa)\\ \end{array}\right)_{j-1/2}^{n}\right]+\Delta t\left(\begin{array}[]{c}(\bar{n}_{\EuFrak{e}}\bar{\nu}^{iz})\\ \bar{n}_{\EuFrak{e}}\bar{\nu}^{iz}\\ \bar{n}_{\EuFrak{e}}\varepsilon^{-1}\partial_{\bar{x}}\bar{\phi}\\ -\bar{n}_{\EuFrak{i}}\partial_{\bar{x}}\bar{\phi}\\ \end{array}\right)_{j}^{n}$$

$${\bar{\varphi }}_{j+1}^{n+1}-2{\bar{\varphi }}_j^{n+1}+{\bar{\varphi }}_{j-1}^{n+1}=\frac{\Delta x^2}{\chi }\left({\bar{n}}_{\text{EuFrak}{e}_j}^n-{\bar{n}}_{\text{EuFrak}{i}_j}^n\right)\text{and}{\partial }_{\bar{x}}{\bar{\varphi }}_j^n=\frac{1}{2\Delta x}\left({\bar{\varphi }}_{j+1}^n-{\bar{\varphi }}_{j-1}^n\right).$$

$${\text{CFL}}^{conv}=\frac{\Delta t|{\lambda }_{max}^{\text{EuFrak}e,\text{EuFrak}i}|}{\Delta x}\text{with}|{\lambda }_{max}^{\text{EuFrak}e,\text{EuFrak}i}|=\max \left(|{\bar{u}}_{\text{EuFrak}e,\text{EuFrak}i}+{\bar{c}}_{\text{EuFrak}e,\text{EuFrak}i}|,|{\bar{u}}_{\text{EuFrak}e,\text{EuFrak}i}-{\bar{c}}_{\text{EuFrak}e,\text{EuFrak}i}|\right).$$

Full text

1 Introduction

Low-temperature plasmas are fundamental to a wide range of technological applications – from lighting to semiconductor manufacturing and electric propulsion. These plasmas are created by an electrical discharge in a gas. The physical conditions where they occur are particularly diverse: pressures varying from below a millitorr to few hundred atmospheres, different excitation methods (e.g., inductively coupled plasma (ICP), capacitively coupled plasma (CCP), dielectric barrier discharge (DBD), magnetron), diverse geometries, and power ranges (see, e.g., [1, 2]). For that reason, a number of different descriptions are used for their study, from kinetic approaches, e.g., Particle-In-Cell/Monte Carlo Collision (PIC-MCC) and Direct Simulation Monte Carlo (DSMC) method, to fluid and global models (see [3, 4] for a review on the modeling and numerical approaches of low-temperature plasmas).

In atmospheric pressure discharges, the drift-diffusion approximation (e.g., [5, 6, 7, 8, 9, 10, 11]) is a valid simplification. However, under low-pressure conditions, the momentum transfer between species is smaller since the number of collisions decreases for decreasing pressures. This results in non-equilibrium conditions that may lead to a differential motion between the heavy species and electrons. This differential motion is responsible for plasma instabilities such as the two-stream instability [12] or the electron drift instability that is observed to cause anomalous transport in Hall thrusters [13, 14, 15, 16, 17, 18, 19]. Even though the electron inertial term is small in the electron momentum equations, these instabilities cannot be explained without it. Similarly, electron inertia cannot be completely neglected to fully explain the physics of sheaths [20], in electron sheaths and presheaths [21, 22], in plasma sheath instabilities [23], Langmuir probes [24], in magnetized plasmas [25, 26, 27], and Hall thrusters [28, 29].

For that reason, in low-pressure discharges, alternative models to the drift-diffusion approach such as kinetic and hybrid models, the plasma moment equations or the multicomponent non-equilibrium one velocity fluid model [30, 31] are needed. Kinetic simulations provide a very accurate description of the state of the plasma, but are computationally very expensive. On the other hand, hybrid models [32, 33, 34] are a cheaper alternative that combines the fluid and the kinetic description. However, both methods need to resolve the Debye length and the electron plasma frequency if the time integration is explicit and the quasi-neutrality hypothesis is not assumed.

In Table 1, we present characteristic conditions of a typical low-temperature low-pressure Argon plasma RF discharge (from [1]). We first note that under these pressure conditions, the Knudsen number is of order one for electrons and ten times smaller for ions. The kinetic phenomena is important and transport models are needed to provide a closure for the plasma fluid equations. The second important feature of Table 1 is related to the smallness of the electron-to-ion mass ratio and the normalized Debye length. Most of the explicit methods (PIC methods, hybrid methods or deterministic discretization of the multi-fluid equations) require to resolve these scales in order to guarantee the stability of the scheme. However, if we consider an explicit 1D simulation of the discharge of Table 1 that resolves the Debye length and the plasma frequency with ten spatial and temporal points, the numerical set-up would need around $10^{4}$ spatial points (in 1D) and $10^{6}$ time steps to simulate one transit time of an ion acoustic wave. Therefore, this set-up involves a significant computational cost for a 1D simulation. Implicit methods can guarantee the numerical stability with larger time-steps. However, the Newton method to solve the linear system can become also very expensive.

In the present work, we consider the isothermal plasma fluid equations under collisionless conditions. This system of equations contains the numerical difficulties associated to the low-temperature moment plasma equations in electrostatic conditions. The model considers the first two moments of Vlasov’s equation for electrons and ions while assuming both temperatures to be constant. The dynamics of the two fluids are coupled through the electric potential, which satisfies a Poisson equation. These equations are widely used in the sheath theory [35, 36, 37] and in the study of plasma waves [12]. Although kinetic phenomena play an important role in the plasma sheath [21], a fluid model that does not assume quasi-neutrality can potentially capture the interaction between the macroscopic scales and the sheaths.

The most general form of the multi-fluid plasma equations considers the electron inertial term in the electron momentum equation and charge separation effects. This is equivalent to considering a finite electron mass and a finite Debye length. As mentioned before, this allows for representing plasma instabilities such as the two-stream or the drift waves and charged regions of the plasma such as the sheaths. Nevertheless, both the electron mass and the Debye length are two very small parameters as compared to the ion and macroscopic scales. Consequently, these small scales impose very restrictive numerical constraints related to the resolution of the Debye length, electron plasma waves, and electron acoustic waves. Consequently, due to these requirements, the time-steps and mesh sizes are not significantly more advantageous than in the kinetic or hybrid approach. Moreover, the multiscale character of the problem can lead to very large discretization errors when the small scales are not properly resolved (see [38]).

Numerical methods for the ideal multi-fluid equations coupled to Maxwell’s equations have been proposed by a number of authors [39, 40, 41, 42, 43, 44, 45, 46], for the study of plasma sheaths [47], and for the study of plasma expansion in vacuum with the isentropic electrostatic approximation [48, 49, 50, 51, 52]. The main difficulties of the low-temperature multi-fluid plasma model and the solutions proposed in the current literature are summarized in the following. First, the stability of explicit schemes for the multi-fluid equations solving for the Poisson equation requires the resolution of the electron plasma wave frequency [52]. Alvarez Laguna et al. [44, 45] have proposed implicit time integration for the previous set of equations coupled to full Maxwell’s equations. However, the inversion of the matrix is computationally costly, which can be improved by the use of GPUs [46]. Alternatively, a more advantageous approach is the asymptotic-preserving (AP) scheme proposed such as the one proposed by Degond et al. [50, 51, 52]. AP schemes [53] preserve the quasi-neutral asymptotic limit without the resolution of the Debye length nor the plasma frequency. These schemes do not need an implicit solver, which results in a significant advantage from the computational point of view. However, the AP scheme from [50, 51] is not able to tackle the asymptotic limit of the small electron-to-ion mass ratio, finding large errors in the velocity of the electrons even with high order schemes. Additionally, the AP scheme proposed in [50] demands the solution of a second order equation for the electric potential that requires the storage of the solution in two different time-steps, which increases the computational stencil.

The second difficulty of the multi-fluid equations is due to the small mass of the electrons, which results in a very large electron speed of sound. However, in most of the cases, the electron fluid travels at much smaller speeds than the electron thermal speed. This corresponds to a low-Mach regime for electrons, which is known to give numerical problems due to an excessive numerical dissipation that restricts the time step of compressible solvers [54]. Different strategies are proposed in order to tackle the low-Mach regime in compressible solvers. One method are the so-called all-regime flux-splitting methods such as the AUSM+-up [55] (applied to the multi-fluid equations in [44, 45, 46]) or the preconditioning methods to remove the stiffness of the low-Mach regime [56]. Similarly, the operator-splitting Lagrange-projection scheme is combined to the preconditioning method by Chalons et al. [57]. Additionally, the AP schemes are also used to tackle the incompressible asymptotic behaviour in the Euler equations [58].

The third difficulty is related to the regime where the temperature of the ions is much lower than the electron temperature. In low-temperature plasmas, the electric potential, in general, scales as the electron temperature (in eV). This results in a Lorentz force that is much larger than the ion pressure flux. The main problem arises when the ion convective flux is treated with an upwind scheme and the Lorentz force (which involves the gradient of the potential) with a cell-centered scheme. The upwind scheme provides accurate non-oscillatory solutions in homogeneous problems, but they can lose accuracy in the presence of stiff source terms [59]. A solution to this problem is to use well-balanced schemes that upwind the source terms in a consistent manner [60, 61, 62].

In the present work, we propose a numerical scheme that addresses these three problems and design three dedicated discriminating test-cases in order to benchmark such a method. Regarding the stiffness introduced by the smallness of the Debye length and the electron mass, we propose a novel operator splitting method based on the all-regime Lagrange-projection [57] coupled to the Poisson equation. We retrieve a numerical scheme that has the AP property of the quasi-neutral limit with massless electrons. Therefore, neither the Debye length nor the electron plasma frequency need to be resolved to be stable and accurate, without the need of an implicit solver. As compared to previous AP schemes for the multi-fluid equations [50, 51], our approach does not need to solve an equivalent second order equation for the electric potential, but it solves the standard Poisson equation. This results in a much simpler algorithm for the potential. Additionally, the main advantage is that the numerical scheme can simultaneously tackle the problem of the small mass of electrons. The solution of these two problems with a unique AP scheme is an original contribution of this work. This can be considered a step forward that could help to reduce the numerical cost of solving the electron fluid dynamics coupled to the one of ions. Finally, a well-balanced discretization of the Lorentz force in the ion momentum equation is proposed. In the results we show that an incorrect discretization of this term can lead to spurious numerical instabilities, which can be avoided with well-balanced schemes.

The numerical scheme is then benchmarked against three numerical set-ups that allows us to assess our strategy in terms of accuracy and computational cost. The first one simulates a quasi-neutral periodic perturbation in a thermal plasma. This case is used to test the asymptotic preserving property of the discretization for a small Debye length and a small electron-to-ion mass ratio. The same case is reproduced in a low-temperature plasma in order to assess the well-balanced discretization of the ions. We demonstrate that the proposed numerical strategy allows for a dramatic reduction of the computational time in cases with quasi-neutral plasma, as compared to a standard discretization. Finally, a low-temperature plasma discharge between two floating walls is simulated. This realistic set-up is able to capture the physics of the electrically charged plasma sheath coexisting with a quasi-neutral bulk and encompasses most of the difficulties of more realistic configurations, while amenable to detailed analysis of the proposed methods.

Outline of the paper

First, we present the set of equations, its normalization, and the asymptotic study when both the mass of electrons and the Debye length are very small as compared to the ion scales. Second, we present a standard discretization of the system in order to illustrate the difficulties of the numerical discretization. Third, we describe the operator splitting scheme and the well-balanced treatment of the ion source-terms while proving that the acoustic step preserves the asymptotic behaviour. Fourth, we simulate a two-stream perturbation in thermal and low-temperature plasmas. Finally, we self-consistently simulate a low-temperature discharge between two floating plates. The numerical scheme proves to be accurate both in the quasi-neutral limit and when the sheaths are included in the computational domain.

2 Set of equations and asymptotic behaviour

We consider the isothermal plasma (electron + ion) equations. The equations are obtained by taking the first two moments of the kinetic equation for electrons and ions (see, e.g., [63] for a derivation). These moments correspond to the mass and momentum balance laws for the two charged species. The system considers ionization reactions and, in the present paper, we neglect the effect of recombination and the elastic collisions. As explained in [1], in most discharges these scales are much smaller than the ionization scale. Nevertheless, the elastic collisions in low-pressure conditions do not impose a major difficulty from the numerical point of view and for that reason are not the focus of the present work. Finally, the system is closed with the Poisson equation for the electric potential.

The equations in dimensional form read

[TABLE]

where $n_{\EuFrak{e}}$ and $n_{\EuFrak{i}}$ stand for the electron and ion number density respectively, and $u_{\EuFrak{e}}$ and $u_{\EuFrak{i}}$ the electron and ion velocities. The electron-impact ionization rate coefficient $\nu^{iz}$ is a function of the electron temperature $\nu^{iz}=n_{\EuFrak{n}}K_{ion}(T_{\EuFrak{e}})$, for instance through Arrhenius’ law, $K_{ion}(T_{\EuFrak{e}})=A\exp[{{-\varepsilon_{ion}}/{(\textup{k}_{\textup{B}}T_{\EuFrak{e}})}}]$, where the quantity $\varepsilon_{ion}$ is the ionization energy. The neutral number density $n_{\EuFrak{n}}$ is assumed to be constant. The partial pressures of the electron and ion fluids are assumed to obey the perfect gas law, $p_{\EuFrak{e}}=n_{\EuFrak{e}}\textup{k}_{\textup{B}}T_{\EuFrak{e}}$ and $p_{\EuFrak{i}}=n_{\EuFrak{i}}\textup{k}_{\textup{B}}T_{\EuFrak{i}}$, where $\textup{k}_{\textup{B}}$ is Boltzmann’s constant, and $T_{\EuFrak{e}}$ and $T_{\EuFrak{i}}$, the electron and ion constant temperatures, respectively. The plasma is considered to be in thermal non-equilibrium, i.e., $T_{\EuFrak{e}}\neq T_{\EuFrak{i}}$, at constant temperatures.

Despite the isothermal assumption, the considered set of equations contains the main difficulties of the low-temperature moment plasma models. These difficulties, as described in the introduction, are related to the mass disparity between electrons and ions, the small Debye length, and the low-temperature of the ions.

2.1 Normalized equations

The set of equations (1a)-(1d) and (1e) are normalized by introducing some reference quantities: $n_{0}$, the characteristic number density common to electrons and ions, $L_{0}=l$, the reference length, $T_{\EuFrak{e}}$, the electron temperature, and $T_{\EuFrak{i}}$, the ion temperature. The rest of the reference variables are calculated as a combination of the previous ones: the reference velocity common to electrons and ions is based on the Bohm velocity $u_{0}\equiv u_{B}=\sqrt{\textup{k}_{\textup{B}}T_{\EuFrak{e}}/m_{\EuFrak{i}}}$, the charateristic time $t_{0}=L_{0}/u_{0}$ is obtained from the reference velocity and reference distance, whereas the reference potential $\phi_{0}=\textup{k}_{\textup{B}}T_{\EuFrak{e}}/e$ is based on the thermal energy of electrons. The normalized set of equations reads

[TABLE]

where the non-dimensional parameters are defined as:

[TABLE]

Note that, for the sake of simplicity of notation, we use the square of the Debye length $\chi$ as a non-dimensional parameter. Additionally, we define the non-dimensional Debye length $\lambda=\sqrt{\chi}$ as it will be used to describe characteristic lengths of the problem.

We highlight the importance of the three nondimensional parameters $\varepsilon$, $\chi$, and $\kappa$ in order to build a numerical scheme. In Table 1, we present the typical values in an Argon RF discharge. As discussed in the introduction, the smallness of the $\varepsilon$ and $\chi$ impose very restrictive numerical constraints. In the following section, we study the asymptotic behaviour when these two parameters tend to zero.

2.2 Asymptotic behavior

We study the multiscale asymptotic behavior [64] with respect to $\varepsilon$ and $\chi$. Previous work [52, 50] performed this study only with respect to the Debye length, $\lambda$. However, the inclusion of $\varepsilon$ in the analysis is fundamental as the electron velocity is generally much smaller than the thermal speed of electrons. Since $\varepsilon$ and $\chi$ are the smallest parameters of the system, we do not include $\kappa$ or $\bar{\nu}^{iz}$ in our study.

From a physics point of view, we can consider three different asymptotic behaviours that correspond to different plasma phenomena, as illustrated in Fig. 1. The complete problem ${}^{\chi}F^{\varepsilon}$ corresponds to the system of eqs. (2). This system considers finite Debye length and electron inertia. This problem resolves all the possible scales, being the fastest one corresponding to the electron plasma waves. The main problem of designing a numerical scheme for these small scales is that it might be very inefficient to represent the macroscopic scales and leading to consistency problems due to an imbalanced numerical dissipation when the two parameters are small.

The regime where the Debye length tends to zero for arbitrarily small $\varepsilon$ corresponds to a quasineutral plasma with the electrons that can move at bulk speed closer to the electron thermal velocity. This problem is of interest specially in the presence of a magnetic field that can produce drift motions at very high speed, such as in Hall effect thrusters. This regime allows for the representation of plasma instabilities such as the two-stream instability or the electron-drift instability. We denote this problem as ${}^{0}F^{\varepsilon}$.

Alternatively, we can consider the asymptotic behaviour of $\varepsilon$ tending to zero for a finite Debye length. In this regime, the electrons move at speeds comparable to the ion sound velocity (Bohm’s velocity) and the Debye length is arbitrarily small. This regime is important, for instance, in the plasma-sheath transition. We denote this problem as ${}^{\chi}F^{0}$.

Finally, we consider the asymptotic limit where the electrons travel at speeds comparable to the ion velocity in a quasi-neutral plasma, i.e., $\varepsilon\rightarrow 0$ and $\chi\rightarrow 0$ This behaviour is present in most of the phenomena occuring at ion scales in cold and thermal plasmas. For that reason, we focus in this problem in the rest of the paper. We denote the problem as ${}^{0}F^{0}$. In the following, we show the set of equations corresponding to the problem ${}^{0}F^{0}$ and we demonstrate that $\lim\limits_{\chi\rightarrow 0}~{}^{\chi}F^{0}=\lim\limits_{\varepsilon\rightarrow 0}~{}^{0}F^{\varepsilon}\equiv~{}^{0}F^{0}$ in the case of periodic boundary conditions.

Let us consider the system of equations (2a)-(2e) with periodic boundary conditions, which we denote as the problem ${}^{\chi}F^{\varepsilon}$. We consider two different asymptotic expansions for the problem ${}^{\chi}F^{\varepsilon}$: (1) In terms of the small parameter $\varepsilon$ and (2) in terms of the small parameter $\chi$. We define the problem ${}^{\chi}F^{0}$ as the system of equations for the zero-th order terms of the expansion in terms of $\varepsilon$. Similarly, we define ${}^{0}F^{\varepsilon}$ as the the system of equations for the zero-th order terms of the expansion in terms of $\chi$.

Proposition 1**.**

The system of equations that corresponds to $\lim\limits_{\chi\rightarrow 0}~{}^{\chi}F^{0}$ with periodic boundary conditions is formally the same as $\lim\limits_{\varepsilon\rightarrow 0}~{}^{0}F^{\varepsilon}$ and it reads

[TABLE]

Proof.

In order to prove Proposition 1, we first propose the following expansion for the variables of the problem ${}^{\chi}F^{\varepsilon}$

[TABLE]

By injecting this expansion in the system of equations (2a)-(2e), we find that the system of equations for the zero-th order terms reads

[TABLE]

We denote this system as $~{}^{\chi}F^{0}$.

Alternatively, we propose an expansion for the variables of the problem ${}^{\chi}F^{\varepsilon}$ in terms of $\chi$ of the form

[TABLE]

The system of equations corresponding to the zero-th order terms reads

[TABLE]

We denote this system as ${}^{0}F^{\varepsilon}$.

We can use the expansion of eq. (7) for the problem ${}^{\chi}F^{0}$ in order to study the problem $F^{(0,0)}\equiv\lim\limits_{\chi\rightarrow 0}~{}^{\chi}F^{0}$. The system of equations for the zero-th order terms is very similar to system ${}^{\chi}F^{0}$, with the difference in the Poisson equation, as shown below

[TABLE]

By using the condition (9e) in eqs (9a) and (9b), we obtain that the fluxes are conserved, i.e., $\partial_{\bar{x}}(\bar{n}_{\EuFrak{e}}^{(0,0)}\bar{n}_{\EuFrak{i}}^{(0,0)})=\partial_{\bar{x}}(\bar{n}_{\EuFrak{i}}^{(0,0)}\bar{u}_{\EuFrak{i}}^{(0,0)})$. This proves that the system $F^{(0,0)}$ is the system proposed as ${}^{0}F^{0}$ in eqs. (4).

Alternatively, we use the expansion of eq. (5) for the problem ${}^{0}F^{\varepsilon}$. The equations for the zero-th order terms are denoted as ${}^{(0,0)}F\equiv\lim\limits_{\varepsilon\rightarrow 0}~{}^{0}F^{\varepsilon}$. In this new system of equations, the electron momentum equation is modified as compared to the one in ${}^{0}F^{\varepsilon}$.

[TABLE]

This system of equations is the same as the problem $F^{(0,0)}$. Consequently, in the case of periodic boundary conditions, both limits result in the same system of equations (4a)-(4e), i.e., ${}^{(0,0)}F=F^{(0,0)}=~{}^{0}F^{0}$.

∎

We highlight that the system $~{}^{0}F^{\varepsilon}$, which corresponds to the limit of the Debye length tending to zero, is the one that was previously considered by the works [52, 50]. In this system, the electric potential becomes a Lagrange multiplier that imposes the charge neutrality. In system (8), we decided to write the set of equations in this form for clarity. Nevertheless, a differential equation for the electric potential can be retrieved from the conservation of electric charge, as done in [50].

The asymptotic behaviour proposed in this paper considers massless electrons as in [65, 66]. It is therefore very different to the one derived in [52, 50]. By including in our analysis the electron-to-ion mass ratio, we find that the quasi-neutrality is found in eq. (4a) and the electron density follows the Boltzmann distribution in eq. (4e). Furthermore, the electric current is conserved by equation (4d) and the ion mass and momentum equations are unchanged. The main difference as compared to [50] is the limit for the zero-th order of the potential in eq. (4e) and the electron momentum equation.

2.3 Study of a bounded low-temperature plasma through the fluid plasma equations

In low-temperature plasma industrial applications, plasmas are confined. As a result, the charged particles that are produced inside the reactor through ionization, are lost through the boundaries when they strike the wall. Since the thermal motion of electrons is larger than that of ions, the surface will charge negatively with respect to the plasma (in electropositive plasmas), forming a charged boundary layer called the plasma sheath.

The analytical models for the sheath and presheath rely on the isothermal multi-fluid equations [67], while assuming the inertia of electrons and temperature of ions to be negligible. Nevertheless, important kinetic phenomena taking place are not included in this model [68].

Let us consider a 1D domain of length $l$ filled with a plasma between two floating walls, with no secondary electron emission, the distribution function of electrons is a Maxwellian and all the electrons that touch the wall are absorbed by the wall. With these assumptions, the flux of electrons collected by the wall (see, e.g., [1]) both in dimensional and dimensionless units read:

[TABLE]

A steady solution is found when the ionization inside the bulk of the plasma balances the particle loss as follows

[TABLE]

As mentioned by Riemann [37], the ionization frequency is an eigenvalue of the problem. Consequently, there is only one ionization frequency that finds a steady state solution for a given distance between plates. In this paper, we propose a numerical methodology that proves to be convergent to find this eigenvalue.

With the previous assumptions, the potential at the pre-sheath $\phi_{p}$ and the wall $\phi_{W}$, in dimensional units [69], as follows

[TABLE]

where $\phi_{p}$ is the potential drop needed to accelerate the ions to Bohm’s speed (neglecting the ion pressure gradient) and $\phi_{W}$ is the potential drop in the sheath. In Fig. 2, we illustrate the steady state solution of a bounded plasma between two floating plates.

3 Standard upwind finite volume discretization

We present a standard discretization of the system (2) in order to illustrate the associated numerical difficulties. An example of a simulation of a low-temperature discharge with this discretization can be found in Alvarez Laguna et al. [38]. Alternatively, a similar discretization is described in [50] in order to illustrate a standard solver of the Euler-Poisson system.

We use a finite volume discretization where the domain $x\in[0,l]$ is divided into $N$ elements of equal length $\Delta x$. We approximate the value of the unknowns as a piecewise function inside the volume $\Omega_{j}$. The flux at the interfaces is approximated by a numerical flux function that is in general a function of the values on the right and left of the cell interface. The source is approximated by a piecewise constant value. After making these assumptions, the first order 1D finite volume discretization for the cell $j$ reads

[TABLE]

with

[TABLE]

The numerical fluxes at the interfaces can be calculated with different Riemann solvers, e.g., Roe as in [38], Lax-Friedrich as in [50] or HLL as in the results of this paper using the standard discretization. We note here that the election of the Riemann solver for this problem has a small impact in the results as the numerical dissipation is dominated by the low-Mach regime of electrons, as it will be shown in the results. TVD reconstruction can improve the results of the standard discretization as shown in [50, 38].

As explained in [50, 52], the stability of the time discretization is restricted by a CFL condition that takes into account the convective scales of both fluids and the characteristic time scales of the source terms. The convective CFL reads

[TABLE]

Where $\bar{c}_{\EuFrak{e},\EuFrak{i}}$ is the dimensionless speed of sound of electrons and ions, i.e., $\varepsilon^{-1/2}$ and $\kappa^{1/2}$, respectively. Note that due to the mass disparity between ions and electrons, the CFL condition of the electrons is typically more restrictive than this of the ions. Therefore, the convective CFL condition of electrons uses the maximum eigenvalues of electrons that are $|\lambda^{\EuFrak{e}}_{max}|=\max\left(|\bar{u}_{\EuFrak{e}}+\varepsilon^{-1/2}|,~{}|\bar{u}_{\EuFrak{e}}-\varepsilon^{-1/2}|\right)$. Similarly, the source terms impose a constraint in the time step. The stability condition for the electrostatic force is related to the resolution of the electron plasma wave [51]

[TABLE]

Finally, the ionization term has the stability constraint as follows

[TABLE]

The stability condition reads

[TABLE]

As it can be seen in Table 1, if the cell size is larger than the Debye length, the most restrictive constraint is the resolution of the electron plasma waves. If the Debye length is resolved, then the convective condition is sufficient to fulfil the condition $\Delta t<\omega^{-1}_{pe}$. As shown in [38, 70], when this scheme does not spatially resolve the Debye length, the simulation leads to large spurious charge separation errors that can excite plasma modes, leading to an erroneous solution. Similarly, the truncation error of the upwind discretization of the fluxes of the electrons leads to a large error in the flux of the electrons. As it will be shown in the results and previously noted in [48, 50], this is due to the low-Mach regime of the electron when the bulk speed is much smaller than the thermal speed.

4 Acoustic/transport operator splitting strategy

We present an alternative to the standard discretization of section 3. We present a novel operator splitting strategy that decouples the acoustic and transport phenomena of electrons. In 1D, this method is analogous to an explicit Lagrange-Projection [57] method. Nevertheless, the present splitting does not need a moving Lagrangian mesh and can be naturally expressed for multi-dimensional problems (see Chalons et al. [57] for the extension of the operator splitting method to multiple dimensions).

We propose to approximate the system of equations (2a)-(2e) in the successive solution of the following systems. The first step solves for the electron acoustic system together with the Poisson equation and the second solves for the electron transport (advection) and the ion equations. It should be noted that the first system contains the small scales related to the parameters $\varepsilon$ and $\lambda$, whereas the second system solves for the slow dynamics, of order $\mathcal{O}(1)$, as explained in Section 2.2.

Electron acoustic and electrostatic system

[TABLE]

Electron transport and ion dynamics

[TABLE]

where the electron pressure $\bar{p}_{\EuFrak{e}}=\varepsilon^{-1}\bar{n}_{\EuFrak{e}}$, in eq. (20b).

Strategy to solve the equations

Given a state at the time $(\bar{n}_{\EuFrak{e}},\bar{n}_{\EuFrak{e}}\bar{u}_{\EuFrak{e}},\bar{n}_{\EuFrak{i}},\bar{n}_{\EuFrak{i}}\bar{u}_{\EuFrak{i}},\bar{\phi})^{n}_{j}$ at the time $t^{n}$ and the cell center $x_{j}$. The scheme is split into

By numerically solving the system (20), we update the state $(\bar{n}_{\EuFrak{e}},\bar{n}_{\EuFrak{e}}\bar{u}_{\EuFrak{e}},\bar{\phi})^{n}_{j}$ to the value at $t^{\ast}$, i.e., $(\bar{n}_{\EuFrak{e}},\bar{n}_{\EuFrak{e}}\bar{u}_{\EuFrak{e}},\bar{\phi})^{\ast}_{j}$. 2. 2.

By numerically solving the system (21), we update the state $(\bar{n}_{\EuFrak{e}},\bar{n}_{\EuFrak{e}}\bar{u}_{\EuFrak{e}},\bar{n}_{\EuFrak{i}},\bar{n}_{\EuFrak{i}}\bar{u}_{\EuFrak{i}},\bar{\phi})^{\ast}_{j}$ to the value $t^{n+1}$, i.e., $(\bar{n}_{\EuFrak{e}},\bar{n}_{\EuFrak{e}}\bar{u}_{\EuFrak{e}},\bar{n}_{\EuFrak{i}},\bar{n}_{\EuFrak{i}}\bar{u}_{\EuFrak{i}},\bar{\phi})^{n+1}_{j}$.

In the following, we present the properties and description of the two steps of the numerical scheme.

4.1 Properties and discretization of the electron acoustic and electrostatic system

The system eq. (20) can be written in conservative form, i.e., the compression terms written as fluxes. To do this, we define the variable $\tau_{\EuFrak{e}}\equiv 1/n_{\EuFrak{e}}$. In this variable, the problem (20) reads

[TABLE]

We approximate $\tau_{\EuFrak{e}}(x,t)\partial_{\bar{x}}$ by the solution at time $t^{n}$, i.e., $\tau_{\EuFrak{e}}(x,t^{n})\partial_{\bar{x}}$. We define $dm=n_{\EuFrak{e}}(x,t^{n})dx$. By using this new variable, we obtain the following system of equations for the electrons

[TABLE]

The left-hand-side of the system is conservative in this new variable. The eigenvalues of the homogeneous system are $\lambda^{acous}=(-\varepsilon^{-1/2},\varepsilon^{-1/2})$. For that reason, the system is called acoustic. Note that the source term in equation (23b) is still with the previous space variable $x$. As this term cannot be written in conservative form in the variable $m$, it will be treated in the discretization as a source term.

4.1.1 Discretization of the electron acoustic and electrostatic system

The system (23) with Poisson’s eq. (22c) will be discretized in time by using a semi-implicit discretization. In this approach, we want to discretize the Lorentz force in the electrons and the electron density in Poisson’s equation implicitly, as follows

[TABLE]

The full discretization of system (24) at $t^{\ast}$ is given by numerical scheme 1.

Numerical Scheme 1**.**

The discretization of the electron acoustic and electrostatic system (20) reads as follows. The discretized electron density is

[TABLE]

with the discretization of the pressure Laplacian as

[TABLE]

and the velocity at the cell interface

[TABLE]

where $f_{a}(M_{\EuFrak{e}})=\mathcal{O}(M_{\EuFrak{e}})$. For this function, we choose the formula proposed by Liou [71], as follows

[TABLE]

where

[TABLE]

The cut-off Mach $M_{co}$ is chosen to avoid having a null numerical dissipation when the electrons have zero velocity.

The discretized potential is solved by the following expression

[TABLE]

The discretized electron velocity is computed as follows

[TABLE]

where the gradient of the electric potencial is discretized as

[TABLE]

and the pressure flux from the expression

[TABLE]

The previous numerical scheme is obtained as follows. In order to build a scheme with an asymptotic-preserving property, we follow a similar approach to this of Dimarco et al. [58]. The divergence of the velocity at time $t^{\ast}$ in eq. (24a) is computed by taking the divergence of eq. (24b). This divergence of the velocity reads

[TABLE]

We discretize the term due to the Lorentz force at the time $t^{\ast}$ by using Poisson’s eq. (24c), as follows

[TABLE]

Note that the ion density is computed at time $t^{n}$ in the Poisson equation. This approximation is justified by the fact that the ions move slower than the electrons in the scales for the acoustic step. The divergence of the velocity is then introduced in the equation for the conservation of mass in the acoustic step. This reads

[TABLE]

We note that the LHS of the equation is the same as in an explicit discretization, whereas the RHS are terms that appear due to the implicitation of the scheme. The first term of the RHS is a diffusion term that is due to the Lorentz force and the second is a diffusion term due to the pressure gradient. We highlight the fact that the diffusion due to the Lorentz force, which will stabilize the scheme, appears without the need of solving Poisson’s equation and additionally, it is linear in $\tau_{\EuFrak{e}}^{\ast}$.

In order to discretize in space eq. (36), we follow the formalism of Chalons et al. [57].

[TABLE]

The pressure diffusion term is discretized as follows. By changing from the variable $m$ to the space variable, the pressure gradient reads $\partial_{m}\bar{p}_{\EuFrak{e}}=\bar{n}_{\EuFrak{e}}^{-1}\partial_{\bar{x}}\bar{p}_{\EuFrak{e}}$. By assuming the electrons to be isothermal, the pressure is $\bar{p}_{\EuFrak{e}}=\bar{n}_{\EuFrak{e}}\varepsilon^{-1}$ and the pressure gradient in the mass variable $\partial_{m}\bar{p}_{\EuFrak{e}}=\varepsilon^{-1}\partial_{\bar{x}}\ln\bar{n}_{\EuFrak{e}}$. Therefore, the Laplacian in the mass variable reads $\partial^{2}_{mm}\bar{p}_{\EuFrak{e}}=(\varepsilon\bar{n}_{\EuFrak{e}})^{-1}\partial^{2}_{\bar{x}\bar{x}}\ln\bar{n}_{\EuFrak{e}}$. The space discretization of this term is based on centered difference as follows

[TABLE]

Similarly, the velocity at the interface is approximated by a Riemann solver that is described in Appendix A. The original velocity at the interface reads

[TABLE]

where the subscripts $R$ and $L$ refer to the values on the right and left of the interface $j+1/2$. As explained in Appendix B, the truncation error is proportional to $\varepsilon^{-1/2}\partial^{2}_{\bar{x}\bar{x}}\tau_{\EuFrak{e}}\Delta x$. As suggested in [57], in order to control this error when $\varepsilon$ is small, we rescale the numerical error produced by an imbalanced numerical dissipation in the low Mach regime. This yields

[TABLE]

where $f_{a}(M_{\EuFrak{e}})=\mathcal{O}(M_{\EuFrak{e}})$. For this function, we choose the formula proposed by Liou [71], as follows

[TABLE]

where

[TABLE]

The cut-off Mach $M_{co}$ is chosen to avoid having a null numerical dissipation when the electrons have zero velocity.

With this equation, the electron density at the time $t^{\ast}$ is calculated as follows

[TABLE]

After solving this equation, the electric potential is computed by solving Poisson’s by the linear system described in eq. (30).

The last part of the acoustic step is to solve the electron velocity equation. For this, we use an approximate Riemann solver in order to discretize the pressure flux in eq. (24b), as follows

[TABLE]

The gradient of the electric potencial is discretized as

[TABLE]

The pressure flux from the Riemann solver in Appendix A reads

[TABLE]

As done previously, the truncation error of eq. (44) is proportional to $\left(\varepsilon^{-1/2}\partial^{2}_{\bar{x}\bar{x}}\bar{u}_{\EuFrak{e}}\Delta x\right)$ (see Appendix B). In this case, we choose to balance the numerical dissipation with the technique proposed by Liou [71].

[TABLE]

As in Chalons et al. [57], by preconditioning the numerical dissipation with the factor $f_{a}$, the truncation error of the numerical system is of order $\mathcal{O}(\Delta x)$.

As the Lorentz force is treated implicitly, the scheme is shown to be unconditionally stable for the plasma wave frequency. Additionally, with the preconditioning of the low-Mach regime, the electrons can have a larger time step, not limited by the sound waves, as it will be shown in the results. The linear stability analysis to determine the CFL condition is not trivial due to the non-linear term $\left[\partial^{2}_{mm}\bar{p}_{\EuFrak{e}}\right]_{j}^{n}$ and hence is left for a future work.

4.1.2 Asymptotic-preserving property of the acoustic and electrostatic system

In this section, we prove the consistency of the acoustic and electrostatic system when $\chi$ and $\varepsilon$ tend to zero with respect to the asymptotic limit of system (4) explained in Section 2.2.

Property 1**.**

In the limits $\varepsilon\rightarrow 0$ and $\chi\rightarrow 0$, the discretized electron density at $t^{\ast}$ reads

[TABLE]

Proof.

First, we may rewrite $f_{a}(M_{\EuFrak{e}})$ at the interface $j+1/2$ in the limit $\varepsilon\rightarrow 0$ as

[TABLE]

Then, injecting (40) and (47) into (43) leads to

[TABLE]

One observes that $\bar{n}_{\EuFrak{e}_{j}}^{\ast}$ is a function of $\varepsilon$ and $\chi$ continuous in $(0,0)$, then the two limits commute and (48) holds. ∎

We highlight that eq. (48) is the quasi-neutral limit that occurs at scales $\mathcal{O}(\chi^{-1})$, as shown in eq. (4a).

Property 2**.**

Assuming $\bar{n}_{\EuFrak{e}_{j}}^{n}=\bar{n}_{\EuFrak{i}_{j}}^{n}\neq 0$, then, in the limits $\varepsilon\rightarrow 0$ and $\chi\rightarrow 0$, the discretized electric potential at $t^{\ast}$ satisfies

[TABLE]

where the operator $[\partial^{2}_{\bar{x}\bar{x}}]_{j}$ stands for the central discretization in space at the cell $j$. This expression is the discretization of the divergence of eq. (4e).

Proof.

By injecting eq. (43) into eq. (30), we obtain the following expression for Poisson’s equation discretization

[TABLE]

where the coefficients $K_{j}^{i}$ are defined in (49).

This is again continuous in $(\varepsilon,\chi)=(0,0)$, thus the two limits commute and provide

[TABLE]

Using property 1, we may replace in this expression $\bar{n}_{\EuFrak{i}_{j}}^{n}$ by $\bar{n}_{\EuFrak{e}_{j}}^{\ast}$ in this expression.

And using the continuity of the logarithm, the hypothesis $\bar{n}_{\EuFrak{e}_{j}}^{n}=\bar{n}_{\EuFrak{i}_{j}}^{n}\neq 0$ and property 1, we may replace $\bar{n}_{\EuFrak{e}_{j}}^{n}$ in this expression to obtain (50).

∎

As seen above, the advantage of this method is that it is AP by solving a standard Poisson equation, i.e., eq. (30). However, when this equation is discretized by a standard explicit discretization, i.e., eq. (15), the electric potential tends to infinity when $\chi\rightarrow 0$ if the $\Delta x$ is not small enough, resulting in a large numerical error.

In the discrete electron momentum equation of the acoustic step in eq. (44), we can find problems of consistency due to the small Mach number ($\varepsilon\to 0$). In the following property we explain how these problems are solved by preconditioning the numerical dissipation.

Property 3**.**

If $f_{a}=\mathcal{O}(\varepsilon^{1/2})$, then, the truncation error of the electron momentum equation in the acoustic step in eq. (22b) is uniform with respect to $\varepsilon<1$.

Proof.

The proof is analogous to the one presented in [57, 72]. We assume that the variables describe a smooth flow so we can expand them in Taylor series to study the truncation error (see Appendix B for the detailed derivation of the truncation error). With this development, we use the numerical pressure flux of (47), and we neglect second order terms, i.e., $\mathcal{O}\left(\Delta x^{2}\right)$, $\mathcal{O}\left(\Delta t^{2}\right)$. Consequently, the discretization is consistent with

[TABLE]

The function $f_{a}$ in eq. (41) is a continuous function of $\varepsilon$, as follows,

[TABLE]

Assuming that $C_{j+1/2}^{1}$ and $C_{j+1/2}^{2}$ are of order one, then, $f^{a}=\mathcal{O}(\varepsilon^{1/2})$ and therefore the truncation error is uniform with respect to $\varepsilon<1$. ∎

4.2 Discretization of the electron transport and ion dynamics system

In this subsection we present the discretization for the system (21). As explained before, this system represents the larger scales related to the transport of electrons and the ion scales.

4.2.1 Discretization of the electron transport system

The transport system for the electrons (eqs.(21a)-(21b)) is a hyperbolic system. By using the upwind scheme proposed in Chalons et al. [57], we can discretize this system by the following numerical scheme.

Numerical Scheme 2**.**

The discretized transport system for the electrons (eqs.(21a)-(21b)) reads

[TABLE]

where the velocity in the interface is defined in eq. (40) and

[TABLE]

As in Chalons et al. [57], the CFL condition associated to the transport system is

[TABLE]

where he superscript stands for $b^{\pm}=\frac{b\pm|b|}{2}$.

4.2.2 Discretization of the ion dynamics system

The low-temperature plasmas have an additional difficulty in the integration of the ion dynamics. As the plasma potential scales with the electron temperature, the Lorentz force on the ions is in general much larger than the ions pressure gradient. Naive discretizations of the source term may create numerical instabilities (see, e.g., Section 5.2 for the present problem or [60, 61, 62, 73, 74] in a general framework). In order to tackle this problem, we propose an upwind discretization [60] of the source term in eqs. (21c)-(21d),as presented in the following numerical scheme.

Numerical Scheme 3**.**

The discretization of the ion dynamics system (21c)-(21d) reads

[TABLE]

where

[TABLE]

The discretization of the convective fluxes ${\mathbf{F}_{\EuFrak{i}}}$ is performed with an approximate Riemann HLL scheme [75]. The discretization of the collisional source term is at the cell center, i.e.,

[TABLE]

Alternatively, as the Lorentz force is dominant in low pressure discharges, we propose a well-balanced scheme as in [60]. We decompose the source term into the contribution at the interfaces. For the sake of simplicity of notation, we do not include the subscript corresponding the time discretization. The discretization reads

[TABLE]

where the functions $\mathcal{G}^{L/R}$ read

[TABLE]

The matrix $\mathbf{A}$ is the Jacobian matrix of the convective flux and $\hat{\mathcal{G}}$ is the source term computed in the interface. They are defined as

[TABLE]

where

[TABLE]

In the particular case of isothermal flow with no mass source term, $\mathcal{G}^{L/R}$ simplify to

[TABLE]

Note that the velocity of the ions is evaluated at the interface between the cells. As $\text{sign}(\delta u_{\EuFrak{i}})$ is a discontinuous function, in numerical experiments, we observe that the $\text{sign}(\delta u_{\EuFrak{i}})$ function introduces discontinuities in the solution that are advected by the fluids, as the model has no physical dissipation. For that reason, we use the following approximation for the $\text{sign}(\delta u_{\EuFrak{i}})$ function

[TABLE]

where $\delta u_{\EuFrak{i}}=\bar{c}_{\EuFrak{i}}+u_{\EuFrak{i}}$ or $\bar{c}_{\EuFrak{i}}-u_{\EuFrak{i}}$ and $u_{\EuFrak{i}_{\infty}}$ is a characteristic speed of the problem.

4.3 Summary of the numerical scheme

In Fig. 5, we summarize the steps followed in order to obtain the numerical results corresponding to the AP scheme. We highlight that the steps in the electron acoustic system need to be done sequentially. The order of steps 4) and 5) can be inverted obtaining the same result.

5 Numerical results

In the present section, we present the numerical results of the proposed scheme. We first simulate a periodic two-stream perturbation in a thermal plasma ($T_{\EuFrak{e}}=T_{\EuFrak{i}}$). In this problem we study: (1) the convergence of our numerical scheme with the time step and the mesh size, (2) the asymptotic behaviour and consistency of the scheme when $\chi\rightarrow 0$, and (3) to assess the low-Mach correction (limit when $\varepsilon\rightarrow 0$) in the electron fluxes.

Second, we simulate the same instability in a low-temperature plasma order to study: (1) the performance of the well-balanced scheme in the ions, (2) assess the spatial convergence in a low-temperature plasma, and (3) to evaluate the computational performance as compared to a standard discretization.

Finally, we propose a numerical set-up to simulate a plasma sheath between two floating conducting walls. In this case, we simulate a physical domain that contains charge separation regions. We show that our numerical set-up captures the physics of the classical plasma sheath [37], computing the potential drop as predicted by the theory. As the numerical domain contains both a quasi-neutral and charged regions, in this test we assess the ability of the scheme to represent both limits.

5.1 Propagation of a two-stream periodic perturbation in a collisionless thermal plasma

The two-stream instability occurs in a uniform plasma of density $n_{0}$ where the electrons have a relative velocity $u_{0}$ with respect to the ions [12]. Perturbations of the equilibrium state can lead to instabilities for sufficiently long wavelengths $k$. In our case, we will restrict our work to study the propagation of a periodic perturbation in a stable plasma. This case has been previously used for the convergence study of the two-fluid model in a number of works [49, 48, 50]. Our set of equations is slightly different to these previous works as we study the isothermal case whereas in Crispel et al. [49, 48, 50] the fluid model assumes an isentropic law. As in [1], in our case, the dispersion relation is found to be

[TABLE]

where the coefficients $\alpha_{i}(k;\varepsilon,\kappa,n_{0},v_{0},\chi)$ read

[TABLE]

Firstly, we consider a perturbation in a thermal plasma, i.e., where the ions have the same temperature as the electrons. We consider a mass ratio of $\varepsilon=10^{-4}$ and a domain of length $L=1$ with $10^{4}$ Debye lengths, i.e., $\lambda=10^{-4}$. The background plasma has a normalized density of $n_{0}=1$ and the electrons travel at a velocity of $u_{0}=1$ that corresponds to Bohm’s speed at these conditions. We fix the wavelength of the perturbation to $k=2\pi$. Under these conditions, the characteristic polynomial, eq. (72), has four real solutions, i.e., the perturbation is stable. There are two high frequency solutions, which correspond to the electron plasma waves, and two low frequency that are quasi-neutral. For our study, we use one of the low frequency solutions, i.e., $\omega=8.8857268$.

The initial field is the analytical solution for a small perturbation at the given frequency, wavelength, and plasma conditions. In the present case, this solution is

[TABLE]

In Fig. 6, we present the solution of the Lagrange-projection AP scheme without the low-Mach correction for the electrons at $t=0.7071$ which corresponds to one period of the wave. We present the results with four different mesh resolutions $N=50$, $100$, $200$, and $400$ points. We highlight that the mesh size is not adapted to the Debye length as there are $200$ Debye lengths inside one cell in the case of lowest resolution. The simulations use a CFL${}^{conv}=0.7$, which corresponds to a time step of $\Delta t=1.4141\times 10^{-4}$, $0.7071\times 10^{-4}$, $0.3535\times 10^{-4}$, and $0.176775\times 10^{-4}$. We also notice that this time step is much larger than the period of a plasma wave, i.e., $\omega_{p\EuFrak{e}}^{-1}=10^{-6}$. A standard simulation that resolved the Debye length with 10 mesh points and the plasma frequency with $10$ time-steps, would have $2000$ more mesh points and a time step $10^{3}$ times smaller.

In Fig. 7, we present the solution for $\bar{n}_{\EuFrak{e}}$, $\bar{n}_{\EuFrak{i}}$, $\bar{u}_{\EuFrak{e}}$, and $\phi$ at $t=0.0626$, i.e., after $900$ iterations of the standard HLL scheme described in Sec. 3 that uses $100$ mesh points and a CFL${}^{conv}=0.7$. These results illustrate that indeed, using the same conditions than in the Lagrange-projection AP scheme with $N=100$, the standard discretization is unstable. It is also important to note that the scheme becomes unstable very fastly due to the rapid unstable response of the electrons.

The same simulation is presented with using the Lagrange-projection AP scheme with the low-Mach correction for the electrons in Fig. 8. We show the previous four discretization resolutions with CFL${}^{conv}=0.7$. With the low-Mach correction, the solution for the electron velocity converges to the analytical solution when the mesh size is decreased, showing a much smaller error.

As explained in [54], when the numerical diffusion is rescaled for the low-Mach regime, the CFL condition is not restricted by the acoustic waves. This allows for the adaptation of the CFL to the eigenvalue $\lambda^{\EuFrak{e}}=\bar{u}_{\EuFrak{e}}$ rather than to $\lambda^{\EuFrak{e}}=\bar{u}_{\EuFrak{e}}\pm\bar{c}_{\EuFrak{e}}$. In Fig. 9, we present the comparison of the solution for CFL${}^{conv}=0.7,~{}7.1,$ and $28.5$ with the low-Mach correction and $100$ cells. We highlight that the largest CFL corresponds to an ion CFL${}_{\EuFrak{i}}^{conv}=0.28$. The results show that the numerical diffusion is lower when the CFL is adapted to the ion scales. In the simulation with the largest CFL, the time step is $\Delta t=1.99\times 10^{4}\omega^{-1}_{p\EuFrak{e}}$, i.e., $10^{4}$ larger than an explicit non-AP discretization.

In Fig. 10, the convergence of the $L_{2}$ error norm is presented for three different cases: the AP scheme without low-Mach correction for the electrons and CFL${}^{conv}=0.7$ (left), the AP scheme with low-Mach correction for the electrons and CFL${}^{conv}=0.7$ (center), and the AP scheme with low-Mach correction for the electrons and CFL${}^{conv}=28.5$ (right). We clearly see that the low-Mach correction of the electrons reduces the error of the electron velocity by one order of magnitude. When the CFL${}^{conv}=0.7$, the convergence of the error corresponds to a first order scheme. In the case of CFL${}^{conv}=28.5$, the error norm is reduced, especially in the ion quantities. However, we observe that the error convergence is slower than a first-order scheme. This is caused by the low-Mach regime of the ions. Since the perturbation of the ion velocity is small, the Mach of the ions varies from $M_{\EuFrak{i}}=0$ to $0.03$. This could be fixed with a low-Mach correction, similar to the one performed in the electrons. Nevertheless, the interest of this paper is low-temperature plasmas, which means that in general, the velocity of the ions is much larger than the ion thermal speed. For that reason, the low-Mach regime of the ions is not treated in the present paper. In the following example, in a low-temperature plasma, we will show how this convergence improves when the Mach of ions is larger.

5.2 Propagation of a two-stream periodic perturbation in a collisionless low-temperature plasma

As mentioned above, in low-temperature plasmas, the bulk velocity of ions is usually much larger than the ion thermal speed. When the ion velocities are comparable or larger than the ion thermal speed, this can produce numerical problems due to the imbalance between the cell-centered source terms and the upwinded fluxes. In order to illustrate the performance of the numerical scheme in a case of low-temperature plasmas, we simulate a two-stream periodic perturbation in a plasma with an ion-to-electron temperature ratio of $\kappa=10^{-4}$.

We use the dispersion relation of eq. (72) with $\varepsilon=\kappa=\lambda=10^{-4}$. We choose $n_{0}=1$ and, in this case, $v_{0}=0.7$. Note that $v_{0}=1$ produces a resonance and therefore the perturbations would be too large to study a linear wave propagation. With these conditions, the initial condition reads

[TABLE]

The period of the oscillations is $T=0.99995$. We run the simulation and compare the results to the initial field after one period.

In Fig. 11, we show the comparison between the scheme with the well-balanced source term in eq. (3) and a cell-centered implementation of the source. Both simulations use CFL${}^{conv}=25$ and $N=100$. The approximation of the sign function used follows eq. (71) with $u_{\EuFrak{i}_{\infty}}=v_{0}$. The results show that the solution without the well-balanced implementation maintain the quasi-neutrality of the solution. However, it develops spurious numerical oscillations. These oscillations are not present with the well-balanced implementation.

Under these conditions, the Mach regime of the maximum velocity of the ions is $M_{\EuFrak{i}}=0.33$ and the electrons $M_{\EuFrak{e}}=0.007$. The convergence for different mesh sizes is presented in Fig. 12. In this case, the ions are in a larger Mach regime and therefore the compressible solver is less dissipative than in the thermal plasma case ($\kappa=1$). As a result, as shown in Fig. 13, the error convergence is indeed first-order accuracy in space for this scheme, despite the extremely low Mach regime of the electrons and the limit of $\lambda\rightarrow 0$.

Finally, we compare the computational performance of the present numerical scheme. We highlight that the Poisson solver with periodic boundaries is solved by adding a Lagrange multiplier to the system that imposes the total charge in the domain to be zero and by solving the linear system with a LU decomposition implemented in the library LAPACK [76]. We choose to use a linear solver instead of a spectral method as we consider it to be more representative of a realistic multidimentional case with non-periodic boundary conditions. This is due to the fact that most of the iterative methods for solving linear systems do not scale linearly with the number of points and therefore the computational performance will be largely penalized by number of points.

In Table 2, we present the computational performance for three different cases: (1) a first-order standard discretization that resolves the Debye length with one cell and the plasma frequency with ten time steps, (2) the same first order standard that a mesh cell that is $\Delta x=100$ and $\Delta t=0.1\omega^{-1}_{p\EuFrak{e}}$ (the scheme needs to resolve $\omega^{-1}_{p\EuFrak{e}}$ for stability reasons), and (3) the AP scheme with $\Delta x=100$ and $\Delta t=2500\omega^{-1}_{p\EuFrak{e}}$. We highlight that the result of the standard case without resolving the Debye length diverges as the numerical error in the electron velocity becomes very large and the perturbation is not linear after a time-step.

In Table 2, we can see that the AP scheme produces a dramatic improvement as compared to a standard case when the plasma is quasi-neutral and the characteristic size of the phenomena is very large as compared to the Debye length. As the AP scheme does not need to resolve the Debye length to preserve the quasi-neutrality limit, the AP discretization needs $100$ times less mesh points and $25000$ times less time steps than a scheme that resolves the scales related to the Debye length. Consequently, the CPU time taken for simulating a period of the wave is of the order of $10^{9}$ times faster than a first-order standard scheme with similar accuracy.

5.3 Collisionless isothermal sheath

The previous numerical experiments simulate a quasi-neutral plasma. Nevertheless, as discussed in Section 2.3, the set of equations allows for capturing charge separation effects, as shown in Fig. 1. For length scales much greater than the Debye length, the plasma behaves in general as quasi-neutral. However, when a surface is in contact with the plasma, a charged boundary layer called the plasma sheath is formed.

We consider a plasma in a 1D domain of unitary size with an electron-to-ion mass ratio $\varepsilon=10^{-5}$, ion-to-electron temperature ratio $\kappa=0.025$, and $\chi=10^{-4}$. The electron flux is imposed on both boundaries by the number of particles crossing the plane with positive velocity component

[TABLE]

The electron and ion density and the ion flux have a Neumann boundary on both sides. Alternatively, the electric potential has a Dirichlet condition on both sides $\bar{\phi}(\bar{x}=0,t)=\bar{\phi}(\bar{x}=1,t)=0$. In the results, the potential is referenced to the potential in the plasma at $\bar{x}=1/2$.

In the sheath theory [37], the ionization is an eigenvalue of the problem. This means that there is only one ionization rate that can produce a steady solution. In our simulation set-up, we find this eigenvalue by changing the ionization frequency such the balance between ionization and ion flux occurs at every time step. Consequently, the ionization frequency at $t^{n}$ is calculated as

[TABLE]

The flow field is initialized as follows

[TABLE]

As the assumptions of the classical theory are not taken by our numerical model, we might expect a result that is slightly different from eq. (13). For this reason, we will use as a reference solution a high-order highly resolved simulation described in Section 2.3. In Fig. 14, we show the steady solution reference solution that uses a standard discretization, as in Section 3, with TVD third-order reconstruction and third-order TVD Runge-Kutta [77] scheme. In the reference solution, the domain is resolved with $N=10^{4}$ mesh points and CFL=$0.9$. By doing this, we ensure that the Debye length is resolved with $100$ points and so is the electron plasma period. Since the scheme is high order and highly resolved, the discretization errors are expected to be much smaller than a first order scheme with lower resolution.

The solution shows two plasma sheaths beside the left and the right boundaries of around $6$ – $7$ Debye lengths width. As shown in the densities, the plasma is quasi-neutral in the rest of the domain, whereas it is positively charged in the sheath. The plasma potential at the wall is very similar to the theory, despite the electron inertial and the ion pressure gradient effects that are neglected in the theory. The flux of ions and electrons are equal in the whole domain. The velocities are also equal in the quasi-neutral region and they differ when the ions reach $|\bar{u}_{\EuFrak{i}}|=1$, i.e., the Bohm’s velocity in our units. This point is the plasma-sheath transition that agrees with the Bohm’s criterion [36].

In Fig. 15, we present the converged results using a first-order standard discretization with CFL $=0.1$ and $N=10^{3}$. The local error is calculated with the reference solution. The main difference with the reference solution can be clearly observed in the electron flux. The error of the electron flux is as large as the actual solution. This is due to the poor behaviour of the HLL scheme under low-Mach conditions, as explained in the previous examples.

The results using the AP Lagrange-projection scheme with CFL ${}^{conv}=0.4$ and $N=10^{3}$ are presented in Fig. 16. We note that, as we want to resolve the plasma sheath, in this case we resolve the Debye length in our discretization. The results show that the plasma sheaths are properly captured with the scheme and the simulation agrees with the reference solution. Additionally, the electron flux is correctly captured, reducing by two orders of magnitude the numerical error. The rest of magnitudes have similar values to the reference solution, with a first-order discretization and ten times lower resolution than the reference solution.

Even though the error in the electron flux seems to not have a large impact in the other variables in Fig. 15. This is only true for the steady state solution. However, if we analyze the evolution of the simulation to steady state, the solutions are radically different even though they have similar steady state solution. In Fig. 17, the solution at $t=75$ is presented for the reference solution, the HLL first order scheme and the AP scheme. We observe that whereas the solution of the AP scheme is very similar to the reference solution, the HLL standard scheme develops a plasma instability. This plasma instability is identified as a two-stream instability that is induced by the numerical error in the electron velocity. The instability is not present in the HLL scheme when the CFL is reduced to CFL$=0.01$ (see [38].

This numerical example illustrates that the AP scheme provides a more stable discretization that avoids spurious numerical instabilities in the quasi-neutral region, while still being able to accurately capture the regions with charge separation. Additionally, the scheme increases dramatically the accuracy of the computation of electron dynamics as well as increases by several orders of magnitude the computational efficiency.

6 Summary and conclusions

In low-temperature plasma applications, the plasma-wall interaction is fundamental to explain the conditions inside the bulk of the plasma. For this reason, the plasma sheaths need to be included in the computational domain to fully explain the characteristics of the discharge. Consequently, the numerical model needs to be able to resolve the Debye length in order to capture the potential drop that forms beside the walls of the device. Similarly, the electron inertia is seen to play a fundamental role in low-pressure conditions. However, in most parts of the domain, the plasma behaves as quasi-neutral with massless electrons. A numerical method that is able to tackle accurately and efficiently these different regimes is a long-standing problem in computational plasmas. In this paper, we have presented a numerical scheme for the isothermal plasma equations coupled to Poisson’s equation that proves to be accurate and stable both in the quasi-neutral and charge regimes, and both when electrons behave as massless and when electron inertia plays a role.

In Section 2.2, we have discussed the asymptotic behaviour of the plasma equations when both the electron mass and the Debye length tend to zero. We find a different asymptotic behaviour as compared with previous work and that the zero-th order solution of the electric potential depends only on the electron dynamics when the electron-to-ion mass ratio is included in the asymptotic expansion. This relation is known in plasma physics as Boltzmann density distribution. Alternatively, the electron momentum equation at $\mathcal{O}(1)$ is similar the low-Mach limit of the Euler equations.

Based on this asymptotic behaviour, we have proposed a novel operator splitting strategy that relies on the Lagrange-projection scheme. Our approach solves for the electron acoustic system coupled to the Poisson equation in a first step and the electron transport together with the ion dynamics in a second step. We show that the first step preserves the asymptotic limit when both the electron mass and the Debye length tends to zero. This means that even when these scales are not resolved, the numerical scheme is able to find the correct asymptotic solution, without the need for an implicit scheme. Consequently, the scheme is able to achieve a great improvement in the accuracy of the solution as the numerical dissipation is properly scaled to the physics. More importantly, the numerical method improves dramatically the computational efficiency when these scales are not needed to be resolved.

The numerical scheme has the advantage of solving a standard Poisson equation, as compared to other AP and semi-implicit methods that solves for more complex parabolic or elliptic equations. Additionally, the acoustic and electrostatic step involves only the electron dynamics and, therefore, the proposed model can be used in hybrid approaches where the electrons are treated as a fluid and the ions with a kinetic approach. Moreover, we have proposed a well-balanced treatment of the Lorentz force in the ions that proves to be fundamental in order to avoid spurious numerical oscillations when the temperature of the ions is much lower than this of electrons.

In the numerical results, we have simulated a quasi-neutral two-stream perturbation in thermal and low-temperature plasmas. Additionally, we have simulated a plasma sheath in a low-temperature plasma. The two-stream perturbation allowed us for assessing the performance of the model when the plasma is quasi-neutral but the electron velocity differs from this of ions. We also have proved that the scheme is able to reproduce the numerical solution even when the Debye length and the electron sound speed are not resolved. This allows us for achieving an efficient numerical scheme as compared to a scheme that needs to resolve the scales related to the Debye length or the electron acoustic waves. In our numerical experiments, the AP scheme is $10^{9}$ faster to simulate a plasma period of a quasineutral wave than the standard discretization that needs to resolve the Debye length and the electron plasma frequency. Additionally, the electron velocity is properly captured, avoiding the problems associated to the low-Mach limit of electrons.

Finally, we have proposed a numerical set-up to simulate the steady solution of a plasma between two floating plates. The set-up is able to find self-consistently a solution that agrees with the analytical plasma sheath theory. This solution contains three different regions: (i) the quasi-neutral plasma where the ions and electrons have the same density and velocity, (ii) the plasma-sheath transition where ions are accelerated to the Bohm’s speed, and (iii) the sheath where the plasma is electrically charged and the electrons are accelerated to speeds comparable to their thermal speed. The proposed AP scheme is able to reproduce accurately these three regions. The standard discretization is not able to reproduce the electron flux as the standard numerical scheme is not able to reproduce the limit when the electron inertia is very small, even when the Debye length is properly resolved. Additionally, the standard scheme develops non-linear instabilities that are triggered by the numerical error in the electron velocity.

In the present paper, we have discussed the fundaments of an asymptotic preserving scheme for a fluid plasma model in low-temperature plasma conditions. We have highlighted that the ideas presented in this paper are the basis for future developments, extending the scheme to multi-dimensions, the inclusion of the energy equation, and the development of high-order schemes including this splitting strategy. Similarly, the scheme presented in the present paper would be specially advantageous for a non-uniform mesh where the quasi-neutral region of the plasma does not need to be resolved, whereas the boundaries resolve the Debye length in order to capture the plasma sheath.

Acknowledgements

The first author is funded by the postdoctoral fellowship from the Fondation Mathematique Jacques Hadamard (FMJH), LabEx Mathématique Hadamard - ANR11-LABX-0056 and LabEX LASIPS - ANR10-LABX-0032. The authors would like to thank Dr. Nagi N. Mansour and Dr. Samuel Kokh for very useful discussions during the 2018 NASA Summer Program at NASA Ames Research Center. T. M.  is supported by a “Chaire Jean d’Alembert” of University Paris Saclay at Ecole Polytechnique.

Appendix A: Approximate Riemann solver of the acoustic step

We can write the system of eq. (23) as

[TABLE]

Note that the pressure is written as function of the conservative variable $\tau_{\EuFrak{e}}$, i.e., $p_{\EuFrak{e}}=\varepsilon^{-1}\bar{n}_{\EuFrak{e}}=\varepsilon^{-1}\tau_{\EuFrak{e}}^{-1}$. We define the conservative variables of the electron acoustic step, the flux term, and the source term as

[TABLE]

The Jacobian matrix of the flux term in the LHS of eq. (79) reads

[TABLE]

The eigenvalues of this matrix are $\lambda=\left(-\bar{n}_{\EuFrak{e}}\varepsilon^{-1/2},~{}\bar{n}_{\EuFrak{e}}\varepsilon^{-1/2}\right)$. Consequently, the Lax-Friedrich numerical flux for the previous Riemann problem reads

[TABLE]

where the density at the interface is computed as $\bar{n}_{\EuFrak{e}_{j+1/2}}=\left(\bar{n}_{\EuFrak{e}_{R}}+\bar{n}_{\EuFrak{e}_{L}}\right)/2$.

Appendix B: Truncation error of the acoustic system

We study the truncation error of the upwind scheme for the acoustic step. We consider that the variables describe a smooth flow so that we can expand them in Taylor series the conservative variables. By using the Lax-Friedrich scheme of eq. (82), we obtain the following approximations in the mass equation

$\frac{\tau_{\EuFrak{e}}^{\ast}-\tau_{\EuFrak{e}}^{n}}{\Delta t}\approx\partial_{\bar{t}}\tau_{\EuFrak{e}}(t^{n},x_{j})+\partial^{2}_{\bar{t}\bar{t}}\tau_{\EuFrak{e}}(t^{n},x_{j})\Delta t+\mathcal{O}(\Delta t^{2})$

• *

$[\partial_{m}\bar{u}_{\EuFrak{e}}]^{n}_{j}\approx\frac{1}{\bar{n}_{\EuFrak{e}}}\partial_{\bar{x}}\bar{u}_{\EuFrak{e}}(t^{n},x_{j})+\varepsilon^{-1/2}\partial^{2}_{\bar{x}\bar{x}}\tau_{\EuFrak{e}}(t^{n},x_{j})\Delta x+\mathcal{O}(\Delta x^{2})+\mathcal{O}(\varepsilon^{-1/2}\Delta x^{3})$

• *

$[\partial^{2}_{mm}\bar{p}_{\EuFrak{e}}]^{n}_{j}\approx\frac{1}{\varepsilon\bar{n}_{\EuFrak{e}}^{2}}\partial^{2}_{\bar{x}\bar{x}}\bar{n}_{\EuFrak{e}}(t^{n},x_{j})+\mathcal{O}(\varepsilon^{-1}\Delta x^{2})$

Therefore, the truncation error of the discretization is

[TABLE]

Similarly, we obtain the following approximations for the equation for the velocity, expanded in Taylor series

$\frac{\bar{u}_{\EuFrak{e}}^{\ast}-\bar{u}_{\EuFrak{e}}^{n}}{\Delta t}\approx\partial_{\bar{t}}\bar{u}_{\EuFrak{e}}(t^{n},x_{j})+\partial^{2}_{\bar{t}\bar{t}}\bar{u}_{\EuFrak{e}}(t^{n},x_{j})\Delta t+\mathcal{O}(\Delta t^{2})$

• *

$[\partial_{m}\bar{p}_{\EuFrak{e}}]^{n}_{j}\approx\frac{\varepsilon^{-1}}{\bar{n}_{\EuFrak{e}}}\partial_{\bar{x}}\bar{n}_{\EuFrak{e}}(t^{n},x_{j})+\varepsilon^{-1/2}\partial^{2}_{\bar{x}\bar{x}}\bar{u}_{\EuFrak{e}}(t^{n},x_{j})\Delta x+\mathcal{O}(\Delta x^{2})+\mathcal{O}(\varepsilon^{-1/2}\Delta x^{3})$

• *

$[\partial_{x}\bar{\phi}]_{j}^{\ast}\approx\frac{1}{\varepsilon}\partial_{\bar{x}}\bar{\phi}(t^{n},x_{j})+\mathcal{O}(\varepsilon^{-1}\Delta x^{2})$

Consequently, the truncation error of the velocity equations reads

[TABLE]