Stability Analysis of Quadrature-based Moment Methods for Kinetic Equations

TL;DR

This paper provides a theoretical stability analysis of quadrature-based moment methods for kinetic equations, showing conditions for hyperbolicity and preserving physical properties, which enhances understanding of their applicability.

Contribution

It offers the first thorough stability analysis of QBMM for Boltzmann equations, identifying conditions for hyperbolicity and the preservation of the H-theorem.

Findings

Gaussian approximation yields hyperbolic moment systems

Delta-function approximation can lead to non-hyperbolic systems

The equilibrium manifold is characterized on the boundary of the state space

Abstract

In this paper, we present a systematic stability analysis of the quadrature-based moment method (QBMM) for the one-dimensional Boltzmann equation with BGK or Shakhov models. As reported in recent literature, the method has revealed its potential for modeling non-equilibrium flows, while a thorough theoretical analysis is largely missing but desirable. We show that the method can yield non-hyperbolic moment systems if the distribution function is approximated by a linear combination of $\delta$-functions. On the other hand, if the $\delta$-functions are replaced by their Gaussian approximations with a common variance, we prove that the moment systems are strictly hyperbolic and preserve the dissipation property (or $H$-theorem) of the kinetic equation. In the proof we also determine the equilibrium manifold that lies on the boundary of the state space. The proofs are quite technical and involve detailed analyses of the characteristic polynomials of the coefficient matrices.

Equations

$$\frac{\partial f}{\partial t}+\xi \frac{\partial f}{\partial x}=Q(f).$$

$${\int }_{\mathbb{R}}Q(f)\varphi (\xi )d\xi =0,\varphi (\xi )=1,\xi ,{\xi }^2,$$

$$f_{eq}=f_{eq}(t,x,\xi )=\frac{\rho }{(2\pi \theta {)}^{1/2}}\exp \left(-\frac{{\left(\xi -U\right)}^2}{2\theta }\right),$$

$$\rho ={\int }_{\mathbb{R}}fd\xi ,\rho U={\int }_{\mathbb{R}}\xi fd\xi ,\rho \theta ={\int }_{\mathbb{R}}{\left(\xi -U\right)}^{2}fd\xi ,$$

$$Q=Q_{BGK}(f)=\nu (f_{eq}-f).$$

$$Q_S(f)=\nu (f_S-f).$$

$$f_S=f_{eq}\times \left(1+\frac{(1-Pr)q(\xi -U)}{3\rho {\theta }^2}\left(\frac{{\left(\xi -U\right)}^2}{\theta }-3\right)\right)$$

$$M_0=\rho ,M_1=\rho U,M_2=\rho (U^2+\theta ),M_3=\rho (U^3+3U\theta )+2q.$$

$${\partial }_t{M}_j+{\partial }_x{M}_{j+1}=\nu \left[\rho {\Delta }_j(U,\theta )-M_j\right].$$

$${\delta }_{\theta }(\xi ;U)=\frac{1}{\sqrt{2\pi \theta }}\exp \left(-\frac{(\xi -U{)}^2}{2\theta }\right).$$

$$f(\xi )=\sum _{i=1}^N{w}_i\delta (\xi -u_i).$$

$$M_j=\sum _{i=1}^N{w}_i{u}_i^j\text{for }j=0,\dots ,2N-1.$$

$${\bar{M}}_{2N}=\sum _{i=1}^N{w}_i{u}_i^{2N}.$$

$${\partial }_{t}M+A(M){\partial }_{x}M=\nu S(M).$$

$$A(M)=\left[\begin{array}{ccccc}0& 1& & & \\ & 0& 1& & \\ & & \ddots & \ddots & \\ & & & 0& 1\\ a_0& a_1& \cdots & a_{2N-2}& a_{2N-1}\end{array}\right],$$

$${\delta }_{{\sigma }^2}(\xi ;u)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{(\xi -u{)}^2}{2{\sigma }^2}\right),$$

$$f(\xi )=\sum _{i=1}^N{w}_i{\delta }_{{\sigma }^2}(\xi ;u_i).$$

$$M_j=\sum _{i=1}^N{w}_i{\Delta }_j(u_i,{\sigma }^2)\text{for }j=0,\dots ,2N,$$

$${\Omega }_W^{open}$$

$${\Omega }_W^{eq}$$

$${\bar{M}}_{2N+1}=\sum _{i=1}^N{w}_i{\Delta }_{2N+1}(u_i,{\sigma }^2).$$

$${\partial }_{t}M+A(M){\partial }_{x}M=\nu S(M)$$

$$A(M)=\left[\begin{array}{ccccc}0& 1& & & \\ & 0& 1& & \\ & & \ddots & \ddots & \\ & & & 0& 1\\ a_0& a_1& \cdots & a_{2N-1}& a_{2N}\end{array}\right]$$

$$\frac{\partial M}{\partial t}+\sum _{d=1}^D{A}_d(M)\frac{\partial M}{\partial x_d}=S(M).$$

$$P(M)S_M(M)=\left[\begin{array}{cc}0& 0\\ 0& \hat{T}(M)\end{array}\right]P(M),\forall M\in \mathcal{E};$$

$$A_0(M)A_d(M)=A_d^T(M)A_0(M)\text{for any }M\in \mathbb{G}\text{ and }d=1,\dots ,D;$$

$$A_0(M)S_M(M)+S_M^T(M)A_0(M)\le -P^T(M)\left[\begin{array}{cc}0& 0\\ 0& I_r\end{array}\right]P(M),\forall M\in \mathcal{E}.$$

$$v_k=\lambda v_{k-1}$$

$$\sum _{k=1}^{2N}{a}_{k-1}{v}_k$$

$$c(\lambda )={\lambda }^{2N}-a_{2N-1}{\lambda }^{2N-1}-\cdots -a_1\lambda -a_0.$$

$$(a_0,a_1,...,a_{2N-1})\frac{\partial M}{\partial W}=\left(\frac{\partial {\bar{M}}_{2N}}{\partial M}\right)\left(\frac{\partial M}{\partial W}\right)=\frac{\partial {\bar{M}}_{2N}}{\partial W}.$$

Full text

{AMS}

35Q79, 76P05, 82-08

1 Introduction

Kinetic theories pioneered by L. Boltzmann arise in a variety of fields beyond the classical rarefield gas dynamics, ranging from multiphase flows [14, 19], aerosol dynamics in atmospheric environments [10, 19], and active matter physics [13], to galactic dynamics in the universe [29]. In the kinetic framework [16], various physical systems are described with a distribution function $f$ which depends on the spatial and other problem-specific microscopic variables and its time evolution is governed by kinetic equations like the Boltzmann equation. Although the kinetic equations have solid physical ground, they are computationally costly and therefore not directly usable in engineering applications.

Because of the above reason, various simplifications or approximations of the kinetic equations have been proposed, including the BGK model [1], discrete velocity models [11, 21], and moment closure systems [12, 17, 19]. All these approximations have their advantages and disadvantages. This work is concerned with moment closure systems, in which the governing equations of several moments of the distribution function are derived from the kinetic equation and an additional procedure must be accompanied to close the moment system [19]. The resultant moment systems consist usually of first-order partial differential equations (PDEs).

To correctly model the observability of physical processes, the derived system of PDEs should be well-posed (or hyperbolic for first-order systems). For instance, the well-known Grad’s closure method yields non-hyperbolic PDEs and produces unphysical results [12, 22]. Its hyperbolic regularization has attracted much attention [2, 3, 17, 20, 26]. A recent work is [4] where the authors introduced a framework to construct hyperbolic moment closure systems.

Furthermore, the moment closure systems derived from the kinetic equations should preserve the key physical properties of the original kinetic equations. For the Boltzmann equation, one of the key properties is the celebrated $H$-theorem characterizing the dissipation property of the mesoscopic system under consideration [16]. In this regard, a paradigm is the widely used BGK model that not only simplifies the collision term in the Boltzmann equation, but also inherits the key conservation and dissipation properties thereof [16]. At this point, an immediate question is how to manifest the $H$-theorem in such moment systems.

It turns out that the structural stability condition proposed in [27] for hyperbolic relaxation systems is a proper counterpart of the $H$-theorem for the kinetic equation. Indeed, this condition has been tacitly respected by many well-developed physical theories [28]. Recently, it was shown in [7] to be satisfied by the hyperbolic regularization models derived in [2, 3, 4]. In contrast, the Biot/squirt (BISQ) model for wave propagation in saturated porous media violates this condition and thus allows exponentially exploding asymptotic solutions [18]. On the other hand, this condition also implies that the resultant moment system is compatible with the classical theories [27]. The implication is important because the lower-order moments are usually associated with the macroscopic parameters of the system [16]. Therefore, we believe that the structural stability condition is a proper criterion to evaluate the moment closure systems.

The objective of this paper is to investigate whether or not the quadrature-based moment method (QBMM, [19]) yields hyperbolic PDEs which satisfy the structural stability condition above. In QBMM, the distribution function $f$ is approximated with a linear combination of $N$ ($N\geq 1$) $\delta$-functions with unknown centers or their Gaussian approximations with unknown variance and centers (named QMOM or EQMOM, respectively) [19]. QBMM has become an effective and popular method in simulating the evolution of fine particulate matter, where the distribution function is independent of the particle velocity and the resultant governing equation is termed population balance equation [19, 23, 30]. However, the QMOM-derived moment system of the Boltzmann equation leads to unphysical shocks in the numerical solution of Riemann problems [9], which is confirmed by our own numerical results (see the Supplementary Material). Thus it is appealing to find the cause for the irregular behaviors and the aforementioned criteria are expected to be useful in clarifying such issues.

This paper deals only with the spatial one-dimensional (1-D) Boltzmann equation with hypothetical collisions (BGK or Shakhov type), just to figure out a road map for further investigations of general cases. We show that the QMOM-derived moment system is not strongly hyperbolic for any number $N$ of nodes, while the Gaussian EQMOM produces strictly hyperbolic moment systems when the variance is positive. For the latter, we further determine their equilibrium manifolds and verify the structural stability condition. The proofs are quite technical and purely analytic. They involve detailed analyses of characteristic polynomials of the coefficient matrices.

Let us remark that for $N=2$, the hyperbolicity of moment systems has been studied in [6] for 1-D QMOM and in [5] for 1-D Gaussian-EQMOM. The proofs rely on direct calculations of the eigenvalues of the coefficient matrix of the moment systems [6, 5] and does not seem generalizable to $N$-node systems. Thus new techniques are needed to handle the general cases. Moreover, the stability of EQMOM has not been analyzed in the existing literature. Given our positive results, EQMOM reveals its potential in solving a wider range of kinetic equations.

The paper is organized as follows. Section 2 presents a brief introduction on QBMM (QMOM and EQMOM) and states our main results. Section 3 is devoted to a proof of non-hyperbolicity of QMOM for $N$-node systems. In Section 4, we verify the structural stability condition for the EQMOM with $N$ nodes. In particular, the hyperbolicity is demonstrated in Section 4.2, the equilibrium states are determined in Section 4.3, and the dissipation property is shown in Sections 4.4 and 4.5. Finally, we conclude our paper in Section 5.

2 Preliminaries

For simplicity, we only consider a hypothetical 1-D ideal gas with the probability density function $f=f(t,x,\xi)$ of time $t\in\mathbb{R}^{+}$, spatial position $x\in\mathbb{R}$ and velocity $\xi\in\mathbb{R}$. The temporal evolution of $f$ is governed by the Boltzmann equation [16]:

[TABLE]

Here the volumetric force is neglected and the right-hand side $Q(f)$ represents the collisions. As a standard assumption [16], $Q=Q(f)$ has only 1, $\xi$ and $\xi^{2}$ as locally conserved quantities:

[TABLE]

and vanishes at a local equilibrium distribution

[TABLE]

where $\rho$, $U$ and $\theta$ are the density, velocity and temperature of the gas, respectively. They are the classical macroscopic parameters related to $f$ as

[TABLE]

In this paper, we mainly consider the BGK model [1], where

[TABLE]

Here $\nu$ is the collision frequency. This simple model has been widely used since it preserves several key properties of the kinetic equation, including Eq. 2 and the $H$-theorem. Because the BGK model results in the Prantle number $Pr=1$, inconsistent with most realistic cases [16], the Shakhov model was proposed [25]:

[TABLE]

Here an alternative equilibrium distribution $f_{S}$ is assumed:

[TABLE]

with $q$ the heat flux defined as $q=\int_{\mathbb{R}}\frac{1}{2}\left(\xi-U\right)^{3}fd\xi$.

Denote by $M_{j}(t,x)=\int_{\mathbb{R}}\xi^{j}fd\xi$ the $j$th velocity-moment of $f$. From Eq. 4 we see that

[TABLE]

The evolution equation for $M_{j}$ can be derived from the Boltzmann equation Eq. 1 with the BGK collision Eq. 5:

[TABLE]

Here $\Delta_{j}(U,\theta)$ denotes the $j$th moment of the normalized Gaussian distribution

[TABLE]

Notice that $\Delta_{0}(U,\theta)=1$ and $\Delta_{1}(U,\theta)=U$.

There are infinitely many equations in Eq. 9. The first $N$ equations for moments $M_{0},\dots,M_{N-1}$ are not closed, because the $M_{N-1}$-equation contains the term $\partial_{x}M_{N}$. Hence a closure method is needed.

In the rest of this section, we introduce the QBMM methods, the structural stability condition for hyperbolic relaxation systems, and our main results of this paper.

2.1 Quadrature-based moment methods

In QBMM, the lower-order moments determine the weights and nodes of the quadrature for the integration $\int f(\xi)g(\xi)d\xi$. Then the unclosed term can be expressed in terms of the lower-order moments and thereby the closure is done [19].

2.1.1 Quadrature method of moment (QMOM)

In QMOM, the distribution function $f$ is assumed to be a sum of $N$ Dirac delta functions

[TABLE]

In order to determine the weights $w_{i}$ and nodes $u_{i}$, the first $2N$ lower-order moments $M_{0},\dots,M_{2N-1}$ are employed:

[TABLE]

These non-linear algebraic equations can be solved to obtain $w_{i}$ and $u_{i}$ as in [19]. Then the next moment $M_{2N}$ can be found as

[TABLE]

Namely, $w_{i}$ and $u_{i}$ are functions of $M_{1},\dots,M_{2N-1}$, and so is $\bar{M}_{2N}$. In this way, we obtain the following system of PDEs:

[TABLE]

Here $M=(M_{0},\dots,M_{2N-1})^{T}\in\mathbb{R}^{2N}$, $S(M)=\rho(\Delta_{0}(U,\theta),\dots,\Delta_{2N-1}(U,\theta))^{T}-M$, and

[TABLE]

with $a_{j}=\frac{\partial\bar{M}_{2N}}{\partial M_{j}}$ for $0\leq j\leq 2N-1$.

2.1.2 Extended-QMOM (EQMOM)

In order to improve QMOM [5], the delta function in Eq. 10 is replaced with its Gaussian approximation

[TABLE]

that is,

[TABLE]

Set $W=(w_{1},u_{1},\dots,w_{N},u_{N},\sigma^{2})^{T}\in\mathbb{R}^{2N+1}$ and $M=(M_{0},\dots,M_{2N})^{T}\in\mathbb{R}^{2N+1}$. They are related with the map $M=\mathcal{M}(W)$:

[TABLE]

defined for $W\in\Omega_{W}=\Omega_{W}^{open}\cup\Omega_{W}^{eq}$, where

[TABLE]

Remark that $\Delta_{j}(u_{i},\sigma^{2})$ is exactly the same as that in Eq. 9. It is shown in Appendix A that the map $M=\mathcal{M}(W)$ is one-to-one for $W\in\Omega_{W}^{open}$. Therefore, $W$ can be uniquely solved from Eq. 16 for $M\in\mathcal{M}(\Omega_{W}^{open})$. In this way, the next moment $M_{2N+1}$ is a function of the lower-order moments $M\in\mathcal{M}(\Omega_{W}^{open})$:

[TABLE]

Therefore, the following moment system is derived:

[TABLE]

for $M\in\mathcal{M}(\Omega_{W}^{open})$. Here $S(M)=\rho\left(\Delta_{0}(U,\theta),\dots,\Delta_{2N}(U,\theta)\right)^{T}-M$ and

[TABLE]

with $a_{j}=\frac{\partial\bar{M}_{2N+1}}{\partial M_{j}}$ for $0\leq j\leq 2N$.

For such systems, the moment set $\mathcal{M}(\Omega_{W})$ and its closure $\overline{\mathcal{M}(\Omega_{W})}$ have been extensively studied as a realizability issue in the literature [5, 19, 23, 24]. A further discussion on this issue is beyond the scope of this paper.

2.2 Structural stability condition

The both QMOM and EQMOM moment systems consist of first-order PDEs derived from the Boltzmann equation. To clarify whether or not these systems inherits the $H$-theorem characterizing the dissipation property of the Boltzmann equation, we recall the structural stability condition proposed in [27] for systems of $D$-dimensional PDEs:

[TABLE]

Here $M$ is the unknown $n$-vector valued function, $A_{d}=A_{d}(M)$ is the $d$th $n\times n$ coefficient matrix, and the source term $S=S(M)$ is a given $n$-vector valued function of $M\in\mathbb{G}\subset\mathbb{R}^{n}$. As in [27], we assume that the equilibrium manifold $\mathcal{E}=\{M\in\mathbb{G}\ |\ S(M)=0\}$ is not empty and denote the Jacobian matrix of $S(M)$ as $S_{M}(M)$. The stability condition reads as

• (i)

There exist an invertible $n\times n$ matrix $P(M)$ and an invertible $r\times r$ ($0<r\leq n$) matrix $\hat{T}(M)$ such that

[TABLE]

• (ii)

There exists a positive definite symmetric matrix $A_{0}(M)$ such that

[TABLE]

• (iii)

The spatial derivative parts and the source are coupled as

[TABLE]

Here $I_{r}$ is the unit matrix of order $r$.

As shown in [28], this set of conditions has been tacitly respected by many well-developed physical theories. Condition (i) is classical for initial value problems of the system of ordinary differential equations (ODE, spatially homogeneous systems), while (ii) means the symmetrizable hyperbolicity of the PDE system. Condition (iii) characterizes a kind of coupling between the ODE and PDE parts. Recently, this structural stability condition is shown in [7] to be proper for certain moment closure systems. On the other hand, this set of conditions implies the existence and stability of the zero relaxation limit of the corresponding initial value problems [27]. Thanks to these, we believe that the structural stability condition is essential for a reasonable moment closure system.

2.3 Main results

For the moment systems derived above, we will establish the following facts as the main result of this paper,

Theorem 2.1** (Non-hyperbolicity of QMOM).**

*The QMOM-derived moment system Eq. 13 is not strongly hyperbolic. *

Theorem 2.2** (Stability of EQMOM).**

*The EQMOM-derived moment system Eq. 19 satisfies the structural stability condition for $M\in\mathcal{M}(\Omega_{W})$. *

A proof of Theorem 2.1 will be presented in the next section. In Section 4, Theorem 2.2 is divided as Theorem 4.6 (hyperbolicity of EQMOM), Theorem 4.17 (equilibrium state), Theorem 4.19 (BGK model) and Theorem 4.22 (Shakhov model), which will be proved in Sections 4.2, 4.3, 4.4, and 4.5, respectively.

3 Non-hyperbolicity of QMOM

This section is devoted to a proof of Theorem 2.1 for the QMOM-derived moment system Eq. 13 with $N\geq 2$. We should mention that this theorem has been proved in [6] but only for $N=2$. For our purpose, we need to consider the $2N\times 2N$ coefficient matrix $A=A(M)$ in Eq. 14.

Proof 3.1** (Proof of Theorem 2.1).**

Let $\lambda$ be an eigenvalue of $A$ and $\mathbf{v}=(v_{1},\dots,v_{2N})^{T}$ the corresponding right eigenvector. A direct calculation indicates that

[TABLE]

Then we have $\mathbf{v}=v_{1}(1,\lambda,\dots,\lambda^{2N-1})^{T}$ and thereby $v_{1}\neq 0$. This shows that the geometric multiplicity of each eigenvalue is 1.

On the other hand, we see from Eq. 22b that the characteristic polynomial of $A$ is

[TABLE]

Note that $(a_{0},a_{1},\dots,a_{2N-1})=\left(\frac{\partial\bar{M}_{2N}}{\partial M_{0}},\frac{\partial\bar{M}_{2N}}{\partial M_{1}},\dots,\frac{\partial\bar{M}_{2N}}{\partial M_{2N-1}}\right)=\frac{\partial\bar{M}_{2N}}{\partial M}$ with $\bar{M}_{2N}$ defined in Eq. 12 and $M=(M_{0},\dots,M_{2N-1})^{T}$. Writing $W=(w_{1},u_{1},\dots,w_{N},u_{N})^{T}\in\mathbb{R}^{2N}$, we have

[TABLE]

In addition, it follows from Eq. 11 that the $2N\times 2N$ Jacobian matrix $\partial M/\partial W$ is

[TABLE]

and from Eq. 12 that

[TABLE]

Substituting the last two relations into Eq. 24, we obtain

[TABLE]

for $k=1,\dots,N$. These mean that $c(u_{k})=0$ and $\left.\frac{dc(\lambda)}{d\lambda}\right|_{\lambda=u_{k}}=0$ for $k=1,\dots,N$. Since $c=c(\lambda)$ is a monic polynomial of order $2N$, there must be

[TABLE]

As a result, the eigenvalues of $A$ are $u_{1},u_{2},\dots,u_{N}$ and each of them has the algebraic multiplicity 2 and the geometric multiplicity 1. In view of its Jordan canonical form, the coefficient matrix $A$ is similar to

[TABLE]

*Hence the moment closure system Eq. 13 is not strongly hyperbolic. *

4 Stability of EQMOM

We prove Theorem 2.2 in this section. In particular, Section 4.2 is devoted to Condition (ii), while Conditions (i) and (iii) are verified in Sections 4.4 and 4.5 for both the BGK and Shakhov collision models.

4.1 Preliminaries

Recall that in Section 2.1, we use the notation

[TABLE]

for the $j$th moment of the Gaussian distribution $\delta_{\sigma^{2}}=\delta_{\sigma^{2}}(\xi;u)=\frac{1}{\sqrt{2\pi}\sigma}\exp{\left(-\frac{(\xi-u)^{2}}{2\sigma^{2}}\right)}$. A direct calculation shows $\Delta_{0}(u,\sigma^{2})=1$ and $\Delta_{1}(u,\sigma^{2})=u$. Moreover, we can show with Lemma 4.1(a) below that $\Delta_{j}(u,\sigma^{2})$ is a bivariate polynomial of $u$ and $\sigma^{2}$.

Lemma 4.1**.**

[TABLE]

Proof 4.2**.**

(a): Note that $d\delta_{\sigma^{2}}/d\xi=-(\xi-u)\delta_{\sigma^{2}}/{\sigma^{2}}$. Then for $j\geq 2$ we have

[TABLE]

This, together with $\Delta_{0}(u,\sigma^{2})=1$ and $\Delta_{1}(u,\sigma^{2})=u$, indicates that $\Delta_{j}=\Delta_{j}(u,\sigma^{2})$ is a polynomial of both $u$ and $\sigma^{2}$.

(b): This can be proven by induction on $j$. It obviously holds for $\Delta_{0}=1$ and $\Delta_{1}=u$. Suppose it is true for $j-1$ and $j$. Then for $j+1$ it follows from (a) that

[TABLE]

Hence the proof is complete.

(c & d): These two follow immediately from (b):

[TABLE]

Remark 4.3**.**

Lemma 4.1*(b) is obviously equivalent to*

[TABLE]

*which was established in [19]. But the former is more convenient for our later use. *

Inspired by Lemma 4.1(b), we introduce a family of linear operators $\mathcal{D}_{\vartheta}$, parameterized with $\vartheta\in\mathbb{R}$, acting on the polynomial algebra $\mathbb{R}[u]$. For $f\in\mathbb{R}[u]$, $\mathcal{D}_{\vartheta}f$ is defined as

[TABLE]

which is a finite sum. Obviously, $\mathcal{D}_{0}$ is an identical operator, $\mathcal{D}_{\sigma^{2}}f$ is a polynomial of $u$ and $\sigma^{2}$, and $\mathcal{D}_{\sigma^{2}}u^{j}=\Delta_{j}(u,\sigma^{2})$. Further useful properties of $\mathcal{D}_{\vartheta}$ are

Lemma 4.4**.**

[TABLE]

Proof 4.5**.**

(a): For the composition, we deduce from the definition that

[TABLE]

(b) follows immediately from (a) and $\mathcal{D}_{0}=id$.

For (c), the first one is obvious, while the second can be shown as Lemma 4.1(d).

*(d): By using $(uf)^{(2k)}=uf^{(2k)}+2kf^{(2k-1)}$, this can be proved as Lemma 4.1(b). Then (e) follows immediately from (d). *

4.2 Hyperbolicity of EQMOM

In this section we prove that the EQMOM-derived moment system Eq. 19 for the 1-D Boltzmann equation is strictly hyperbolic, which will be shown to be sufficient for the structural stability condition (ii). The conclusion can be stated as

Theorem 4.6**.**

*For $M\in\mathcal{M}(\Omega_{W})$, the $(2N+1)\times(2N+1)$ coefficient matrix $A=A(M)$ in Eq. 20 has $(2N+1)$ distinct real eigenvalues. Namely, the EQMOM-derived moment system Eq. 19 is strictly hyperbolic. *

We should mention that this theorem was already established in [5] for $N=2$ (the two-node system) but the proof does not seem to work for $N>2$.

Our proof of this theorem needs some preparations. First of all, the characteristic polynomial of $A$ in Eq. 20 reads as

[TABLE]

Here the coefficient $a_{j}=a_{j}(W)=\frac{\partial\bar{M}_{2N+1}}{\partial M_{j}}$ ($j=0,1,\dots,2N$), with $\bar{M}_{2N+1}$ defined in Eq. 18, is a function of $W$. To show that $c(u;W)$, as a polynomial of $u$, has $(2N+1)$ distinct real roots for $M\in\mathcal{M}(\Omega_{W})$, we introduce an auxiliary function

[TABLE]

By Lemma 4.4(b), we have

[TABLE]

Set $a_{2N+1}=-1$. Then $c(u;W)$ can be rewritten as $-\sum_{j=0}^{2N+1}a_{j}u^{j}$ and from the linearity of $\mathcal{D}_{\sigma^{2}}$ it follows that

[TABLE]

Moreover, from Eqs. 29 and 30 we see that $g(u;W)$ is a $u$-polynomial of degree $(2N+1)$:

[TABLE]

with $g_{2N+1}=a_{2N+1}=-1$. Further relations between the coefficients of $g(u;W)$ and $c(u;W)$ are

[TABLE]

This can be shown as

[TABLE]

Furthermore, $g=g(u)=g(u;W)$ has the following elegant expression.

Lemma 4.7**.**

[TABLE]

where $u_{1},\dots,u_{N}$ are the nodes solved from Eq. 16, and

[TABLE]

*for $W\in\Omega_{W}^{open}$. *

Remark 4.8**.**

*This lemma shows that for $W\in\Omega_{W}^{open}$, $\tilde{U}$ is a convex combination of the $u_{i}$’s. Moreover, for $W=(w_{1},U,w_{2},U,\dots,w_{N},U,\sigma^{2})\in\Omega_{W}^{eq}$ and any sequence $\{W_{k}\}\subset\Omega_{W}^{open}$ approaching $W$, $\tilde{U}(W_{k})$ converges to $U$. Because of this, for $W\in\Omega_{W}^{eq}$ we define $\tilde{U}(W)=U$ ($=M_{1}/M_{0}$) and thereby $g(u;W)=(u-U)^{2N+1}$. *

Proof 4.9**.**

By Lemma 4.1(c&d), the Jacobian matrix of the map $M=\mathcal{M}(W)$ defined in Eq. 16 is

[TABLE]

Note that the dependence of $\Delta_{j}$ on $\sigma^{2}$ has been omitted here for clarity. Moreover, from Eq. 18, $\partial\bar{M}_{2N+1}/\partial W$ reads as

[TABLE]

Then from the simple relation

[TABLE]

we obtain

[TABLE]

Eqs. 36a* and 36b, together with Eq. 32 and Lemma 4.1(c), imply that $g(u_{k})=0$ and $\left.\frac{dg(u)}{du}\right|_{u=u_{k}}=0$ for $k=1,\dots,N$. Thus we see the expected expression of $g(u;W)$ from Eq. 33 with $\tilde{U}$ to be determined.*

Next, we use Eq. 36c to determine $\tilde{U}$. Recall that $a_{2N+1}=-1$. We use Eq. 16 to rewrite Eq. 36c as

[TABLE]

On the other hand, we deduce from Eqs. 32 and 33 that

[TABLE]

Then we see from Eq. 37 that

[TABLE]

Now we define $\tilde{g}(u)=(u-u_{1})^{2}\cdots(u-u_{N})^{2}=-\sum_{j=0}^{2N}\tilde{g}_{j}u^{j}$ with $\tilde{g}_{2N}=-1$. Then $g(u)=(u-\tilde{U})\tilde{g}(u)$ and the coefficients are related with

[TABLE]

for $0\leq j\leq 2N+1$ ($\tilde{g}_{-1}=\tilde{g}_{2N+1}=0$). Substituting this relation into Eq. 38, we obtain

[TABLE]

where $M_{j}^{*}=\sum_{i=1}^{N}w_{i}u_{i}^{j}$. It remains to show

[TABLE]

These two follow from the obvious relations

[TABLE]

*for any $1\leq i\leq N$. This completes the proof. *

Remark 4.10**.**

Lemma 4.7* indicates that the coefficients $g_{j}$ of $g(u;W)$ in Eq. 33 are independent of $\sigma^{2}$. From Eq. 31 we see that $c(u;W)$ is a bivariate polynomial of $u$ and $\sigma^{2}$, the coefficients $a_{j}$ of $c(u;W)$ are polynomials of $\sigma^{2}$, and $c(u;W)=g(u)$ for $\sigma^{2}=0$. Furthermore, the $j$th derivative $c^{(j)}(u;W)$ of $c(u;W)$ with respect to $u$ can be viewed as a perturbation of $g^{(j)}(u)$ with the single parameter $\sigma^{2}\geq 0$ for $0\leq j\leq 2N+1$. *

By Lemma 4.7, $g(u)$ has $(2N+1)$ ($=$ the degree of $g$) real roots (including multiplicity). This fact can be further generalized as follows.

Lemma 4.11**.**

*For any $0\leq j\leq 2N$, $g^{(j)}(u)$ has $(2N+1-j)$ real roots (including multiplicity). Hence any local minimum (maximum) value of $g^{(j)}(u)$ is non-positive (non-negative). *

Proof 4.12**.**

We prove by induction on $j$. As discussed above, the conclusion holds for $j=0$. Namely, $g$ has $(2N+1)$ roots. Suppose it holds for $0,\dots,j$. Then we have $g^{(j)}(u)=C(u-\tilde{u}_{1})^{k_{1}}\cdots(u-\tilde{u}_{m})^{k_{m}}$, where $m\geq 1$, $\tilde{u}_{1}<\tilde{u}_{2}<\dots<\tilde{u}_{m}$, $k_{i}\geq 1$ and $k_{1}+\dots+k_{m}=2N+1-j$. Thus $(u-\tilde{u}_{i})^{k_{i}-1}$ is a factor of $g^{(j+1)}(u)$ for any $1\leq i\leq m$. Besides, Rolle’s theorem implies the existence of at least one root of $g^{(j+1)}(u)$ in each open interval $(\tilde{u}_{i},\tilde{u}_{i+1})$ for $1\leq i\leq m-1$. Therefore, the number of roots of $g^{(j+1)}$ is no less than

[TABLE]

*Since $g^{(j+1)}$ is of degree $(2N-j)$, it must have $(2N-j)$ roots (including multiplicity). This also indicates that $g^{(j)}$ has only one extreme point in each open interval above. Hence any local minimum (maximum) value of $g^{(j)}(u)$ is non-positive (non-negative). *

With the preparations above, we are in a position to prove Theorem 4.6.

Proof 4.13** (Proof of Theorem 4.6).**

We will prove the following stronger statement: for $0\leq j\leq 2N+1$, $c^{(2N+1-j)}(u;W)$ has $j$ distinct roots for any $W=(w_{1},u_{1},\dots,w_{N},u_{N},\sigma^{2})$ $\in\Omega_{W}$ with $\sigma^{2}>0$. This will be done with induction on $j$. For $j=0,\ 1$, the statement is obvious because $c^{(2N+1)}(u;W)=(2N+1)!$ and $c^{(2N)}(u;W)$ is of degree 1.

Suppose the conclusion holds for $j\leq k(\leq 2N+1)$. From Remark 4.10 we know that $c^{(2N+1-k)}(u;W)$ is a bivariate polynomial of $u$ and $\sigma^{2}$ on $\mathbb{R}\times[0,\infty)$. Denote $u^{*}(\sigma^{2})\in\mathbb{R}$ to be one root of $c^{(2N+1-k)}(u;W)$. Thus $u^{*}(\sigma^{2})$ is an extreme point of $c^{(2N-k)}(u;W)$ and $u^{*}(0)$ is a root of $g^{(2N+1-k)}(u)$. Moreover, $u^{*}(\sigma^{2})$ is continuous on $\sigma^{2}\in[0,\infty)$ and differentiable on $(0,\infty)$ because the roots are distinct [15].

Next we consider the extreme values of $c^{(2N-k)}(u;W)$ at $u=u^{*}(\sigma^{2})$. Since $c^{(2N-k)}(u;W)$ is a polynomial of $u$ and $\sigma^{2}$, the composite $h_{k}(\sigma^{2}):=c^{(2N-k)}(u^{*}(\sigma^{2});W)$ is continuous on $[0,\infty)$ and differentiable on $(0,\infty)$. And $h_{k}(0)$ is the extreme value of $g^{(2N-k)}(u)$. According to Lemma 4.11, $h_{k}(0)\geq 0$ if it is a local maximum and $h_{k}(0)\leq 0$ if it is a local minimum. For $\sigma^{2}>0$, because $c^{(2N+1-k)}(u^{*}(\sigma^{2});W)=0$, the derivative of $h_{k}(\sigma^{2})$ reads as

[TABLE]

Thus, if $c^{(2N-k+2)}(u^{*}(\sigma^{2});W)<0$ (that is, $u^{*}(\sigma^{2})$ is a local maximum point of $c^{(2N-k)}(u;W)$), then the local maximum value $h_{k}(\sigma^{2})$ strictly increases on $\sigma^{2}\in(0,\infty)$. Since $h_{k}(0)\geq 0$ and $h_{k}(\sigma^{2})$ is continuous at $\sigma^{2}=0$, we conclude that $h_{k}(\sigma^{2})>0$ for all $\sigma^{2}>0$. Similarly, if $c^{(2N-k+2)}(u^{*}(\sigma^{2});W)>0$, we have $h_{k}(\sigma^{2})<0$ for all $\sigma^{2}>0$.

*In summary, the above arguments show that each local maximum value of the ($k+1$)th oder polynomial $c^{(2N-k)}(u;W)$ is positive and each local minimum value is negative. On the other hand, by the induction assumption $c^{(2N+1-k)}(u;W)$ has $k$ distinct real roots, which are naturally extreme points of $c^{(2N-k)}(u;W)$ for $\sigma^{2}>0$. Therefore, $c^{(2N-k)}(u;W)$ has $(k-1)$ distinct real roots among the extreme points. Moreover, the induction assumption implies that $c^{(2N-k)}(u;W)$ has one root larger and another one less than all the extreme points. Thus, for each $\sigma^{2}>0$, $c^{(2N-k)}(u;W)$ has $(k+1)$ distinct real roots. By the induction principle this completes the proof. *

Remark 4.14**.**

*By Theorem 4.6, the coefficient matrix $A=A(M)$ of the 1-D moment system Eq. 19 has $n=(2N+1)$ distinct real eigenvalues $\lambda_{i}$ ($1\leq i\leq n$) for $\sigma^{2}>0$. Denote by $r_{i}$ the corresponding left eigenvectors. Set $L=(r_{1}^{T},\dots,r_{n}^{T})^{T}$. It is clear that $A_{0}(M)=L^{T}\Lambda L$ with $\Lambda$ an arbitrary positive diagonal matrix is a symmetrizer in the structural stability condition (ii). As a matter of fact, it is straightforward to show that such a symmetrizer can only be of the form $L^{T}\Lambda L$ . *

4.3 Equilibrium state

As stated in Section 2.2, (i) and (iii) of the structural stability condition should be examined on the equilibrium manifold $\mathcal{E}$ where $S(M(W))=0$. In this section we determine the equilibrium manifold.

For the BGK model, $S(M(W))=0$ is equivalent to

[TABLE]

(see Section 2.1.2). Thus, the equilibrium state $W=(w_{1},u_{1},\dots,w_{N},u_{N},\sigma^{2})^{T}$ is determined by the three macroscopic parameters $\rho$, $U$ and $\theta$. And we need to find $W$ from Eq. 39 for $1\leq i\leq N$.

For this purpose, we recall from Section 4.1 that $\Delta_{0}(u,\sigma^{2})=1$, $\Delta_{1}(u,\sigma^{2})=u$ and $\Delta_{2}(u,\sigma^{2})=u^{2}+\sigma^{2}$. Thus, for $j=0,1,2$, Eq. 39 is just

[TABLE]

Then we deduce from the inequality $\left(\sum_{i=1}^{N}w_{i}\right)\left(\sum_{i=1}^{N}w_{i}u_{i}^{2}\right)\geq\left(\sum_{i=1}^{N}w_{i}u_{i}\right)^{2}$ that

[TABLE]

For further discussions, we need the following fact.

Proposition 4.15**.**

[TABLE]

Proof 4.16**.**

Recall that $\Delta_{j}(u,\sigma^{2})=\mathcal{D}_{\sigma^{2}}u^{j}$. From Lemma 4.4(b) and Lemma 4.1(c) we deduce that

[TABLE]

*Then taking the weighted summation $\sum_{i=1}^{N}w_{i}$ and using Eq. 16 give the proposition. *

Next we define $\zeta^{2}=\theta-\sigma^{2}\geq 0$ and show that Eq. 39 is equivalent to

[TABLE]

Indeed, if Eq. 39 holds (i.e. $M(W)\in\mathcal{E}$), the last proposition implies that

[TABLE]

Here the expression $\mathcal{D}_{\vartheta}f(U)$ denotes $\mathcal{D}_{\vartheta}f(u)|_{u=U}$ for arbitrary polynomial $f$ and the last step is due to Lemma 4.4(a). This is just Eq. 41. The deduction of Eq. 39 from Eq. 41 is similar.

Now we are in a position to state the central result of this section.

Theorem 4.17**.**

The equilibrium state belongs to $\Omega_{W}^{eq}$, that is,

[TABLE]

*Hence, at equilibrium $\bar{M}_{2N+1}=\rho\Delta_{2N+1}(U,\theta)$. *

Proof 4.18**.**

Thanks to Eq. 40, it suffices to show that $\zeta^{2}:=\theta-\sigma^{2}=0$. Otherwise, the $u_{i}$’s must take $N^{\prime}$ different values ($1<N^{\prime}\leq N$). Then, by redefining $w_{i}$, the summation $\sum_{i=1}^{N}w_{i}u_{i}^{j}$ in the left-hand side of Eq. 41 is reduced to $\sum_{k=1}^{N^{\prime}}w_{k}u_{k}^{j}$ where the $u_{k}$’s are distinct ($1\leq k\leq N^{\prime}$). Thus, we may as well assume that all the $u_{i}$’s are distinct and $\zeta^{2}>0$. Then we will derive a contradiction in three steps, where the abbreviation

[TABLE]

will be frequently used.

Step I. Because $\Delta_{j}(u,\zeta^{2})=\mathcal{D}_{\zeta^{2}}u^{j}$, the first $N$ equations ($j=0,\dots,N-1$) in Eq. 41 can be rewritten as a system of linear algebraic equations:

[TABLE]

Since all the $u_{i}$’s are distinct, this gives a unique $(w_{1},\dots,w_{N})$ in terms of $(u_{1},\dots,u_{N})$, $\rho$, $U$ and $\zeta^{2}$. We claim that for $1\leq i\leq N$,

[TABLE]

To see this, we use the uniqueness and only need to show that the $w_{i}$’s solve the system of equations above. Indeed, thanks to the Lagrange interpolating polynomial

[TABLE]

and the linearity of the operator $\mathcal{D}_{\zeta^{2}},$Eq. 42 implies that for $0\leq j\leq N-1$,

[TABLE]

Namely, the $w_{i}$’s defined in Eq. 42 solve the system of linear algebraic equations above.

Step II. With the $w_{i}$’s defined in Eq. 42, we turn to the next $N$ equations ($j=N,\dots,2N-1$) in Eq. 41 to solve $u_{i}$:

[TABLE]

Again, the solution $w_{i}u_{i}^{N}$ is unique. As in Step I, we can show that

[TABLE]

for $1\leq i\leq N$.

Substituting Eq. 42 into Eq. 43, we obtain

[TABLE]

for $1\leq i\leq N$. By the linearity of $\mathcal{D}_{\zeta^{2}}$, Eq. 44 is equivalent to

[TABLE]

which can be rewritten as

[TABLE]

with $F=F(U)=\prod_{k=1}^{N}(U-u_{k})$. Since all the $u_{i}$’s are distinct, this says

[TABLE]

Having this, in Lemma 4.4(e) we take $u_{0}=U$ and $f(u)=u^{j}F(u)$ ($0\leq j\leq N-2$) and deduce from Eq. 45 that

[TABLE]

Hence $\mathcal{D}_{\zeta^{2}}(U^{j}F^{\prime})=0$ for $0\leq j\leq N-2$. This procedure can be repeated for the derivative of $F^{\prime}(u)$ to yield $\mathcal{D}_{\zeta^{2}}(U^{j}F^{\prime\prime})=0$ for $0\leq j\leq N-3$. Moreover, we have

[TABLE]

for $0\leq k\leq N-1$ and $0\leq j\leq N-1-k$.

Step III. In this step, we use Eq. 43 and the Lagrange interpolating polynomial

[TABLE]

to deduce that

[TABLE]

Thus, the last equation in Eq. 41 is equivalent to $\mathcal{D}_{\zeta^{2}}(U^{N}(U^{N}-F))=\mathcal{D}_{\zeta^{2}}(U^{2N})$ or

[TABLE]

Then we use Lemma 4.4(d) and Eq. 46 to see that

[TABLE]

*which implies that $\zeta^{2}=0$. This contradicts the assumption that all the $u_{i}$’s are distinct and $\zeta^{2}>0$. Hence the proof is complete. *

4.4 BGK model

In this subsection we show that the EQMOM moment system Eq. 19 with BGK source term

[TABLE]

satisfies the structural stability condition (i)–(iii). Indeed, (ii) has been verified in Remark 4.14.

To see Condition (i), we compute the Jacobian matrix of $S=S(M)$. Notice that the first three components of $S$ vanish identically and $\rho,U,\theta$ depend only on $M_{0}$, $M_{1}$ and $M_{2}$. Then the Jacobian matrix can be written as

[TABLE]

where $\hat{S}_{M}$ is a $(2N-2)\times 3$ matrix with

[TABLE]

for $3\leq i\leq 2N$ and $j=0,1,2$. Now we take

[TABLE]

and see that $PS_{M}=\begin{bmatrix}0_{3\times 3}&\\ &-I_{2N-2}\end{bmatrix}P$, which justifies Condition (i).

The rest of this subsection is to show Condition (iii). To this end, we need to choose the symmetrizer $A_{0}=A_{0}(M)$. As pointed out in Remark 4.14, such a symmetrizer $A_{0}$ can only be of the form $L^{T}\Lambda L$ with $\Lambda$ a diagonal positive definite matrix to be determined.

Firstly, we specify the matrix $L=(r_{1}^{T},\dots,r_{2N+1}^{T})^{T}$ with $r_{i}$ a left eigenvector of the coefficient matrix $A=A(M)$ corresponding to the eigenvalues $\lambda_{i}$ for $1\leq i\leq 2N+1$. Let $r_{i}=\left(r_{i}^{(1)},\dots,r_{i}^{(2N+1)}\right)$. From $r_{i}A=\lambda_{i}r_{i}$ we have

[TABLE]

Here we have assumed $r_{i}^{(0)}=0$ for simplicity. From the last equation we see that $r_{i}^{(2N+1)}\neq 0$; otherwise the eigenvector $r_{i}=0$. Thus we may as well assume $r_{i}^{(2N+1)}=1$. Recall that $a_{2N+1}=-1$. Then we can easily obtain

[TABLE]

for $0\leq j\leq 2N$. Therefore, we have

[TABLE]

With this $L$, we can state our main result of this subsection.

Theorem 4.19**.**

*For the EQMOM moment system Eq. 19, the inequality in the structural stability condition (iii) holds with $A_{0}=L^{T}L$ and $P$ defined in Eq. 49. *

Proof 4.20**.**

According to Theorem 2.1 in [27], it suffices to show that at equilibrium states $M$,

[TABLE]

is of the block-diagonal form $\operatorname{diag}(K_{1},K_{2})$, in which $K_{1}$ and $K_{2}$ are $3\times 3$ and $(2N-2)\times(2N-2)$ matrices, respectively. Namely, the first three columns of $LP^{-1}$ are orthogonal to its other columns. In what follows all the states $M$ are in equilibrium.

To show the orthogonality, we compute the $(2N+1)$-matrix $LP^{-1}:=(b_{il})$. From Eq. 49 we see that

[TABLE]

This, together with Eq. 50, gives

[TABLE]

The expression above indicates that the last $(2N-2)$ columns of $LP^{-1}$ are linear combinations of $(\lambda_{1}^{\beta},\dots,\lambda_{2N+1}^{\beta})^{T}\in\mathbb{R}^{2N+1}$ for $0\leq\beta\leq 2N-3$. Thus it reduces to show that

[TABLE]

Set $p_{k}=\sum_{i=1}^{2N+1}\lambda_{i}^{k}$. By using the above expression of $b_{il}$ for $1\leq l\leq 3$, the last equation is equivalent to

[TABLE]

for $l=0,1,2$ and $0\leq\beta\leq 2N-3$.

To prove Eq. 52, it suffices to show that

[TABLE]

where $\mathcal{H}_{j}$ can be replaced by any of $\Delta_{j}=\Delta_{j}(U,\theta)$, $\partial_{U}\Delta_{j}$ and $\partial_{\theta}\Delta_{j}$. Indeed, Eq. 52 follows immediately from Eq. 53 and Eq. 48 which says

[TABLE]

for $0\leq j\leq 2N$ and $l=0,1,2$.

Before proceeding, two tools are needed. The first one is Newton’s power sum formulas for $p_{k}$ [8]:

[TABLE]

The second tool is the following relation

[TABLE]

for $0\leq k\leq 2N$ and $0\leq j\leq 2N-k$, where $c^{(j)}$ denotes the $j$th derivative of the characteristic polynomial $c=c(u;W)$ with respect to $u$. This relation can be proved as below. Lemma 4.7 tells $g(u)=(u-U)^{2N+1}$ in equilibrium and therefore $g^{(j)}(U)=0$ for $0\leq j\leq 2N$. From Lemma 4.4(c) we see that $g^{(j)}(u)=\mathcal{D}_{\theta}c^{(j)}$ and thereby $\left.\mathcal{D}_{\theta}c^{(j)}\right|_{u=U}=0$ for $0\leq j\leq 2N$. This is just the case for $k=0$ in Eq. 55. Then using Lemma 4.4(d) we have

[TABLE]

for $j=0,\dots,2N-1$, which validates the case for $k=1$. This procedure can be repeated to show Eq. 55 for other $k\leq 2N$.

With these preparations, we only need to prove Eq. 53 for the following two cases.

Case I: $\beta=0$. Noting that $p_{0}=2N+1$, we deduce from Eq. 54b that

[TABLE]

Thus Eq. 53 in this case is equivalent to

[TABLE]

When taking $\mathcal{H}_{j}$ to be $\Delta_{j}$, $\partial_{U}\Delta_{j}$ or $\partial_{\theta}\Delta_{j}$, the left-hand side of the last equation is equivalent to $\left.\mathcal{D}_{\theta}c^{\prime}\right|_{u=U}$, $\left.\mathcal{D}_{\theta}c^{\prime\prime}\right|_{u=U}$ or $\left.\mathcal{D}_{\theta}c^{\prime\prime\prime}\right|_{u=U}$, respectively. They are all equal to zero due to Eq. 55 and hence Eq. 53 with $\beta=0$ is proved.

Case II: $\beta\geq 1$. As in Case I, we first simplify the coefficients $\sum_{k=\beta}^{2N-j}a_{j+k+1}p_{k}$ of $\mathcal{H}_{j+\beta}$ in Eq. 53 by using Newton’s power sum formulas Eqs. 54b and 54a. They can be rewritten as

[TABLE]

With these two relations, Eq. 53 is equivalent to

[TABLE]

for $1\leq\beta\leq 2N-3$.

Eq. 56* can be further simplified by using the following relations*

[TABLE]

Indeed, replacing $\mathcal{H}_{j}$ by $\Delta_{j}$, $\partial_{U}\Delta_{j}$ or $\partial_{\theta}\Delta_{j}$, the sum is just $\left.\mathcal{D}_{\theta}\left(u^{\beta-1-k}c\right)\right|_{u=U}$, $\left.\mathcal{D}_{\theta}\left(u^{\beta-1-k}c^{\prime}\right)\right|_{u=U}$ or $\left.\mathcal{D}_{\theta}\left(u^{\beta-1-k}c^{\prime\prime}\right)\right|_{u=U}$, respectively. They are all equal to zero due to Eq. 55 and the fact that $0\leq\beta-1-k\leq\beta-2\leq 2N-5$. With the last relation, the first term in the right-hand side of Eq. 56 is reduced to

[TABLE]

The last step resorts again to Eq. 54b for $j\geq 2N+1-\beta\geq 4$. With this, Eq. 56 is equivalent to

[TABLE]

This is our final task.

In Eq. 57, we take $\mathcal{H}_{j}$ to be $\Delta_{j}$, $\partial_{U}\Delta_{j}$ or $\partial_{\theta}\Delta_{j}$ and arrive at

[TABLE]

*for $1\leq\beta\leq 2N-3$. Here $k=0,1,2$ correspond to $\mathcal{H}_{j}=\Delta_{j},\ \partial_{U}\Delta_{j}\text{ or }\partial_{\theta}\Delta_{j}$, respectively. The last relations hold due to Eq. 55 and the linearity of the operator $\mathcal{D}_{\theta}$. Hence the orthogonality is validated and the proof is completed. *

Remark 4.21**.**

*It is worth pointing out that the structural stability condition still holds if the collision frequency $\nu=\nu(M)$ in the BGK model depends on $M$, because in equilibrium $S=S(M)=0$ and thus $\partial_{M}(\nu S)=\nu S_{M}(M)+S\partial_{M}\nu=\nu S_{M}(M)$. Hence all the analyses above are valid. *

4.5 Shakhov model

This subsection is devoted to the EQMOM moment system of the 1-D Boltzmann equation with Shakhov source term Eq. 6. We first introduce the notation

[TABLE]

with the equilibrium distribution $f_{S}=f_{S}(t,x,\xi)$ defined in Eq. 7. In this situation the moment system has the form Eq. 19 but the source term is different:

[TABLE]

We need to investigate whether this source term satisfies the structural stability condition (i) & (iii). For this purpose, some basic properties of $\Delta_{j}^{S}$ are required. A direct calculation shows that $\Delta_{0}^{S}=\Delta_{0}=1$, $\Delta_{1}^{S}=\Delta_{1}=U$ and $\Delta_{2}^{S}=\Delta_{2}=U^{2}+\theta$. Moreover, for $j\geq 3$ we have

[TABLE]

Here Lemma 4.1(a) is used for the integration and the last step is due to the definition of $q$.

As in Section 4.3, the equilibrium state $W$ needs to be determined. From $S^{Sh}(M)=0$ we see that $\rho\Delta_{j}^{S}=M_{j}$ for $0\leq j\leq 2N$. With Eq. 58, the equation $\rho\Delta_{3}^{S}=M_{3}$ clearly implies that $M_{3}=\rho\Delta_{3}$ for $Pr\neq 1$. Therefore, we have $\rho\Delta_{j}^{S}=\rho\Delta_{j}$ for any $0\leq j\leq 2N$ and the equilibrium manifold $\mathcal{E}$ is determined by $M_{j}=\rho\Delta_{j}$ for any $0\leq j\leq 2N$. This is exactly the same as that of the BGK model, which has already been determined in Theorem 4.17 to be $W\in\Omega_{W}^{eq}$.

At equilibrium, the Jacobian matrix of $S^{Sh}$ can be computed with Eq. 58:

[TABLE]

where $S_{M}$ is the Jacobian matrix Eq. 47 for the BGK model and the $(2N+1)$-matrix $E_{ij}=\left(e_{ij}\right)$ with $e_{ij}=1$ and all the other entities being zero. $S_{M}^{Sh}$ is diagonalizable by an invertible matrix $P^{S}$ such that $P^{S}S_{M}^{Sh}=-\operatorname{diag}(0,0,0,Pr,1,\dots,1)P^{S}$, and

[TABLE]

where $P^{-1}$ is defined in Eq. 51. Hence the structural stability condition (i) is justified.

For Condition (iii), we take the same symmetrizer $A_{0}=L^{T}L$ as that for the BGK model. This is reasonable since the equilibrium state is the same. It then suffices to show that the first three columns of $L\left(P^{S}\right)^{-1}$ are orthogonal to its other columns in equilbrium. From Eq. 60 we see that the only difference between $L\left(P^{S}\right)^{-1}$ and $LP^{-1}$ is the fourth column. For $L\left(P^{S}\right)^{-1}$, its fourth column is a linear combination of the last $(2N-2)$ columns of $L$, while the last $(2N-2)$ columns of $LP^{-1}$ are exactly those of $L$. Since the first three columns of $LP^{-1}$ are orthogonal to its other columns, the fourth column of $L\left(P^{S}\right)^{-1}$ is also orthogonal to its first three columns. This has validated Condition (iii). In this way, we have the main result of this subsection:

Theorem 4.22**.**

*For the 1-D Boltzmann equation with the Shakhov model, the EQMOM moment system satisfies the structural stability condition. *

5 Conclusions

This paper presents a rigorous stability analysis of the quadrature based moment methods (QBMM) for the Boltzmann equation. To figure out a road map for more general cases, only the spatial one-dimensional (1-D) Boltzmann equation with hypothetical collisions (BGK or Shakhov type) is considered here. In the QBMM, the distribution function $f$ is approximated with a linear combination of $N$ ($N\geq 1$) $\delta$-functions with unknown centers or their Gaussian approximations with unknown variance and centers (named QMOM or EQMOM, respectively). For QMOM, we show purely analytically that the resulting moment systems of first-order PDEs are not strongly hyperbolic for any $N$. Furthermore, we prove that the moment systems produced by the Gaussian EQMOM are strictly hyperbolic, when the variance is positive, and preserve the dissipation property of the kinetic equation. As a step in the proof, we also determine the equilibrium manifold that lies on the boundary of the state space for the parameters $(w_{i},u_{i},\sigma^{2})$ ($1\leq i\leq N$). These conclusions explain why the EQMOM gives reasonable numerical results while QMOM does not.

The proofs are quite technical and involve detailed analyses of the characteristic polynomial of the coefficient matrices. They offer a guideline to investigate the multidimensional cases with multiple nodes, which is underway.

Appendix A Injectivity of EQMOM

In this appendix we show

Proposition A.1**.**

*For EQMOM, the map $M=\mathcal{M}(W)$ in Eq. 16 is injective for $W\in\Omega_{W}^{open}$ defined in Eq. 17a. *

Proof A.2**.**

It suffices to demonstrate that the Jacobian matrix $\frac{\partial\mathcal{M}}{\partial W}$ in Eq. 35 is invertible for $W\in\Omega_{W}^{open}$. In fact, we can show that

[TABLE]

for multiple nodes $N\geq 2$. To this end, we set $\mathcal{F}(u)=\left(\Delta_{0}(u,\sigma^{2}),\Delta_{1}(u,\sigma^{2}),\cdots,\Delta_{2N}(u,\sigma^{2})\right)^{T}$ and see from Eq. 35 that

[TABLE]

Denote

[TABLE]

And we see that for each $1\leq i\leq N$,

[TABLE]

and thereby

[TABLE]

Thus, it remains to show

[TABLE]

and

[TABLE]

Note that $f_{N}(u_{1},\dots,u_{N};\sigma)$ is a homogeneous polynomial of $u_{1},\dots,u_{N},\sigma$ with degree $2(N^{2}-1)$. This can be seen from the definition of determinant and the fact that $\Delta_{j}(u,\sigma^{2})$ is a homogeneous polynomial of $u$ and $\sigma$ with degree $j$ (see Lemma 4.1(a)). On the other hand, the right-hand side of Eq. 61 is also a homogeneous polynomial of $u_{1},\dots,u_{N},\sigma$ with degree $2(N^{2}-1)$. Thus, to prove Eq. 61, we need to show that $(u_{j}-u_{1})^{6}$ and $(u_{j}-u_{i})^{4}$ are factors of $f(u_{1},\dots,u_{N};\sigma)$ for any $2\leq j\neq i\leq N$.

From the definition of $f_{N}=f_{N}(u_{1},\dots,u_{N};\sigma)$, it is not difficult to compute that for $j\neq 1$,

[TABLE]

Thus it follows that for $j\neq 1$,

[TABLE]

for any $1\leq i\neq j\leq N$ and that

[TABLE]

This justifies Eq. 61 and $C(N)$ is a constant.

Then all we need is to prove Eq. 62. A direct calculation for $N=2$ indicates that $f_{2}(u_{1},u_{2};\sigma)=(u_{1}-u_{2})^{4}$ and thus $C(2)=1$. For $N>2$, we deduce from Eq. 61 that the leading coefficient of $u_{N}$ (with degree $(4N-2)$) is

[TABLE]

On the other hand, the determinant definition of $f_{N}(u_{1},\dots,u_{N};\sigma)$ implies that the leading term of $u_{N}$ (with degree $(4N-2)$) is included in the following part:

[TABLE]

Thus, using Eq. 61, the leading coefficient of $u_{N}$ is

[TABLE]

*By equating the leading coefficients of $u_{N}$ in the above two expressions, we see immediately that $C(N)=C(N-1)$ for $N>2$. Since $C(2)=1$, this justifies Eq. 62 and hence completes the proof. *