# Comparison of high-order moment models for the ion dynamics in a bounded low-temperature plasma

Anatole Berger, N. Lequette, Thierry E. Magin, Anne Bourdon, A. Alvarez

Laguna

► To cite this version:

Anatole Berger, N. Lequette, Thierry E. Magin, Anne Bourdon, A. Alvarez Laguna. Comparison of high-order moment models for the ion dynamics in a bounded low-temperature plasma. Physics of Plasmas, 2025, 32 (10), 10.48550/arXiv.2505.10456. [https://doi.org/10.50689/2025-32-10.48550]

HAL Id: hal-05068896

https://hal.science/hal-05068896v2

Submitted on 31 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Comparison of high-order moment models for the ion dynamics in a bounded
low-temperature plasma

A. Berger,¹,² N. Lequette,¹ T. Magin,² A. Bourdon,¹ and A. Alvarez Laguna¹

1)Laboratoire de Physique des Plasmas, Centre National de la Recherche Scientifique,
 Sorbonne Université, École polytechnique, Institut Polytechnique de Paris,
 route de Saclay 91128 Palaiseau, France

2) Aeronautics and aerospace department, von Karman Institute for Fluid Dynamics,
Waterloosesteenweg 72, B-1640 Sint-Genesius-Rode, Belgium

(*Electronic mail: anatole.berger@lpp.polytechnique.fr)
Low-temperature plasmas often present non-equilibrium ion distribution functions due to
the collisions with the background gas and the presence of strong electric fields. This
non-equilibrium is beyond classical fluid models, often requiring computationally-intensive
kinetic simulations. In our work, we study high-order moment models in order to capture the
non-equilibrium state with a macroscopic set of equations, which is more computationally
efficient than kinetic simulations. We compare numerical simulations of different moment
closures: Grad's closure, the hyperbolic quadrature method of moments (HyQMOM), the
extended quadrature method of moments, and a method based on entropy maximization.
We assess the different closures for plasma applications and propose efficient numerical
discretizations. The numerical solution of the high-order moment models is compared to
kinetic simulations of an argon plasma between two floating walls at different pressure
regimes, from nearly collisionless to collisionally-dominated. In general, all the high-order
moment closures capture the ion transport with high fidelity as compared to the kinetic
simulations, providing an improvement as compared to classical fluid models. Classical
fluid closures such as the Fourier law for the heat flux is shown to be not suitable to capture
the sheath or the low pressure regime. In addition, the ability of each moment method
to reconstruct the velocity distribution function from the moments is assessed. The high-
order moment models are able to capture the non-equilibrium distributions in the bulk and
sheath with remarkable fidelity, dramatically improving classical fluid models while having
comparable computational cost. In particular, the HyQMOM shows to be a robust method
that provides an excellent comparison with the kinetic simulations of both the moments
and the distribution function in the bulk and the sheath.
# I. INTRODUCTION

Low-temperature plasmas are used in a wide variety of applications such as electric propulsion¹, plasma etching for microelectronics²,³, plasma-assisted combustion⁴, and biomedical applications⁵. In general, these plasmas present a low ionization level that results in conditions far from thermodynamical equilibrium, where the charged species collide more often with the neutral species of the gas than among themselves. In addition, low-temperature plasmas are often in contact with the walls of the reactor and other surfaces. The presence of the walls creates a plasma sheath with strong gradients and electric fields, driving the plasma species even further from thermodynamic equilibrium. For these reasons, the charged species present strongly non-Maxwellian distribution functions in the sheath regions, which plays an important role for some applications such as plasma etching or surface treatment, where the ion energy distribution has a large impact in the surface processes⁶.

Classical fluid models solve the particle number, momentum, and energy balance equations, i.e., 3 moments. The validity of classical fluid models for modeling these type of non-equilibrium plasmas is limited since, strictly speaking, they only consider Maxwellian velocity distribution functions (VDF) or small perturbations from the thermodynamic equilibrium, and are therefore unable to capture phenomena such as strong heat fluxes, pressure anisotropies or the non-Maxwellian VDF in the sheath. Despite these limitations, different theoretical models based on the fluid formalism have been proposed⁷,⁸ to describe bounded low-temperature plasmas (a summary can be found in Refs.²,³). More recent works⁹,¹⁰, have revisited these fluid models based on kinetic simulations. The kinetic simulations (such as particle-in-cell simulations¹¹,¹², semi-Lagrangian methods¹³, Boltzmann solvers¹⁴,¹⁵ or Fokker-Plank solvers¹⁶) are valid in the whole range of pressures but computationally much more expensive than fluid models. However, most of the fluid models used in the theory of bounded plasmas still consider a reduced number of moments (e.g., often assuming isothermal conditions or at most following empirical polytropic laws⁹,¹⁷). In addition, since they are based on a Maxwellian VDF, they are not able to provide a better approximation of the non-equilibrium distribution function of the ions at the wall.

A potential alternative to kinetic models is the method of moments. This method extends the fluid variables by considering higher-order moment equations, e.g., pressure tensor, heat-flux, etc. As a result the method of moments can capture more general VDFs than classical
fluid models, i.e. non-Maxwellian distributions, hence extending the validity of the model to non-equilibrium conditions. In addition, moment methods are, in general, computationally more efficient than kinetic models and without the statistical noise of particle-based kinetic methods. The difficulty associated to these models is in providing a consistent, mathematically stable, and numerically efficient model for the closing flux that appears in the last equation of the moment hierarchy as well as for the reconstruction of the VDF from the moments. Various closures for the moment equations have been proposed for different applications (see e.g. Ref.¹⁸ for a recent review). In particular, in partially ionized plasmas, Zhdanov¹⁹ used Grad’s moment method to derive general transport coefficients of neutral and charged species, in the linearized moment regime (i.e., close to thermodynamic equilibrium), also recently revisited by Hunana²⁰. Under conditions further from equilibrium we find some recent numerical works. Bocelli et al.²¹,²² have applied an interpolative maximum-entropy closure to electrons and ions in a magnetized low-temperature plasma. Taunay et al.²³ considered different quadrature methods of moments (QMOM) for a collisionless Vlasov-Poisson system. Alvarez Laguna et al.²⁴,²⁵ developed a regularized Grad’s model for electrons in a reacting partially-ionized plasma, later generalized in Ref.²⁶. Kuldinow et al.²⁷ have proposed a ten-moment model with a heat-flux closure. However, despite these efforts, the realistic modeling of the multicomponent collisional processes in the moment equations as well as the behavior of these models at different pressures in the presence of sheaths is still an open question for the moment closures. Similarly, the comparison with kinetic simulations of a bounded low-temperature plasma, representing the plasma-sheath transition, is still a challenge due to the presence of strong electric fields in the plasma sheath²⁸,²⁹. In addition, as mentioned above, the ability of moment models to capture the ion VDF inside the sheath has not been studied so far.

In this paper, we focus on the ion dynamics in a one-dimensional argon plasma between two floating walls in order to represent the plasma-sheath transition. This work follows a similar methodology to the work of Laplante and Groth³⁰, where different high-order moment models were compared for the representation of a shock in a monoatomic rarefied gas. However, in low-temperature plasmas, the sources of non-equilibrium and the numerical challenges for the moment models are different to a shock in monoatomic rarefied gases. In the shock test case of Laplante and Groth, the non-equilibrium is due to a gradient that is smaller than the characteristic particle-like collisional mean-free path, which leads to
Comparison of high-order moment models for ions

smooth shock jumps. Furthermore, the hypersonic velocities are a challenge for the moment
models in order to correctly capture the smooth jump. Alternatively, in low-temperature
plasmas, due to the weakly-ionized state, the collisions between charged particles are often
negligible as compared to the collisions of the charged particles with the background gas.
As a result, even in collisionally-dominated plasmas at high pressures, the ions can be out
of thermodynamic equilibrium when the plasma is weakly ionized. Another source of non-
equilibrium are the charge-exchange collisions, which is often dominant in the ion collisions,
and introduces an important anisotropy in the collision geometry, not present in the collision
model of Laplante and Groth. Finally, the electric field plays a fundamental role in the non-
equilibrium of low-temperature plasmas. In particular, it is very strong in the sheath region,
and leads to non-Maxwellian distributions with distributions that strongly depends on the
pressure regime. These challenges are very specific to partially-ionized plasmas and will be
studied in this paper.

In our paper, we will compare different high-order moment closures, some of them never
studied for low-temperature plasmas. We will consider: a regularized Grad's moment
closure³¹, an interpolative maximum entropy closure³²⁻³⁴, and two QMOM models: the hy-
perbolic QMOM (HyQMOM)³⁵ and the extended QMOM (EQMOM)³⁶. For each method,
we will solve up to the fourth-order moment, i.e., five scalar moments (5M) in 1D. As some
of the closures do not allow for an analytical integration of the collisional processes, we
will study a simplified collisional model, considering only ionization and charge-exchange
collisions. We will propose an efficient numerical resolution of the moment equations that is
based on analytical formulas or approximations of the information required by the numeri-
cal scheme (such as the eigenvalues of the flux Jacobian or the closing flux), minimizing the
computational cost, which remains comparable to classical 3M models.

The paper is organized as follows. In Section II, we will introduce the setup representing a
low-pressure plasma plasma between two floating walls. Then, we will present the underlying
kinetic model for the ions. After, we will present the moment equations as well as the four
considered high-order moment closures, focusing on the theoretical aspects and requirements
for a robust moment set of equations. In Section III, we will present the numerical scheme
for the different set of equations. In particular, we will provide analytical and approximated
expressions for the numerical fluxes, which allows to have very efficient numerical solvers.
Finally, in Section IV, we will present the numerical results of the different moment closures
Comparison of high-order moment models for ions

representing the ion dynamics in an argon plasma between two floating walls. We will
compare the solutions of the moment models to kinetic simulations, comparing both the
representation of the moments and the reconstructed VDF in the bulk and the sheath. The
paper will finish with a discussion of the findings and a summary of the conclusions in
Section V.

II. KINETIC AND MOMENT MODELS

A. Numerical setup and kinetic simulations

We consider a one-dimensional partially-ionized plasma between two floating walls.
The numerical set-up mimics an electropositive noble gas discharge at different pressure
regimes<sup>2,3</sup>. A scheme of the setup is presented in Fig. 1. We show the ion and electron
densities (in black and orange, respectively), along with the electrical potential (in green).
The system exhibit two regions with very different dynamics: the bulk, in the center of the
discharge, where the potential is mostly flat and densities of ions and electrons are nearly
equal; and the sheath, in the vicinity of the wall, where the plasma is non-neutral (positively
charged), leading to a sharp drop of the electrical potential that will accelerate the ions to
very high velocity towards the wall.

<figure><img src="image_2.png" /><figcaption>FIG. 1. Scheme of the numerical set-up.</figcaption></figure>

The kinetic simulations that will be used as a reference for the high-order moment model
are based on the particle-in-cell Monte-Carlo collisions (PIC-MCC) method. We study a
domain that extends from $x \in [-L, L]$, where the length $L = 5$ cm. The electrons and
Comparison of high-order moment models for ions

singly-ionized ions are simulated with the PIC-MCC method whereas the neutral gas is
considered to be a spatially-homogeneous background at different pressures with constant
temperature of T_g = 300 K. The plasma species are initialized with a constant profile at
densities n_e(t = 0) = n_i(t = 0) = 10^15 m^{-3} and T_e = 5 eV and T_i = T_g, where the subscripts
e, i and g refer to electrons, ions, and gas, respectively. The charged species are absorbed
at the boundaries of the simulation domain, where the electric potential is imposed to be
φ(x = -L) = φ(x = L) = 0. In order to balance the loss of charged species at the walls, a
number of electrons and ions that is equal to the number of ions lost at the walls is injected
at every time step with a probability that is proportional to the electron local density,
as explained in Ref.²⁸. The electrons are injected at T_inj (given in Table I) and the ions
at the neutral temperature T_g, both species injected following a Maxwellian distribution.
This particle injection method mimics both the ionization process for the creation of ions
and the electron heating mechanism that sustains the discharge. In this paper, we will
focus on the ion kinetics and the considered injection model is effectively equivalent to an
electron-impact ionization for the ion population. Similar one-dimensional set-ups have
been previously considered by different authors to study the plasma-sheath transition in gas
discharges⁹,¹⁰,²⁸,³⁷.

In our study, we will compare different high-order moment closures for ions to the kinetic
simulations. Some of the moment closures do not allow for analytical integration of the
collisional terms. For this reason we will consider a simplified model for the collisional cross
sections in an argon plasma. The ion collisions with the neutrals will be modeled with
a back-scattering collision with a polarization potential model where the cross section is
proportional to the inverse of the relative velocity, as described in the next section. Coulomb
collisions are neglected in the model. The PIC code used in this work is 1D-3V but, as it will
be shown in the following section, the ion kinetic equation of the above-mentioned set-up
can be reduced to a 1D-1V problem. The electrons carry out elastic collisions with the gas
(ionization and excitation collisions are not included in the model), with the elastic cross
section as found in the Phelps database from LXCat³⁸.

As explained by classical works on the plasma-sheath formation, based on the isothermal
fluid equations with Boltzmann distributed electrons<sup>3,39,40</sup>, the studied plasma equilibrium
depends mainly on the following non-dimensional parameters: The ratio between the Debye
length and the domain length λ<sub>D</sub>/L, with λ<sub>D</sub> = &radic;<span style="text-decoration: overline;">ε<sub>0</sub>T<sub>e</sub>/(n<sub>e</sub>e)</span> where ε<sub>0</sub> is the vacuum
Comparison of high-order moment models for ions

<table><caption>TABLE I. Parameters of the simulations. $N_{ppc}$ is the number of particle per cell, $N_x$ the number of cells, $T_g$ the gas temperature (considered constant and uniform during the simulation), $T_{inj}$ the injection temperature of the electrons, $\lambda_{mpf}$ the mean-free-path of the ions, $L$ the characteristic length of the system (corresponding to the length of half of the domain), $\lambda_D$ the Debye length, and $T_i$ and $T_e$ the ion and electron temperatures.</caption><thead><tr><th>pressure (Pa)</th><th>N<sub>ppc</sub></th><th>N<sub>x</sub></th><th>T<sub>g</sub> (eV)</th><th>T<sub>inj</sub> (eV)</th><th>λ<sub>mpf</sub>/L</th><th>λ<sub>D</sub>/L</th><th>T<sub>i</sub>/T<sub>e</sub></th></tr></thead><tbody><tr><td>0.01</td><td>100</td><td>2000</td><td>0.025</td><td>5</td><td>10</td><td>2·10<sup>-3</sup></td><td>0.04</td></tr><tr><td>0.1</td><td>100</td><td>2000</td><td>0.025</td><td>5</td><td>1</td><td>2·10<sup>-3</sup></td><td>0.06</td></tr><tr><td>1</td><td>100</td><td>2000</td><td>0.025</td><td>5</td><td>0.1</td><td>2·10<sup>-3</sup></td><td>0.05</td></tr><tr><td>10</td><td>200</td><td>2000</td><td>0.025</td><td>10</td><td>0.01</td><td>2·10<sup>-3</sup></td><td>0.07</td></tr></tbody></table>

permittivity, $T_e$ is the electron temperature, $n_e$ is the electron density, and $e$ is the elementary charge; the ratio between the ion mean free path and the domain length $\lambda_{mpf}/L$, with $\lambda_{mpf} = (n_g\sigma_0)^{-1}$ where $n_g$ is the gas density, $\sigma_0 = 10^{-18} \text{ m}^2$ is a reference ion-neutral collision cross section for argon³; and the ratio between the ion and electron temperatures $T_i/T_e$. A summary of the numerical parameters as well as the above-mentioned normalized parameters for the different cases of study are summarized in Table I.

## B. Ion kinetic equation

As mentioned previously, we focus on the ion kinetics. To simplify notation, we will remove the ion index in the rest of the paper. The ion kinetic model is based on the Boltzmann equation, which describes the evolution in time and space of a VDF $f(t, \vec{r}, \vec{v})$ of a particle of mass $m$ and charge $e$, under an external electric field $\vec{E}$, as follows,

$$ \frac{\partial f}{\partial t} + \vec{v} \cdot \nabla_r f + \frac{e\vec{E}}{m} \cdot \nabla_v f = \left. \frac{\delta f}{\delta t} \right|_c = \left. \frac{\delta f}{\delta t} \right|_c^{\text{el}} + \left. \frac{\delta f}{\delta t} \right|_c^{\text{iz}}, \quad (1) $$

where the right-hand side term accounts for the collisional processes, divided into ion-gas elastic collisions and ionization collisions (the self collisions and recombination collisions are neglected as the plasma is considered to be weakly-ionized and therefore these collisions are negligible as compared to the ion-gas collisions). The ion-gas elastic collisions are modeled
Comparison of high-order moment models for ions

with the Boltzmann collision term:

$$ \left. \frac{\delta f}{\delta t} \right|_{\text{c}}^{\text{el}} = \iint [f(\vec{v}') f_{\text{g}}(\vec{v}_{\text{g}}') - f(\vec{v}) f_{\text{g}}(\vec{v}_{\text{g}})] \sigma(|\vec{v} - \vec{v}_{\text{g}}|, \chi) |\vec{v} - \vec{v}_{\text{g}}| \, d^2\Omega \, d^3 v_{\text{g}}, \quad (2) $$

where the primed velocities refer to the velocities of the restitution collision (which conserves momentum, energy and follows micro-reversibility of the collision), $\sigma$ is the differential cross section, $\chi$ is the scattering angle, and $d^2\Omega$ is the unit sphere element of the collision angles.

As explained by Robson et al.⁴¹, the charge-exchange collision can be approximated as purely backscattering, i.e., the scattering angle is $\chi \simeq \pi$, and the differential cross section will have the following expression:

$$ \sigma(|\vec{v} - \vec{v}_{\text{g}}|, \chi) = \frac{\sigma^{(0)}(|\vec{v} - \vec{v}_{\text{g}}|)}{2\pi} \delta(\cos \chi - 1), \quad (3) $$

where $\delta$ is the Dirac delta function and $\sigma^{(0)}(|\vec{v} - \vec{v}_{\text{g}}|)$ is the total cross section.

If we consider a Langevin potential, i.e., a polarization potential, the cross section is proportional to the inverse of the relative speed, and the integrals over the scattering angle and the gas velocity can be performed, as follows,

$$ \left. \frac{\delta f}{\delta t} \right|_{\text{c}}^{\text{el}} = K^{(0)} \int [f(\vec{v}_g) f_g(\vec{v}) - f(\vec{v}) f_g(\vec{v}_g)] \, d^3 v_g = n_g K^{(0)} [n(\vec{r}, t) w_g(\vec{v}) - f(\vec{v})] , \quad (4) $$

where $K^{(0)}$ is a collision rate, such that $\sigma^{(0)}(|\vec{v} - \vec{v}_{\text{g}}|) = \frac{K^{(0)}}{|v - v_{\text{g}}|}$. In this work, we choose $K^{(0)} = 4.6 \cdot 10^{-16} \, \text{m}^3 \cdot \text{s}^{-1}$ so that the simplified cross section has a similar magnitude to the actual argon cross section. In Eq. (4), $n = \int f(\vec{v}) d^3 v$ is the ion density, and $w_g(\vec{v}) = f_g(\vec{v})/n_g$. As seen in Eq. (4), the resulting collision operator is equivalent to a Bhatnagar-Gross-Krook (BGK) operator.

The ionization collision is expressed as follows:

$$ \left. \frac{\delta f}{\delta t} \right|_{\text{c}}^{\text{iz}} = n_g K_{\text{iz}}^{(0)} n_e(\vec{r}, t) w_g(\vec{v}), \quad (5) $$

where the ionization rate is such that, at each time step, the number of particles added to the domain is equal to the number of particles that left the domain at the boundaries. In other words,

$$ K_{\text{iz}}^{(0)} = \frac{\Gamma_{\text{i,wall}}(x = -L) + \Gamma_{\text{i,wall}}(x = L)}{n_g \int_{-L}^{L} n_e(x) \, dx}, \quad (6) $$

where $\Gamma_{\text{i,wall}}$ is the absolute value of the ion flux at each wall.

In the resulting kinetic model, due to the angular geometry of the back-scattering collision and the 1D geometry, the velocities of the ions in the perpendicular direction are decorrelated
and hence the kinetic model can be reduced to a 1D-1V problem, as the distribution function is in equilibrium at the gas temperature in the transverse direction. In the following, to simplify notation, we will denote $\vec{v} = v\vec{e}_x + v_\perp\vec{e}_\perp$. In that same spirit, we define the one-dimensional VDF that depends on the velocity along $x$, $f(v)$, defined as

$$f(v) = \int_{\mathbb{R}^2} f^{3D}(\vec{v}) dv_y dv_z. \qquad (7)$$

### C. One-dimensional high-order moment equations

The moment equations are obtained by taking averages over the velocity space of the kinetic equation, i.e., Eq. (1). In the high-order moment theory, it is useful to use different definitions for the moments, as it helps for the notation and analysis of the results. First, we define the velocity moments of order $n$ as,

$$M^{(n)} = \int_{-\infty}^{\infty} mv^n f(v) dv. \qquad (8)$$

The first moments are: $M^{(0)} = \rho = mn$ the mass density, and $M^{(1)} = \rho u$ the mass density flux with $u$ the drift velocity. Secondly, we define the centered moments as

$$P^{(n)} = \int_{-\infty}^{\infty} mc^n f(c) dc \quad \text{with} \quad c = v - u. \qquad (9)$$

Finally, in the moment theory it is useful to define the normalized (or standarized) moments as follows:

$$P_{\star}^{(n)} = \frac{mv_{\text{th}}}{\rho} \int_{-\infty}^{\infty} c_{\star}^{n} f(c_{\star}) \mathrm{d}c_{\star} = \frac{P^{(n)}}{\rho v_{\text{th}}^{n}} \quad \text{with} \quad c_{\star} = \frac{v-u}{v_{\text{th}}}, \qquad (10)$$

where the thermal velocity is defined as $v_{\text{th}} = \sqrt{p/\rho}$, with $p = P^{(2)}$ being the pressure.

In our study, we will work with 5M closures, which solve the equations up to the fourth-order moment (for $n = \{0, \dots, 4\}$ in Eq.(8)) and the fifth-order moment is the closing flux, needed to truncate the moment hierarchy. Indeed, we found that the fifth-order moment as seen in PIC-MCC simulations depends on the third and fourth order moment in a manner that was independent of the pressure$^{42}$. Hence, there was a possibility to express the fifth-order moment as a function of the previous moments, which motivates one to consider it as the closing moment.

We define the 5M centered moments as:

$$P^{(2)} \equiv p, \quad P^{(3)} \equiv q, \quad P^{(4)} \equiv r \quad \text{and} \quad P^{(5)} \equiv s, \qquad (11)$$
Comparison of high-order moment models for ions

respectively, the pressure, the heat flux (or skewness), the kurtosis, and the hyperskewness$^{43}$. Note that as we are considering a 1D-1V problem, the pressure represents the xx component of the pressure tensor, the heat-flux the xxx component of the third-rank heat-flux tensor and so forth.

The standardized first six moments are denoted as:

$$P_{\star}^{(0)} = 1, \quad P_{\star}^{(1)} = 0, \quad P_{\star}^{(2)} = 1, \quad P_{\star}^{(3)} = q_{\star}, \quad P_{\star}^{(4)} = r_{\star}, \quad P_{\star}^{(5)} = s_{\star}. \quad (12)$$

Finally, the velocity moments are related to the centered moments by the following relations,

$$M^{(2)} = \rho u^2 + p, \quad (13a)$$

$$M^{(3)} = \rho u^3 + 3up + q, \quad (13b)$$

$$M^{(4)} = \rho u^4 + 6u^2p + 4uq + r, \quad (13c)$$

$$M^{(5)} = \rho u^5 + 10u^3p + 10u^2q + 5ru + s. \quad (13d)$$

The moment equations are obtained by taking moments of the kinetic equation, Eq. (1), which yields the following non-linear set of equations,

$$\frac{\partial \rho}{\partial t} + \frac{\partial(\rho u)}{\partial x} = S_{iz} \quad (14a)$$

$$\frac{\partial(\rho u)}{\partial t} + \frac{\partial M^{(2)}}{\partial x} = \frac{eE}{m}\rho + C^{(1)} \quad (14b)$$

$$\frac{\partial M^{(2)}}{\partial t} + \frac{\partial M^{(3)}}{\partial x} = 2\frac{eE}{m}\rho u + C^{(2)} + S_{iz}\frac{k_B T_g}{m} \quad (14c)$$

$$\frac{\partial M^{(3)}}{\partial t} + \frac{\partial M^{(4)}}{\partial x} = 3\frac{eE}{m}M^{(2)} + C^{(3)} \quad (14d)$$

$$\frac{\partial M^{(4)}}{\partial t} + \frac{\partial M^{(5)}}{\partial x} = 4\frac{eE}{m}M^{(3)} + C^{(4)} + 3S_{iz}\left(\frac{k_B T_g}{m}\right)^2 \quad (14e)$$

where $C^{(n)}$ are the collision terms due to elastic collisions and $S_{iz} = mn_g n_e K_{iz}^{(0)}$ is the rate of ion density creation by ionization.

For the elastic collision terms, we take the moments of Eq. (4), which yields,

$$C^{(1)} = -\frac{K^{(0)}}{m}\rho_{i}\rho_{g}u_{i}, \quad C^{(2)} = -\frac{K^{(0)}}{m}\left(\rho_{g}M_{i}^{(2)} - \rho_{i}p_{g}\right),$$

$$C^{(3)} = -\frac{K^{(0)}}{m}\rho_{g}M_{i}^{(3)}, \quad C^{(4)} = -\frac{K^{(0)}}{m}\left(\rho_{g}M_{i}^{(4)} - 3\rho_{i}p_{g}^{2}\right). \quad (15)$$
## D. Analytical and interpolative 5M closures

The system of equations (14) is not closed as the moment $M^{(5)}$ (or more precisely the centered moment $s$, also known as hyper skewness) needs to be written as a function of the other "known" moments. This is the so-called closure problem in the moment theory. The usual way to do so is to consider a VDF expression that depends on five parameters. In general, these five parameters depend non-linearly on the first five moments and require a numerical inversion of a system of equations. As a result, once the numerical system has been inverted, the VDF is obtained and used to compute the closing moment. The mathematical expression for the VDF as a function of the moments should respect several properties in order to be used as a moment closure: (i) since the VDF represents a probability density, it should be positive in all the velocity domain (here $\mathbb{R}$) for all sets of parameters in their definition space; (ii) the moment system inversion should have a solution for any set of moments in the realizability space (defined as the space of moments that provide a positive distribution function, which is independent of the closure); and (iii) the existence of an entropy inequality of the closed set of equations is important for the stability of the system of equations and the uniqueness of the solution in discontinuities$^{44}$. From a numerical point of view, the system should be hyperbolic (which is related to the existence of an entropy$^{44}$), otherwise the system will have complex (non-real) eigenvalues, which can lead to nonphysical behavior. Note that all these properties are not necessarily satisfied by all the common moment closures.

In this study, we consider four of the most common closures for 5M equations, where we provide analytical expressions or interpolative approximations of the closure.

### 1. Globally hyperbolic regularized Grad closure

The distribution function of Grad's method of moment is a Maxwellian perturbed by a polynomial in the velocity, which can be conveniently expressed as function of the Hermite polynomials$^{45}$, as follows,

$$f_{\text{Grad}}(c) = \mathcal{M}(c; \rho, T) \left( 1 + \sum_{k=0}^{N-1} h_k H^{(k)}(c) \right), \quad (16)$$
Comparison of high-order moment models for ions

where $M(c; \rho, T)$ is a centered Maxwellian of density $\rho$ and temperature $T$, $H^{(k)}$ is the Hermite polynomial of order $k$, and $h_k$ is a coefficient that is a linear combinations of the standarized moments. The inversion is therefore analytical and straightforward, due to the orthogonality of the Hermite polynomials$^{46}$. In the 5M case, the Grad distribution function reads

$$f_{\text{Grad}}^{5\text{M}}(c_{\star}) = \frac{\rho}{mv_{\text{th}}} M(c_{\star}) \left( 1 - \frac{q_{\star}}{6} (3c_{\star} - c_{\star}^3) + \frac{r_{\star} - 3}{24} (3 - 6c_{\star}^2 + c_{\star}^4) \right), \quad (17)$$

where the normalized Maxwellian distribution (or Gaussian distribution) is defined as,

$$M_{\star}(c_{\star}) = \frac{e^{-c_{\star}^2/2}}{\sqrt{2\pi}}. \quad (18)$$

The standardized closing flux can be easily computed as

$$s_{\star} = \frac{mv_{\text{th}}}{\rho} \int_{-\infty}^{\infty} c_{\star}^{5} f_{\text{Grad}}^{5\text{M}}(c_{\star}) dc_{\star} = 10q_{\star}. \quad (19)$$

In Fig. 2, the values of the closing flux is presented in the realizability space and compared to other closures. Despite the simplicity of Grad’s distribution function, it has major drawbacks in its applicability to conditions far from equilibrium. First, as seen in Eq. (17), the distribution function is not guaranteed to be non-negative, which can lead to negative populations in the tail of the VDF. In addition, when the closing flux of Eq. (19) is used in the system of equations (14), the resulting system has mathematical problems related to the hyperbolicity of the equations. In this paper, we will use a fix to this problem, proposed by Cai et al.$^{31}$, which consists in adding an additional term to Eq. (14e), allowing to recover a hyperbolic system of equations. As a result, in the globally hyperbolic regularized Grad model, Eq. (14e) becomes,

$$\frac{\partial M^{(4)}}{\partial t} + \frac{\partial M^{(5)}}{\partial x} + 10 \frac{pq}{\rho^2} \frac{\partial \rho}{\partial x} - 5K \frac{\partial u}{\partial x} - 10 \frac{q}{\rho} \frac{\partial p}{\partial x} = 4 \frac{eE}{m} M^{(3)} + C^{(4)} + 3S_{iz} \left( \frac{k_B T_g}{m} \right)^2, \quad (20)$$

where $K = r - 3p^2/\rho$. Note that the regularization does not change the closure (that remains the regular closure of Grad). This is simply a fix for solving the global system (Eq. (14)) without loss of hyperbolicity.

## 2. Interpolative Maximum Entropy

The maximum entropy closure$^{47,48}$ is based on the thermodynamical principle of entropy maximization. We underline that the maximum-entropy principle differs in multi-component
Comparison of high-order moment models for ions

mixtures as compared to monoatomic non-reacting gases, as explained in Ref.⁴⁹. Nevertheless, in this work, we take the simple monoatomic gas maximum-entropy principle, as done previously²¹,²². For a given number of moments, the VDF that maximizes the entropy is an exponential of a polynomial in the velocity. In our 5M case, the VDF will thus take the form

$$f_{ME}^{1D}(c_{\star}) = \exp\left(\sum_{i=0}^{4} k_i c_{\star}^i\right). \qquad (21)$$

The moment inversion of this closure has no analytical expression and its numerical resolution is stiff and, in general, badly conditioned when far from equilibrium. In the results (Sec. IV), we will show a comparison between the VDF as reconstructed by the maximum entropy closure and the PIC simulations. In order to compute the maximum entropy distribution from the moments, we solve an optimization problem, as explained in Appendix A.

For the computation of the closing flux, however, McDonald and Torrilhon³² developed an interpolative formula that approximates the exact inversion with very high accuracy. This interpolative formula can be written as

$$s_{\star} = \frac{q_{\star}^{3}}{\beta^{2}} + \left(10 - 8\sqrt{\beta}\right) q_{\star}, \qquad (22)$$

with

$$\beta = \frac{1}{4} \left( 3 - r_{\star} + \sqrt{(3 - r_{\star})^2 + 8q_{\star}^2} \right). \qquad (23)$$

The variation of $s_{\star}$ as a function of $q_{\star}$ and $r_{\star}$ is shown in Fig. 2 along with the other closures.

This closure has the disadvantage of presenting a singular subspace in the realizability space⁵⁰, referred to as the Junk line. Indeed, for $r_{\star} > 3$, if $q_{\star} = 0$ (which corresponds to the vertical semi-line above the equilibrium in Fig. 2), the closing moment and the VDF are not defined. Furthermore, on each side of this singular line, the closing hyperskewness has opposite signs and are diverging towards $\pm\infty$ when approaching the singularity, and so will do the eigenvalues. Numerically, as proposed by Refs.²²,³², one can prevent the closing flux of Eq. (22), becoming singular as $\beta \to 0$. These references propose to limit the value of $\beta$ around the singularity. In our simulations, we fixed it to $\beta_{min} = 10^{-4}$, similar to what was done in Refs.²²,³². The value of the cutoff $\beta_{min}$ did not seem to have a significant impact on the profiles, if small enough, but the time step becomes extremely small for smaller values of
Comparison of high-order moment models for ions

$\beta_{min}$ if some points of the simulation are close to the Junk line. As will be seen in section IV, under low pressure conditions the presence of the singular Junk line becomes a problem.

### 3. Hyperbolic quadrature method of moments (HyQMOM)

The quadrature method of moments (QMOM) considers a VDF that is a sum of Dirac delta distributions (that represent the quadrature points to approximate integrals over the distribution)$^{51}$, formally,

$$f_{QMOM}(c) = \sum_{k=0}^{N} w_k \delta(c - c_k). \qquad (24)$$

where $w_k$ are the weights of each Dirac delta distribution and $c_k$ is the abscissa of the Dirac delta distributions.

The Hyperbolic QMOM (HyQMOM) closure$^{35}$ is a variant of QMOM that adds a condition that ensures the hyperbolicity of the system of equations. This closure will therefore consider $2N - 1$ moments for $N$ Diracs.

In the case of a 5-moments system, we impose the fifth order moment to respect the following relation,

$$s_{\star} = 2r_{\star}q_{\star} - q_{\star}^{3}.$$

The closing flux is presented in Fig. 2 along with the other considered closures.

The VDF will therefore consist of two Dirac distributions with abscissae that depend on the moments, and one Dirac fixed at the center of the VDF, namely at the drift velocity (i.e., $c = 0$). Therefore, in our case, we will consider a VDF of the form

$$f_{\text{HyQMOM}}^{5\text{M}}(c_{\star}) = \frac{\rho}{mv_{\text{th}}} [w_0\delta(c_{\star}) + w_1\delta(c_{\star} - c_1) + w_2\delta(c_{\star} - c_2)]. \qquad (25)$$

where the weights and abscissa can be analytically computed from the moments$^{35}$, as follows,

$$ \begin{aligned} w_0 &= 1 - (w_1 + w_2), & w_1 &= \alpha w_2, & w_2 &= \frac{\alpha}{(1 + \alpha)c_1^2} \\ c_1 &= \sqrt{\frac{r_{\star}}{1 - \alpha + \alpha^2}}, & c_2 &= -\alpha c_1 \\ \alpha &= \frac{2r_{\star} - q_{\star}^2 - |q_{\star}|\sqrt{4r_{\star} - 3q_{\star}^2}}{2(r_{\star} - q_{\star}^2)}. \end{aligned} \qquad (26) $$

Note that, in general, the QMOM closure for more than 5M does not allow for analytical expressions of the closing flux nor the VDF and they need to be computed through an eigenvalue problem with the Chebychev algorithm$^{52}$.
4. Extended QMOM (EQMOM)

The Extended QMOM (EQMOM) closure³⁶ is similar to the QMOM closure where the Dirac distributions are replaced by continuous functions, such as Gaussians with a single temperature.

In the 5M case, the EQMOM VDF will therefore have the form,

$$f_{\text{EQMOM}}^{\text{5M}}(c_{\star}) = \frac{n}{v_{\text{th}}} \left[ \rho_{1}^{*} \mathcal{M}_{\star} \left( \frac{c_{\star} - c_{1}^{*}}{a_{\star}} \right) + \rho_{2}^{*} \mathcal{M}_{\star} \left( \frac{c_{\star} - c_{2}^{*}}{a_{\star}} \right) \right], \quad (27)$$

where the normalized Maxwellians are defined in Eq. (18) and the normalized parameters are computed from the standardized moments as follows:

Denoting $b_{\star} = 1 - a_{\star}^2$, the parameters are defined as

$$\begin{align}
c_1^* &= -\sqrt{\frac{\rho_2^*}{\rho_1^*} b_*}, & c_2^* &= \sqrt{\frac{\rho_1^*}{\rho_2^*} b_*}, \nonumber \\
\rho_1^* &= \frac{1}{2} + \frac{q_*}{2\sqrt{q_*^2 + 4b_*^3}}, & \rho_2^* &= 1 - \rho_1^*. \tag{28}
\end{align}$$

and $b_{\star}$ is a solution of the equation

$$b_{\star}^{3} + \frac{r_{\star} - 3}{2} b_{\star} - \frac{(q_{\star})^{2}}{2} = 0, \quad (29)$$

that gives, using the third order polynomial formula

$$b_{\star} = C_{3} - \frac{C_{2}}{C_{3}} \quad (30)$$

$$\text{with } C_3 = (C_1 + (C_1^2 + C_2^3)^{1/2})^{1/3}, \quad C_1 = \frac{q_*^2}{4} \quad \text{and} \quad C_2 = \frac{r_* - 3}{6}, \quad (31)$$

where the roots must be treated as complex principal root.

Finally, this closure allows for an analytical closing flux, written as,

$$s_{\star} = \frac{q_{\star}^{3}}{b_{\star}^{2}} + (10 - 8b_{\star})q_{\star}. \quad (32)$$

Like the maximum entropy closure, EQMOM has the same singular subspace, as shown in Fig. 2. Both closures are in fact very similar. The advantage of EQMOM is its mathematical simplicity compared to maximum entropy that requires the inversion of an integral system. However, as will be seen in Sec. IV, the latter is better at capturing a wide variety of VDF forms. The singular subspace problem is addressed using the same approach as in maximum entropy.
<figure><img src="image_3.png" /><figcaption>FIG. 2. Variation of the closing flux $s_{\star}$ as the heat flux $q_{\star}$ and the kurtosis $r_{\star}$. The hatched zoned represent the non-realizable region where no positive function can have such heat-flux and kurtosis, and the dashed line represents the Junk line and the purple cross the thermodynamic equilibrium.</figcaption></figure>

In Fig. 2, we compare the standarized hyperskewness for the different models in the realizability space of the 5M closures, where we represented the singular subspace for EQMOM and maximum entropy (Max. Entr.). We can clearly see that these two closures are indeed almost identical. Alternatively, Grad closure shows a great simplicity as the hyper skewness only depends on the heat flux. Finally, the HyQMOM closure has no singularities and a behavior that resembles Grad around the equilibrium, but that also depends on the kurtosis far from equilibrium.

### III. NUMERICAL METHODS

The presented moment closures are non-linear system of equations that can be written in the form:

$$ \frac{\partial \mathbf{M}}{\partial t} + \frac{\partial \mathbf{F}(\mathbf{M})}{\partial x} + \mathbf{B}(\mathbf{M}) \frac{\partial \mathbf{M}}{\partial x} = \mathbf{S}(\mathbf{M}), \qquad (33) $$

where $\mathbf{M}$ are the velocity-moment variables, $\mathbf{F}(\mathbf{M})$ are the fluxes, $\mathbf{B}(\mathbf{M})$ is the non-conservative product matrix (which is $\mathbf{B}(\mathbf{M}) = 0$ for all the closures except for regularized Grad), and $\mathbf{S}(\mathbf{M})$ is the source term that contains the electric field and the collisional terms, the notations for the numerical scheme are given in Appendix B. The numerical resolution of Eq. (33) presents two main difficulties: (i) in high-order moment closures it is often difficult to have analytical approximations for the numerical fluxes, as required for computationally-efficient schemes, and (ii) the numerical treatment of this non-conservative product is technically complex and requires specific numerical schemes. In the following,
Comparison of high-order moment models for ions

we present a second-order numerical scheme that aims at providing efficient and stable solutions.

## A. Second-order finite volume scheme for Maximum Entropy, HyQMOM, and EQMOM closures

The system of equations of Eq. (33), with $\mathbf{B}(\mathbf{M}) = 0$, can be discretized with classical finite volume schemes$^{53,54}$. In the finite volume discretization, the domain is divided into $N_C$ cells with size $\Delta x_i$ for $i \in \{1, 2, \dots, N_C\}$, were the time evolution of the cell-averaged conservative variables in the cell i $\mathbf{M}_i$ follows the equation,

$$ \frac{d\mathbf{M}_i}{dt} + \frac{1}{\Delta x_i} (\mathbf{F}_{i+1/2} - \mathbf{F}_{i-1/2}) = \mathbf{S}_i, \qquad (34) $$

where the fluxes $\mathbf{F}_{i\pm 1/2}$ are evaluated in the cell interface $i \pm 1/2$. For the time discretization, in this work we use a simple forward Euler scheme, as we are interested in steady-state solutions and the convergence times are relatively fast. As a result, the discretized equation reads

$$ \mathbf{M}_i^{n+1} = \mathbf{M}_i^n - \frac{\Delta t}{\Delta x_i} (\mathbf{F}_{i+1/2}^n - \mathbf{F}_{i-1/2}^n) + \Delta t \mathbf{S}_i^n. \qquad (35) $$

Here, the cell-averaged source term is approximated with the cell averages variables, i.e., $\mathbf{S}_i^n = \mathbf{S}(\mathbf{M}_i^n)$. Alternatively, in order to compute the fluxes at the cell interface, the considered closures have the advantage that the homogeneous system of equation (i.e., with $\mathbf{S}(\mathbf{M}) = 0$), is a hyperbolic non-linear system of equations (i.e., the Jacobian matrix $\partial \mathbf{F}/\partial \mathbf{M}$ is diagonizable with real eigenvalues). As a result, Godunov schemes and approximated Riemann solvers can be applied to the discretization of the fluxes. In the present work, we use a Rusanov scheme that reads as follows,

$$ \mathbf{F}_{i+1/2}^n = \frac{1}{2} [\mathbf{F}(\mathbf{M}_L^n) + \mathbf{F}(\mathbf{M}_R^n)] - \frac{|\lambda_{i+1/2}^{max}|}{2} (\mathbf{M}_R^n - \mathbf{M}_L^n), \qquad (36) $$

where the sub-index $R/L$ represents the reconstructed values on the right and left of the interface and $|\lambda_{i+1/2}^{max}|$ is the local maximum spectral radius of the Jacobian matrices with the left and right states. The reconstructed values of the conservative variables are computed with a Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL) scheme$^{55}$ with a minmod limiter. We reconstruct the primitive variables $\mathbf{P}$ (defined in Appendix B) and compute the conservative variables at the interface from the reconstructed primitive
Comparison of high-order moment models for ions

variables. We stress that this choice in the reconstructed variables plays a very important role in the stability of the solution. We also tested the reconstruction of the conservative variables, which provides unstable solutions, in particular in the sheath region.

The numerical scheme presented in Eq. (35) with (36) has the following stability condition

$$CFL = \min \left( \Delta t \max_{i \in N_c} \left( \frac{|\lambda_i^{max}|}{\Delta x_i} \right), \Delta t \nu_{ig} \right) < 1. \quad (37)$$

## B. Second-order finite volume scheme for the regularized Grad closure

The regularization of Grad's model is proposed by Cai et al.³¹ such that the matrix $\partial F/\partial M + B$ is hyperbolic. This method is an elegant solution to the loss of hyperbolicity of original Grad's moment method but it comes with the drawback of adding a non-conservative product to the equations. A consequence of the non-conservative character of the fluxes is that classical Godunov's schemes cannot be applied as there are no general Rankine-Hugoniot relations. A solution to this difficulty was proposed by Dal Maso et al.⁵⁶ with a generalization of the weak solution introducing the notion of a “path” in the phase space of the variables⁵⁷. This “path” adds a degree of freedom to the numerical method. In the present work, we adapt the work of Dumbster & Balsara⁵⁸ to our system of equations. The discretization of Eq. (33) is done by writing the fluctuation form of the discretized system, as follows,

$$M_i^{n+1} = M_i^n - \frac{\Delta t}{\Delta x_i} (D_{i+1/2}^{-} + D_{i-1/2}^{+}) - \frac{\Delta t}{\Delta x_i} (F_{i+1/2}^{-} - F_{i-1/2}^{+}) - \frac{\Delta t}{\Delta x_i} \hat{B}(P_i^n) \Delta P_i^n + \Delta t S_i^n, \quad (38)$$

where the details of the discretization of the fluxes and the non-conservative matrix are given in Appendix B. The stability condition of this scheme is the same as in Eq. (37).

## C. Analytical and approximated expressions of the eigenvalues of the Jacobian of the flux

In order to estimate $|\lambda_{i+1/2}^{max}|$, needed for the discretization of the fluxes, we need to evaluate the eigenvalues of the Jacobian matrix of the flux. This is one of the major difficulties in high-order moment models as they, in general, do not have analytical expressions. The numerical computation of $|\lambda_{i+1/2}^{max}|$ may result in computationally expensive numerical calculation of the Jacobian matrix and its eigenvalues. In order to avoid this, in this work,
Comparison of high-order moment models for ions

we provide analytical expressions (when possible) and interpolative expressions in order to
have computationally efficient and accurate algorithms.

a. HyQMOM: The eigenvalues of the HyQMOM Jacobian flux can be analytically computed and read as follows³⁵,

$$
\lambda_{0,1,3,4}^{\text{HyQMOM}} = u + \sqrt{\frac{p}{\rho}} \left( \frac{q_*}{2} \pm \sqrt{1 + \Upsilon + \frac{q_*^2}{4} \pm \sqrt{\Upsilon(1+\Upsilon)}} \right), \quad \text{and} \quad \lambda_2^{\text{HyQMOM}} = u \quad (39)
$$

where $\Upsilon = r_{\star} - q_{\star}^{2} - 1$ can be proven to be a non-negative quantity for any realizable set of moments and, hence, the eigenvalues are indeed always real.

b. *Maximum entropy*: The eigenvalues of the maximum entropy closure cannot be computed analytically. In this work, we use the approximation proposed by Baradaran³³ (c.f. Boccelli³⁴ for the 3D case), which we give for completeness, as follows,

$$
\lambda_m^{\text{Max.Ent.}} \approx u + v_{\text{th}} \lambda_m^* \quad \text{with} \quad m \in \{0, 1, 2, 3, 4\} \qquad (40)
$$

and

$$
\begin{align}
\lambda_{0,4}^* &= \frac{q_* \pm \sqrt{q_*^2 \mp \frac{4}{5} q_*^2 \beta C_* + 4 \beta^2 Y_*}}{2\beta} \pm \frac{8}{10} C_*, & \lambda_{1,3}^* &= \frac{q_* \pm \sqrt{q_*^2 \mp \frac{6}{5} q_*^2 \beta C_* + 4 \beta^2 X_*}}{2\beta} \mp \frac{3}{10} C_* \notag \\
\lambda_2^* &= \frac{2q_*(q_*^2\sqrt{\beta+2\beta^3})}{\sqrt{\beta^3(q_*^2+2\beta^2)}} - \lambda_0^* - \lambda_1^* - \lambda_3^* - \lambda_4^*, & & \tag{41}
\end{align}
$$

where

$$
\begin{align}
X_{\star} &= A_{\star} + \left(\frac{3}{10}\right)^2 C_{\star}^2 + \frac{3}{5} C_{\star} \sqrt{A_{\star}}, & Y_{\star} &= B_{\star} + \left(\frac{8}{10}\right)^2 C_{\star}^2 - \frac{8}{5} C_{\star} \sqrt{B_{\star}}, & C_{\star} &= \sqrt{3-3\beta}, \nonumber \\
A_{\star} &= 5 - 4\sqrt{\beta} - \sqrt{10-16\sqrt{\beta}+6\beta}, & B_{\star} &= 5 + 4\sqrt{\beta} - \sqrt{10-16\sqrt{\beta}+6\beta}. \tag{42}
\end{align}
$$

c. **EQMOM**: The eigenvalues of EQMOM cannot be computed analytically. In the present work, we provide an interpolative approximation of the spectral radius, similar to the method of Baradaran³³ for maximum entropy shown above. The asymptotes of the highest eigenvalues at positive and negative high heat-flux are respectively $q_*/b_* + \sqrt{3(1-b_*)}$ and $\sqrt{3(1-b_*)}$. The asymptotes of the lowest eigenvalues at positive and negative high heat-flux are respectively $-\sqrt{3(1-b_*)}$ and $q_*/b_* - \sqrt{3(1-b_*)}$. As a result, the spectral radius can be estimated with an interpolation between these two asymptotes, as follows,

$$
|\lambda_{\text{max}}^{\star, \text{HyQMOM}}| \approx \max \left( \frac{1}{2} \left| \frac{q_{\star}}{b_{\star}} + \sqrt{4 + \left(\frac{q_{\star}}{b_{\star}}\right)^2} + \sqrt{3(1-b_{\star})} \right|, \frac{1}{2} \left| \frac{q_{\star}}{b_{\star}} - \sqrt{4 + \left(\frac{q_{\star}}{b_{\star}}\right)^2} - \sqrt{3(1-b_{\star})} \right| \right). \quad (43)
$$

This fit is able to approximate the spectral radius $\lambda_{\text{max}}^{\text{HyQMOM}} \approx u + v_{\text{th}}|\lambda_{\text{max}}^{\star, \text{HyQMOM}}|$ within
a 10% accuracy in the worst case, which is shown to be valid for the studied cases.
d. Globally hyperbolic regularized Grad The advantage of the regularization proposed by Cai et al.³¹ is that the system has an analytical solution for the eigenvalues at arbitrary moment order. In our case, they read as follows,

$$ \lambda_{0,1,3,4}^{\text{Reg. Grad}} = u \pm v_{\text{th}} \sqrt{5 \pm \sqrt{10}} \quad \text{and} \quad \lambda_{2}^{\text{Reg. Grad}} = u. \qquad (44) $$

## D. Ionization and electric field

In this work, we focus on the influence of the different moment models on the ion dynamics at different pressure regimes. Following the methodology proposed by Bocelli²¹,⁵⁹, we will impose the quantities that depend on the electron dynamics. These are the ionization rate $S_{iz}$ and the electric field $E$, that will be taken from the PIC results and the values will be interpolated to the mesh used for the moment models. In this way, the accuracy of the different ion moment models can be assessed, independently of the model used for the electrons. A self-consistent simulation where the two charged species are coupled is also possible, as previously done in Refs.²⁸,³⁷, and is left for future work comparing electron moment closures.

## E. Numerical grid, boundary and initial conditions

The simulations are run on a non-uniform grid with a cell size of 0.5 mm in the bulk that decrease to 10 µm in the sheath. The simulations are initialized with a constant background of density $n = 10^{15}$ m³ with a Maxwellian distribution at rest with the gas temperature $T = 0.025$ eV. The boundary conditions are imposed with a ghost cell strategy with outflow conditions, where we assume that the ion flow is supersonic at the boundary. The CFL was set to 0.1 at the beginning due to the large variations during the transient and can then be increased to nearly 1.

## IV. RESULTS

We present here the results of the high-order moment simulations representing the ions in the numerical setup of an argon plasma between two floating walls (described in Section II A) and compare them with the PIC-MCC simulations with the conditions summarized
Comparison of high-order moment models for ions

in Table I. As discussed earlier, the simulations study different pressure regimes. We vary the pressure of the background gas from $p_g = 10$ Pa, that corresponds to an ion-neutral mean free path much smaller than the domain, i.e., $\lambda_{mfp}/L = 0.01$, to $p_g = 0.01$ Pa, that corresponds to an ion-neutral mean free path larger than the domain, i.e., $\lambda_{mfp}/L = 10$. We will therefore refer to these two pressures as respectively high and low pressure in the rest of the paper, and the other two pressures as intermediate pressures.

In order to decrease the noise of the PIC data, the results have been averaged over $3.6 \cdot 10^6$ time steps, after convergence which takes of the order of $10^7$ time steps. The high-order moment equations are solved in time until convergence, which requires of the order of $10^5$ to $10^6$ time steps (depending on the pressure regime and the closure). Furthermore, moment results do not require to be time averaged as they are free of statistical noise. The high-order moments models will also be compared to several classical fluid models. One based on the isothermal fluid equations (i.e. 2M) and two non-isothermal fluid models (with 3M) based on different closures: the isotropic Maxwellian and an anisotropic Maxwellian (i.e. with two temperatures, one in the direction $x$ and another in the perpendicular direction considered to be equal to the neutral gas temperature).

A. Moments profiles

We first present the moment profiles for every model considered at the studied gas pressures. Due to the symmetry of the system, we show only half of the domain for each profile. In Fig. 3, we present the density profiles for the 5M models as well as the 3M model with anisotropic Maxwellian compared to the PIC simulations as well as with the isothermal and 3M isotropic Maxwellian models. As it is shown in the figure, the density profile varies from a dome-like profile (Shottky collisionally dominated regime) at high-pressure to a flat profile with a sharp drop in the sheath (Tonks-Langmuir nearly collisionless regime), at low-pressure. In the density profiles, we can see that, in general, regularized Grad (Reg. Grad) and HyQMOM are able to capture almost perfectly the density profile at all pressures. However, both EQMOM and maximum entropy (Max. Entr.) present a spurious density gradient in the center of the domain at low pressure. This feature is associated with the singularity of the closing flux (Junk line) when the kurtosis is positive and the heat-flux is null. As compared to the isothermal 2M and isotropic Maxwellian 3M moments, the mo-
<figure><img src="image_4.png" /><figcaption>FIG. 3. Density profiles of all the simulations for four pressures.</figcaption></figure>

ment models improve the representation of the density profile in the low and intermediate pressures, whereas they compare similarly at high pressures. Alternatively, the anisotropic Maxwellian 3M model is able to largely improve the density profile, although, as it will be shown later, is not able to correctly predict the higher-order moments, the heat-flux at the wall nor the VDF.

The temperature (the *xx* component of the temperature tensor) profiles are shown in Fig. 4. Note that a zoom in the sheath is presented on the right of each profile. At high pressure, the temperature is constant at *T*<sub>g</sub> except in the sheath region, where it increases because of the effect of the strong electric field and the collisions with the gas. Alternatively, the temperature at low and intermediate pressures is lower than the gas temperature in the
<figure><img src="image_5.png" /><figcaption>FIG. 4. Temperature profiles of all the simulations for four pressures. For each plot we show a zoom of the sheath next to it.</figcaption></figure>

center (because of the strong advection) and larger in the sheath (because of the electric field). As in the density profiles, the high-order moment models, in particular HyQMOM, are able to capture with very high fidelity the temperature profiles both in the bulk and in the sheath at all pressures. We also see the effect of the Junk line at low pressure for EQMOM and maximum entropy, where the temperature decreases slightly more than expected in the center of the domain. Note that the large increase at the wall at low pressures for these closure is a consequence of lack of resolution. This is due to the lack of stability and robustness of these models (because of the Junk line) that prevents to run these simulations with better resolution (because the time step tends to zero due to the presence of the Junk line). Finally, the 3M and 2M classical fluid models are less accurate in most of the domain at all pressures. The 3M anisotropic Maxwellian model is as good as the moment models at high pressure, but loses its accuracy as the pressure decreases. This is because of the lack of heat-flux, which becomes important at low pressures, as it will be shown below. The
<figure><img src="image_6.png" /><figcaption>FIG. 5. Heat flux profiles of all the simulations for four pressures. We represent the Fourier heat flux of the non-isothermal Maxwellian cases in dashed lines. For each plot we show a zoom of the sheath next to it.</figcaption></figure>

isotropic 3M model has the right qualitative behavior but lacks of precision at all pressures.

The real advantage of the high-order moment models is their ability to capture higher-order moments without statistical noise that is characteristic of particle-based kinetic simulations. We show the heat flux profile in Fig. 5. At high pressure, the heat-flux is almost zero in the bulk and it presents a sharp drop in the sheath, whereas it is large in the bulk for low pressures. The four moment models overall capture the kinetic behavior with a high level of accuracy. However, we still see the problem of the Junk line in the center of the domain for EQMOM and maximum entropy. We also see small discrepancies in the sheath for Reg. Grad, which can be explained due to the lack of precision of the moment of order 4. For comparison, we also present the heat-flux obtained with the Fourier approximation$^{60}$

$$Q_{xxx} = \frac{6}{5}q_x = -3\frac{pk_B}{m\nu}\partial_x T, \qquad (45)$$

applied to 3M fluid models with the isotropic Maxwellian and the anisotropic Maxwellian.
Comparison of high-order moment models for ions

We clearly see that the Fourier approximation is not well-suited here and that it almost
always overestimates (in norm) the heat-flux, especially at low pressure. This justifies a
high-order moment approach in order to capture the heat-flux to the wall in low-temperature
plasmas: Fourier law is not able to represent the heat-flux inside the sheath due to the large
deviation from thermodynamic equilibrium.

In figure 6 we show the profiles of the standardize heat-flux and kurtosis represented in the phase-space (the so-called realizability space). At high pressures (upper plots), the VDF is closer to a Maxwellian in the bulk, so the heat flux and kurtosis stay mostly around the equilibrium point ($q_*=0$ and $r_*=3$, represented with a pink cross). At these high pressures, all the moment models capture the kinetic behavior with great accuracy in the bulk. In the sheath, while the moments deviate greatly from equilibrium, all except Reg. Grad still match the kinetic behavior with a very high level of accuracy. In Reg. Grad's simulations, the moments even leave the realizability region at the boundaries, in the upper left plot, as a consequence of a VDF that is not ensured to be positive by this method. On the other hand, at lower pressures (lower plots), HyQMOM is the most accurate overall, despite an offset of kurtosis in the bulk. For EQMOM and Maximum Entropy, we see clearly the problem of the crossing of Junk's line in the center of the domain. These two models exhibit an interesting behavior, where the moments peak in the Junk line after being deviated towards the equilibrium. A similar behavior can be seen in the results presented by Boccelli²². This nonphysical behavior is due to the artificial cutoff of the closing moment around Junk's line to avoid its divergence. Despite this problem in the center, EQMOM and maximum entropy are still highly accurate far from the center of the domain except in the last mesh points of the domain close to the wall, due to the lack of resolution of these two simulations at low pressures as mentioned in the previous paragraph. Finally, Reg. Grad highly overestimates the kurtosis everywhere. This overestimation is due to the regularization terms in the kurtosis equation, as mentioned in the previous paragraph. An interesting point to stress on is that despite the very large discrepancy of Reg. Grad in the kurtosis, this has a minor impact in the previous moments (only small inaccuracies in heat-flux, as discussed above).
<figure><img src="image_7.png" /><figcaption>FIG. 6. High-order moments profiles plotted in the realizability space (heat-flux - kurtosis space) for four pressures.</figcaption></figure>

## B. VDF reconstruction

We will now assess the accuracy of the different closures to reconstruct the VDF from the simulated moments. We stress the fact that, in order to obtain low-noise VDF from PIC simulations, the data were averaged over $3.6 \cdot 10^6$ time steps. Alternatively, the VDFs of the high-order moment models are instantaneous values, obtained with a method that is computationally almost as efficient as classical fluid models.

First, since the HyQMOM VDF consists of a sum of Dirac delta distributions, this model will not recover a continuous VDF. This is a major drawback if one is interested in the VDF
Comparison of high-order moment models for ions

shape, in particular for an accurate computation of the collisions or wall processes (e.g., sputtering). However, the three Diracs represent in fact quadrature points used to calculate integrals of a continuous distribution. In that sense, Fox et al.⁶¹ developed a method, called generalized QMOM (GQMOM), to dramatically increase the number of Diracs while keeping the same first moments. The VDF can therefore be refined as much as desired, making it quasi-continuous. This method will be used in order to represent the reconstructed VDFs of the HyMOM results.

The results of the VDF reconstruction at high pressure ($p_g = 10$ Pa) are shown in Fig. 7, in three different points of the domain, including the bulk, the entrance of the sheath and close to the wall. The VDF is very close to Maxwellian in most of the bulk and can therefore be reproduced with great accuracy by the four 5M models in the bulk as well as the 3M model with anisotropic Maxwellian model. In the entrance of the sheath, the distribution has a positive skewness, which is captured by all the 5M models, and it is beyond the 3M representation. Close to the wall, the VDF is strongly non-Maxwellian, and the global shape of the VDF is recovered by all the models with remarkable accuracy, despite some small differences for EQMOM. We note small differences with the kinetic simulation in order to capture the fast decreasing tail on the high-velocities side.

The intermediate pressure regime ($p_g = 1$ Pa) is presented in Fig. 8. We present the VDF in four points: in the center, in the presheath, in the entrance of the sheath, and close to the wall. The kinetic simulation shows that the VDF is nearly Maxwellian at the center, whereas it is very elongated in the presheath and it is a beam with a long tail in the sheath and close to the wall. As seen in the results, these VDFs are beyond the 3M description with an anisotropic Maxwellian, which provides a particular wrong representation in the sheath and close to the wall. Alternatively, all the 5M models largely improve the 3M results. We can see that the regularized Grad model has positivity problems (i.e., negative tails) in the sheath. Alternatively, maximum entropy provides an excellent representation at this pressure. EQMOM and GQMOM capture the right trend while presenting some differences, but still a remarkable accuracy.

Finally, the low pressure ($p_g = 0.01$ Pa) is presented in Fig. 9. As in the previous case, we present the VDF in four points: in the center, in the presheath, in the entrance of the sheath, and close to the wall. In this case, the VDF in not Maxwellian in any of them. In the center of the domain, it presents heavy tails with a symmetric distribution, whereas in
Comparison of high-order moment models for ions

presheath, sheath and close to the wall, the distribution is very skewed, with a very long tail
towards low energies. As a result, these very non-equilibrium distributions can be seen as a
challenge for the different high-order moment models. As previously discussed, due to the
presence of the Junk line, the moment inversion of maximum entropy and the EQMOM do
not necessarily converge around the center of the domain and for visualization purposes are
set to be Maxwellian. In the rest of the domain, these two closures still manage to capture
the very long tail on the low-energy side of the VDF, but they lack precision in capturing
the shape of the peak of the VDF. Nevertheless, they still give a better approximation
than the anisotropic Maxwellian. Reg. Grad is clearly far from its region of validity in the
whole domain with VDFs that exhibit large negative parts. Alternatively, GQMOM is the
only method that manages to capture a very accurate shape of the VDF in the center of the
domain as well as a good approximation of the VDF in the presheath, sheath and close to the
wall. In summary, the 5M HyQMOM (combined with GQMOM for the VDF reconstruction)
shows a great ability to represent non-equilibrium distributions at all pressure both in the
bulk and the sheath.

V. CONCLUSION AND DISCUSSION

We compared the numerical resolution of some of the most common high-order moment
closures to kinetic simulations in a one-dimensional bounded low-temperature argon plasma.
For each moment closure, we used expressions that allow for an efficient computation of the
closing flux and a stable and accurate numerical scheme that is able to have a computational
time that is comparable to classical 3M fluid models, but with a largely improved accuracy.
The accuracy of each method to recover the profiles of fluid quantities (i.e., the density, the
temperature) is similar for the four 5M models and always more accurate than classical 3M
fluid models, except in the center of the domain at low pressure for EQMOM and maximum
entropy, due to the singular behavior of these closures. In particular, the high-order moment
closures are able to capture the ion temperature evolution in the sheath better than models
assuming anisotropic Maxwellians, showing that the heat-flux in the sheath has an impact
that cannot be modeled with simplified a Fourier-type heath-flux. Indeed, the high-order
moment models are able to capture the heat-flux profile with high accuracy both in the bulk
and the sheath, largely improving the classical Fourier approximation. We also compared
<figure><img src="image_8.png" /><figcaption>FIG. 7. VDF from the kinetic simulations and reconstructed from the 5-moment simulation results for a pressure of 10 Pa. The upper plot shows the ion positions at which the VDF are taken.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 8. VDF from the kinetic simulations and reconstructed from the 5-moment simulation results for a pressure of 1 Pa. The upper plot shows the positions at which the VDF are taken.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 9. VDF from the kinetic simulations and reconstructed from the 5-moment simulation results for a pressure of 0.01 Pa. The upper plot shows the positions at which the VDF are taken.</figcaption></figure>
Comparison of high-order moment models for ions

the ability of each method to reconstruct the VDF from the first five moments, showing
promising results. At high pressure, all 5M closures recover the kinetic VDF with high
fidelity. At low pressures, Reg. Grad presents strong negative tails, EQMOM and maximum
entropy have an improved fidelity as compared to a Maxwellian approximation and GQMOM
presents a good accuracy both in the bulk and the sheath of the domain at all pressures.

We propose a summary of the advantages, drawbacks, and possibility of extension of each
of the considered closures:

* Reg. Grad is the simplest method despite the difficulty associated to the numerical discretization, as shown in this work. In our comparison to the kinetic results, Reg. Grad's closure presents a good accuracy in the representation of moments up to the heat-flux in both the sheath and the bulk. However, it presents problems to capture the kurtosis, due to the error in the closing flux associated to the non-conservative terms of the regularization. As shown in the VDF reconstruction, at low pressures the tail of the VDF presents strong negative parts. In our work we present the 5M model, but extending Reg. Grad to a larger number of moments is straightforward, with a moderate increase of the computational cost. In addition, Grad's closure allows for an analytical integration of the collisional terms for general collisional interactions.

* The maximum entropy closure shows some clear advantages in accuracy at high and intermediate pressures, capturing with great accuracy the transport and the kinetic VDF in all the domain. In the 5M case considered here, we used an interpolative method that shows to be computationally very efficient. However, the presence of the singular Junk line is a real drawback at low pressure, both for the accuracy and the stability of the method. The extension to higher-order moment closures is not straightforward and would require specific interpolative expressions. In addition, there will always be a singular subspace (due to the fact that the last coefficient of the maximum entropy distribution function must necessarily be negative), whose impact is difficult to foresee. Furthermore, another drawback is the high computational cost of the VDF reconstruction. As a result, the computation of the collisional terms for general interaction potentials adds a layer of complexity to this problem since there is no analytical simplification of the collision terms and a numerical integration would be required in a general case, losing numerical efficiency.
* The main advantage of EQMOM is that its closure is very similar to maximum entropy but with a much lower computational cost to reconstruct the VDF. In our results, we show that the reconstructed VDF is, however, slightly less accurate than the maximum entropy closure. Nevertheless, an increase in the number of moments considered is rather simple and one can use some efficient moment inversion algorithms to limit the computational cost. In addition, the computation of the collisional terms might be possible in an analytical manner as the model is based on Maxwellian distributions.

* Finally, HyQMOM was found to be the most efficient, robust, accurate, and stable method that captures with very high fidelity all the simulated moments both in the sheath and the bulk at all the studied pressures. The VDF reconstruction is more challenging, but the GQMOM shows to be a great tool, allowing to capture VDFs with high fidelity to the PIC ones at all pressure regimes and both in the bulk and the sheath. In addition, as compared to PIC they do not have any statistical noise. Furthermore, as in the case of EQMOM, an increase in the number of moments is rather simple and there exists efficient algorithms, such as the Wheeler algorithm$^{62}$, for the VDF reconstruction. In addition, the computation of the collisional terms for general collisional cross sections has analytical expressions.

These results prove that the use of high-order moment hierarchies to simulate ions in a low-temperature plasma is a sound approach and could allow an accurate modeling at a much lower cost than kinetic simulations. In the 5M moment cases shown here, the computational cost is comparable to this of classical 3M fluid models. The collision modeling with general collision cross sections is still a challenge in high-order moment closures. Some of the moment methods studied in this work (e.g., Grad, EQMOM, HyQMOM) allow for affordable collision computations that will be studied in a future work.

## ACKNOWLEDGMENTS

This work is part of a PhD funded by the Institut Polytechnique de Paris (IPP) and EUR PlasmaScience.
# Appendix A: Moment inversion in maximum entropy closures

We solve the moment inversion system by performing the Newton-Raphson algorithm: writing the system as

$$F(U) = M, \qquad (A1)$$

where $U = (k_0, ..., k_4)$ is the vector of parameters of the VDF and $M = (1, 0, 1, q_\star, r_\star)$ the vector of normalized moments, one performs iteratively until convergence the following operation:

$$U_{n+1} = U_n - J^{-1}(F(U_n) - M), \qquad (A2)$$

where $J = \partial F / \partial U$ is the Jacobian of the function $F$.

However, the convergence of this method is not ensured and requires an initial guess very close to the solution for it to convergence in a reasonable number of iterations. For this reason, instead of solving directly equation A1, we solve iteratively the equation

$$F(U_k) = M_0 + (M - M_0) \frac{k}{N}, \qquad (A3)$$

where $M_0 = (1, 0, 1, 0, 3)$ is the equilibrium moment, and $N$ is chosen so that the convergence is fast enough. We found $N = \lceil 10 \max(q_\star, r_\star - 3) \rceil + 1$ to work most of the time.

# Appendix B: Notation for the numerical schemes

## 1. Definition of vectors and matrices

The velocity moments are

$$M = \begin{pmatrix} M^{(0)} \\ M^{(1)} \\ M^{(2)} \\ M^{(3)} \\ M^{(4)} \end{pmatrix} = \begin{pmatrix} \rho \\ \rho u \\ \rho u^2 + p \\ \rho u^3 + 3up + q \\ \rho u^4 + 6u^2p + 4uq + r \end{pmatrix}. \qquad (B1)$$
Comparison of high-order moment models for ions

We define the array of primitive variables, as follows,

$$ \mathbf{P} = \begin{pmatrix} \rho \\ u \\ p \\ q \\ \mathcal{K} \end{pmatrix}, \qquad (B2) $$

where the kurtosis is defined as $\mathcal{K} = r - 3p^2/\rho$.

The conservative flux reads

$$ \mathbf{F}(\mathbf{M}) = \begin{pmatrix} M^{(1)} \\ M^{(2)} \\ M^{(3)} \\ M^{(4)} \\ \rho u^5 + 10u^3p + 10u^2q + 5ru + s \end{pmatrix}, \qquad (B3) $$

where the closing flux $s$ depends on the particular model. The Jacobian matrix of the transformation from conservative to primitive variables reads

$$ \frac{\partial \mathbf{M}}{\partial \mathbf{P}} = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 \\ u & \rho & 0 & 0 & 0 \\ u^2 & 2\rho u & 1 & 0 & 0 \\ u^3 & 3\rho u^2 + 3p & 3u & 1 & 0 \\ u^4 - \frac{3p^2}{\rho^2} & 4\rho u^3 + 12pu + 4q & 6u^2 + 6\frac{p}{\rho} & 4u & 1 \end{pmatrix}. \qquad (B4) $$

The non-linear source term contains the electric field and the collisional terms, as follows,

$$ \mathbf{S}(\mathbf{M}) = \begin{pmatrix} S_{iz} \\ M^{(0)} \frac{eE}{m} + \mathcal{C}^{(1)} \\ 2M^{(1)} \frac{eE}{m} + \mathcal{C}^{(2)} + S_{iz} \frac{k_B T_g}{m} \\ 3M^{(2)} \frac{eE}{m} + \mathcal{C}^{(3)} \\ 4M^{(3)} \frac{eE}{m} + \mathcal{C}^{(4)} + 3S_{iz} \left( \frac{k_B T_g}{m} \right)^2 \end{pmatrix}. \qquad (B5) $$

Finally, the non-conservative product matrix $\mathbf{B}$ is zero for all the models except for the
Comparison of high-order moment models for ions

regularized Grad. In that case, the matrix in conservative variables reads,

$$ \mathbf{B}(\mathbf{M}) = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ \frac{5Ku}{\rho} + \frac{10pu}{\rho^2} - \frac{10qu_x^2}{\rho} - \frac{5K}{\rho} + \frac{20qu}{\rho} - \frac{10q}{\rho} & 0 & 0 & 0 & 0 \end{pmatrix}. \quad (B6) $$

For the numerical method, instead, we use the non-conservative product in primitive variables, where the matrix reads

$$ \hat{\mathbf{B}}(\mathbf{P}) = \mathbf{B}(\mathbf{P}) \frac{\partial \mathbf{M}}{\partial \mathbf{P}} = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ \frac{10pq}{\rho^2} - 5K - \frac{10q}{\rho} & 0 & 0 & 0 & 0 \end{pmatrix}. \quad (B7) $$

## 2. Numerical fluxes in fluctuation form

For the vectors and matrices of Eq. (38), we use the Rusanov-type scheme for the approximate Riemann solver, as follows,

$$ \mathbf{D}^{\pm}(\mathbf{M}_{L}^{n}, \mathbf{M}_{R}^{n}) = \frac{1}{2} \left[ \mathbf{F}_{R}^{n} - \mathbf{F}_{L}^{n} + \tilde{\mathbf{B}}(\mathbf{P}_{L}^{n}, \mathbf{P}_{\star}^{n})(\mathbf{P}_{\star}^{n} - \mathbf{P}_{L}^{n}) + \tilde{\mathbf{B}}(\mathbf{P}_{\star}^{n}, \mathbf{P}_{R}^{n})(\mathbf{P}_{R}^{n} - \mathbf{P}_{\star}^{n}) \right] \pm \frac{|\lambda^{\mathrm{max}}|}{2} (\mathbf{M}_{R}^{n} - \mathbf{M}_{L}^{n}) \quad (B8) $$

where the primitive variables $\mathbf{P}$ and the Roe-type matrix $\tilde{\mathbf{B}}$ are given in Appendix B 1. The intermediate state (with the Rusanov scheme) is computed as follows,

$$ \mathbf{M}_{\star}^{n} = \frac{1}{2|\lambda^{\mathrm{max}}|} \left[ |\lambda^{\mathrm{max}}| (\mathbf{M}_{L}^{n} + \mathbf{M}_{R}^{n}) - (\mathbf{F}_{R}^{n} - \mathbf{F}_{L}^{n}) - \tilde{\mathbf{B}}(\mathbf{M}_{L}^{n}, \mathbf{M}_{R}^{n})(\mathbf{M}_{R}^{n} - \mathbf{M}_{L}^{n}) \right], \quad (B9) $$

and the averaged value of the primitive variables $\mathbf{P}_{\star}^{n}$ are computed from the conservative variables $\mathbf{M}_{\star}^{n}$. Finally, the fluxes at the interface are computed with the reconstructed variables, as

$$ \mathbf{F}_{i+1/2}^{-} = \mathbf{F}(\mathbf{M}_{L}^{n}) \quad \text{and} \quad \mathbf{F}_{i-1/2}^{+} = \mathbf{F}(\mathbf{M}_{R}^{n}). \quad (B10) $$

The slopes of the extrapolated values of the primitive variables are computed as

$$ \Delta \mathbf{P}_{i}^{n} = \minmod (\mathbf{P}_{i+1}^{n} - \mathbf{P}_{i}^{n}, \mathbf{P}_{i}^{n} - \mathbf{P}_{i-1}^{n}). \quad (B11) $$
Comparison of high-order moment models for ions

The left and right states at the $i+1/2$ wall interface are computed with this slope, as in the MUSCL scheme, i.e.,

$$ \mathbf{P}_L^n = \mathbf{P}_i^n + \frac{1}{2} \Delta \mathbf{P}_i^n \quad \text{and} \quad \mathbf{P}_R^n = \mathbf{P}_{i+1}^n - \frac{1}{2} \Delta \mathbf{P}_{i+1}^n. \qquad (B12) $$

The conservative variables at the cell interface $M_{L/R}^n$ are computed from the reconstructed primitive variables of Eq. (B12). The Roe-type matrix of the non-conservative product is easily obtained by doing the arithmetic mean between the left and right states, as follows:

$$ \tilde{\mathbf{B}}(\mathbf{P}_R, \mathbf{P}_L) = \hat{\mathbf{B}}(\bar{\mathbf{P}}) \quad \text{with} \quad \bar{\mathbf{P}} = \frac{1}{2}(\mathbf{P}_R + \mathbf{P}_L). \qquad (B13) $$

## REFERENCES

¹D. Goebel and I. Katz, *Fundamentals of Electric Propulsion: Ion and Hall Thrusters*, JPL Space Science and Technology Series (Wiley, 2008).

²M. A. Lieberman and A. Lichtenberg, *Principles of Plasma Discharges and Materials Processing*, 2nd ed. (Wiley-Interscience, 2005) Chap. 6.

³P. Chabert and N. Braithwaite, *Physics of radio-frequency plasmas* (Cambridge University press, 2011) Chap. 3.

⁴S. M. Starikovskaia, “Plasma assisted ignition and combustion,” *Journal of Physics D: Applied Physics* **39**, R265 (2006).

⁵M. G. Kong, G. Kroesen, G. Morfill, T. Nosenko, T. Shimizu, J. van Dijk, and J. L. Zimmermann, “Plasma medicine: an introductory review,” *New Journal of Physics* **11**, 115012 (2009).

⁶E. Kawamura, V. Vahedi, M. A. Lieberman, and C. K. Birdsall, “Ion energy distributions in rf sheaths; review, analysis and simulation,” *Plasma Sources Science and Technology* **8**, R45 (1999).

⁷W. Schottky, “Diffusionstheorie der positiven säule,” *Phys. Zeits* **25**, 635 (1924).

⁸V. A. Godyak, *Soviet radio frequency discharge research* (Delphic Associates, 1986).

⁹A. Tavant, R. Lucken, A. Bourdon, and P. Chabert, “Non-isothermal sheath model for low pressure plasmas,” *Plasma Sources Science and Technology* **28**, 075007 (2019).

¹⁰L. P. Beving, M. M. Hopkins, and S. D. Baalrud, “How sheath properties change with gas pressure: modeling and simulation,” *Plasma Sources Science and Technology* **31**, 084009 (2022).
<footnote n="11">¹¹V. Vahedi and M. Surendra, “A monte carlo collision model for the particle-in-cell method: applications to argon and oxygen discharges,” Computer Physics Communications 87, 179–198 (1995), particle Simulation Methods.</footnote>

<footnote n="12">¹². L. A. Birdsall, C.K., *Plasma Physics via Computer Simulation* (1st ed.) (CRC Press., 1991).</footnote>

<footnote n="13">¹³E. Sonnendrücker, J. Roche, P. Bertrand, and A. Ghizzo, “The semi-lagrangian method for the numerical resolution of the vlasov equation,” *Journal of Computational Physics* **149**, 201–220 (1999).</footnote>

<footnote n="14">¹⁴I. M. Gamba and S. Rjasanow, “Galerkin–petrov approach for the boltzmann equation,” *Journal of Computational Physics* **366**, 341–365 (2018).</footnote>

<footnote n="15">¹⁵G. Oblapenko, D. Goldstein, P. Varghese, and C. Moore, “A velocity space hybridization-based boltzmann equation solver,” *Journal of Computational Physics* **408**, 1093–102 (2020).</footnote>

<footnote n="16">¹⁶V. I. Kolobov, “Fokker–planck modeling of electron kinetics in plasmas and semiconductors,” *Computational Materials Science* **28**, 302–320 (2003), proceedings of the Symposium on Software Development for Process and Materials Design.</footnote>

<footnote n="17">¹⁷S. Kuhn, M. Kamran, N. Jelić, L. Kos, J. Tskhakaya, D., and S. Tskhakaya, D. D., “Closure of the hierarchy of fluid equations by means of the polytropic-coefficient function (pcf),” *AIP Conference Proceedings* **1306**, 216–231 (2010).</footnote>

<footnote n="18">¹⁸T. Pichard, “Some recent advances on the method of moments in kinetic theory,” *ESAIM: ProcS* **75**, 86–95 (2023).</footnote>

<footnote n="19">¹⁹V. Zhdanov and G. Tirskii, “The use of the moment method to derive the gas and plasma transport equations with transport coefficients in higher-order approximations,” *Journal of Applied Mathematics and Mechanics* **67**, 365–388 (2003).</footnote>

<footnote n="20">²⁰P. Hunana, “Generalized fluid models of the braginskii type. part 2. the boltzmann operator,” (2025).</footnote>

<footnote n="21">²¹S. Boccelli, F. Giroux, T. E. Magin, C. P. T. Groth, and J. G. McDonald, “A 14-moment maximum-entropy description of electrons in crossed electric and magnetic fields,” *Physics of Plasmas* **27**, 123506 (2020).</footnote>

<footnote n="22">²²S. Boccelli, J. G. McDonald, and T. E. Magin, “14-moment maximum-entropy modeling of collisionless ions for Hall thruster discharges,” *Physics of Plasmas* **29**, 083903 (2022).</footnote>

<footnote n="23">²³P.-Y. C. Taunay and M. E. Mueller, “Quadrature-based moment methods for kinetic plasma simulations,” *Journal of Computational Physics* **473**, 111700 (2023).</footnote>
<footnote n="24">²⁴A. Alvarez Laguna, B. Esteves, A. Bourdon, and P. Chabert, “A regularized high-order moment model to capture non-Maxwellian electron energy distribution function effects in partially ionized plasmas,” *Physics of Plasmas* **29**, 083507 (2022).</footnote>

<footnote n="25">²⁵A. Alvarez Laguna, B. Esteves, J.-L. Raimbault, A. Bourdon, and P. Chabert, “Discussion on the transport processes in electrons with non-maxwellian energy distribution function in partially-ionized plasmas,” *Plasma Phys. and Controlled Fusion* **65**, 054002 (2023).</footnote>

<footnote n="26">²⁶A. Alvarez Laguna and T. Pichard, “Kinetic theory and moment models of electrons in a reactive weakly-ionized non-equilibrium plasma,” *Kinetic and Related Models* (2025), 10.3934/krm.2025007.</footnote>

<footnote n="27">²⁷D. A. Kuldinow, Y. Yamashita, A. R. Mansour, and K. Hara, “Ten-moment fluid model with heat flux closure for gasdynamic flows,” *Journal of Computational Physics* **508**, 113030 (2024).</footnote>

<footnote n="28">²⁸A. Alvarez Laguna, T. Pichard, T. Magin, P. Chabert, A. Bourdon, and M. Massot, “An asymptotic preserving well-balanced scheme for the isothermal fluid equations in low-temperature plasmas at low-pressure,” *Journal of Computational Physics* **419**, 109634 (2020).</footnote>

<footnote n="29">²⁹R. Sahu, A. R. Mansour, and K. Hara, “Full fluid moment model for low temperature magnetized plasmas,” *Phys. Plasmas* **27**, 113505 (2020).</footnote>

<footnote n="30">³⁰J. Laplante and C. P. T. Groth, “Comparison of maximum entropy and quadrature-based moment closures for shock transitions prediction in one-dimensional gaskinetic theory,” *AIP Conference Proceedings* **1786**, 140010 (2016).</footnote>

<footnote n="31">³¹Z. Cai, Y. Fan, and R. Li, “Globally hyperbolic regularization of grad’s moment system,” *Communications on Pure and Applied Mathematics* **67**, 464–518 (2014).</footnote>

<footnote n="32">³²J. McDonald and M. Torrilhon, “Affordable robust moment closures for cfd based on the maximum-entropy hierarchy,” *Journal of Computational Physics* **251**, 500–523 (2013).</footnote>

<footnote n="33">³³A. R. Baradaran, *Development and implementation of a preconditioner for a five-moment one-dimensional moment closure*, Ph.D. thesis, Université d’Ottawa/University of Ottawa (2015).</footnote>

<footnote n="34">³⁴S. Boccelli, W. Kaufmann, T. E. Magin, and J. G. McDonald, “Numerical simulation of rarefied supersonic flows using a fourth-order maximum-entropy moment method with interpolative closure,” *J. Comput. Phys.* **497** (2024), 10.1016/j.jcp.2023.112631.</footnote>
<footnote n="35"><sup>35</sup>R. O. Fox, F. Laurent, and A. Vié, “Conditional hyperbolic quadrature method of moments for kinetic equations,” Journal of Computational Physics **365**, 269–293 (2018).</footnote>

<footnote n="36"><sup>36</sup>C. Chalons, R. Fox, and M. Massot, “A multi-gaussian quadrature method of moments for gas-particle flows in a les framework,” in *Proceedings of the Summer Program* (2010) pp. 347–358.</footnote>

<footnote n="37"><sup>37</sup>G. M. Gangemi, A. Alvarez Laguna, M. Massot, K. Hillewaert, and T. Magin, “Bridging multifluid and drift-diffusion models for bounded plasmas,” *Physics of Plasmas* **32**, 023502 (2025).</footnote>

<footnote n="38"><sup>38</sup>“Phelps database, www.lxcat.net, retrieved on april 14, 2021.”.</footnote>

<footnote n="39"><sup>39</sup>K.-U. Riemann, “Plasma—sheath transition and bohm criterion,” *Contributions to Plasma Physics* **32**, 231–236 (1992).</footnote>

<footnote n="40"><sup>40</sup>K.-U. Riemann, J. Seebacher, D. D. Tskhakaya, and S. Kuhn, “The plasma-sheath matching problem,” *Plasma Physics and Controlled Fusion* **47**, 1949 (2005).</footnote>

<footnote n="41"><sup>41</sup>R. Robson, R. White, and M. Hildebrandt, *Fundamentals of Charged Particle Transport in Gases and Condensed Matter*, Monograph Series in Physical Sciences (CRC Press, 2017).</footnote>

<footnote n="42"><sup>42</sup>A. Berger, N. Lequette, T. Magin, A. Bourdon, and A. Alvarez Laguna, “Suitability of moment methods to model the ion dynamics in a low-temperature plasma,” (2025), submitted to *Proceedings of the 33rd International Symposium on Rarefied Gas Dynamics* (RGD 2024).</footnote>

<footnote n="43"><sup>43</sup>P. Hunana, T. Passot, E. Khomenko, D. Martínez-Gómez, M. Collados, A. Tenerani, G. P. Zank, Y. Maneva, M. L. Goldstein, and G. M. Webb, “Generalized fluid models of the braginskii type,” *The Astrophysical Journal Supplement Series* **260**, 26 (2022).</footnote>

<footnote n="44"><sup>44</sup>P. Lax, *Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves*, CBMS-NSF Regional Conference Series in Applied Mathematics No. n°s 11 à 16 (Society for Industrial and Applied Mathematics, 1973).</footnote>

<footnote n="45"><sup>45</sup>H. Grad, “On the kinetic theory of rarefied gases,” *Communications on Pure and Applied Mathematics* **2**, 331–407 (1949).</footnote>

<footnote n="46"><sup>46</sup>R. Balescu, “4 - the hermitian moment representation,” in *Classical Transport*, Transport Processes in Plasmas, edited by R. Balescu (North-Holland, Amsterdam, 1988) pp. 163–209.</footnote>

<footnote n="47"><sup>47</sup>W. Dreyer, “Maximisation of the entropy in non-equilibrium,” *Journal of Physics A: Mathematical and General* **20**, 6505 (1987).</footnote>
⁴⁸C. D. Levermore, “Moment closure hierarchies for kinetic theories,” Journal of Statistical Physics **83**, 1021–1065 (1996).

⁴⁹M. Pavić, T. Ruggeri, and S. Simić, “Maximum entropy principle for rarefied polyatomic gases,” Physica A: Statistical Mechanics and its Applications **392**, 1302–1317 (2013).

⁵⁰M. Junk, “Domain of definition of levermore’s five-moment system,” Journal of Statistical Physics **93**, 1143–1167 (1998).

⁵¹R. McGraw, “Description of aerosol dynamics by the quadrature method of moments,” Aerosol Science and Technology **27**, 255–265 (1997).

⁵²D. Marchisio and R. Fox, *Computational Models for Polydisperse Particulate and Multi-phase Systems*, Cambridge Series in Chemical Engineering (Cambridge University Press, 2013).

⁵³R. LeVeque, *Finite Volume Methods for Hyperbolic Problems*, Cambridge Texts in Applied Mathematics (Cambridge University Press, 2002).

⁵⁴E. Toro, *Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction* (Springer Berlin Heidelberg, 2009).

⁵⁵B. van Leer, “Towards the ultimate conservative difference scheme. v. a second-order sequel to godunov’s method,” Journal of Computational Physics **32**, 101–136 (1979).

⁵⁶G. Dal Maso, P. Le Floch, and F. Murat, “Definition and weak stability of nonconservative products,” Journal de mathématiques pures et appliquées **74**, 483–548 (1995).

⁵⁷C. Chalons and F. Coquel, “A new comment on the computation of non-conservative products using roe-type path conservative schemes,” Journal of Computational Physics **335**, 592–604 (2017).

⁵⁸M. Dumbser and D. S. Balsara, “A new efficient formulation of the hllem riemann solver for general conservative and non-conservative hyperbolic systems,” Journal of Computational Physics **304**, 275–319 (2016).

⁵⁹S. Boccelli, T. Charoy, A. Alvarez Laguna, P. Chabert, A. Bourdon, and T. E. Magin, “Collisionless ion modeling in Hall thrusters: Analytical axial velocity distribution function and heat flux closures,” Physics of Plasmas **27**, 073506 (2020).

⁶⁰G. M. Kremer, *An Introduction to the Boltzmann Equation and Transport Processes in Gases*, Interaction of Mechanics and Mathematics (Springer Berlin, Heidelberg, 2010).

⁶¹R. O. Fox, F. Laurent, and A. Passalacqua, “The generalized quadrature method of moments,” Journal of Aerosol Science **167**, 106096 (2023).
Comparison of high-order moment models for ions

<footnote n="62"><sup>62</sup>J. C. Wheeler, “Modified moments and gaussian quadratures,” The Rocky Mountain Journal of Mathematics **4**, 287–296 (1974).</footnote>