# Stability analysis of an extended quadrature method of moments for kinetic
equations

Ruixi Zhang  Department of Energy and Power Engineering, Tsinghua University
Beijing, 100084, China<EMAIL_ADDRESS>,  Qian Huang  Department of Energy
and Power Engineering, Tsinghua University
Beijing, 100084, China<EMAIL_ADDRESS><EMAIL_ADDRESS>and  Wen-An
Yong  Department of Mathematical Sciences, Tsinghua University
Beijing, 100084, China
Yanqi Lake Beijing Institute of Mathematical Sciences and Applications
Beijing 101408, China<EMAIL_ADDRESS>
###### Abstract.

This paper performs a stability analysis of a class of moment closure systems
derived with an extended quadrature method of moments (EQMOM) for the one-
dimensional BGK equation. The class is characterized with a kernel function. A
sufficient condition on the kernel is identified for the EQMOM-derived moment
systems to be strictly hyperbolic. We also investigate the realizability of
the moment method. Moreover, sufficient and necessary conditions are
established for the two-node systems to be well-defined and strictly
hyperbolic, and to preserve the dissipation property of the kinetic equation.

###### Key words and phrases:

kinetic equation, extended quadrature method of moments, BGK model,
hyperbolicity, structural stability condition

* Corresponding author

##  1\. Introduction

By describing the evolution of problem-specific distribution functions, the
kinetic models are founded with a solid basis in a wide range of complex
interacting systems. For instance, the celebrated Boltzmann equation governs
the distribution of the molecular velocity and is believed to better
characterize the rarefied flows where their hydrodynamic counterparts, being
Euler and/or Navier-Stokes equations, become less reliable [12]. However,
solving the kinetic equation is challenging. For one thing, the binary
collision operator of the Boltzmann equation causes quadratic costs while
treating the velocity dependence. A widely-accepted measure is to apply the
BGK operator which models the collision as a relaxation process towards the
local equilibrium [2]. This model not only reduces the computational costs,
but also has the desired conservation laws and an $H$-theorem characterizing
the dissipation properties due to collisions [12].

On the other hand, the high dimension of the phase space raises significant
difficulties in computation, even for the BGK equation. Several deterministic
methods have thus been developed [9, 18], including the discrete velocity
model, the spectral methods, and various methods of moments, to remove the
dependence upon the molecular velocity and deduce spatial-time models of
macroscopic variables from the kinetic equation. This paper focuses on the
method of moments. In the moment system, only the lower-order moments have
clear physical interpretations after being related to the density, mean
velocity and temperature (internal energy) of the system [12]. The higher-
order moments may provide additional information beyond the classical
hydrodynamic models.

All moment systems need to be closed, which is mostly done by reconstructing
the distribution from the transported moments [16]. The well-known Grad’s
13-moment theory was established based on a linear expansion of the Maxwellian
[11], but far away from equilibrium, this ansatz can lead to negative values
of distributions. In contrast, the quadrature-based method of moments [17] was
proposed with nonlinear reconstruction of the distribution function, and such
a treatment seems to be suitable for non-equilibrium flows. However,
approximating the distribution with a linear combination of multiple Dirac
$\delta$-functions with unknown centers [17], the method (called QMOM) fails
in the simulation of BGK equation due to the occurrence of singularity [10].

This highlights the importance of understanding the mathematical properties of
the moment closure systems, which are usually first-order PDEs. For real-world
physical models, the PDE is expected to be hyperbolic so that the system is
robust against small perturbations of the initial data [22]. Indeed, the
unphysical behaviors of the Grad’s theory and the QMOM approach can both be
attributed to the lack of hyperbolicity [3, 5]. Let us mention some efforts to
achieve hyperbolic regulation of the Grad’s theory [3, 23] and other moment
closure systems [15]. Furthermore, a thorough stability analysis should as
well account for the source term of the model, as the Boltzmann and BGK
collisions are both featured with $H$-theorems [12]. For the moment closure
system, it is believed that the structural stability condition proposed in
[25] for hyperbolic relaxation systems is a proper characterization of the
dissipation property. The condition specifies how the source term should be
coupled with the hyperbolic part in the vicinity of the equilibrium. Admitting
such a structure, the resultant moment system is compatible with the classical
theories [25]. Recently, it was shown that the structural stability condition
is fulfilled by many moment closure systems, including the hyperbolic
regularization models of rarefied gases [8, 28] and a series of hyperbolic
shallow water moment models [14].

The objective of this paper is to investigate the stability properties of the
extended QMOM (termed EQMOM) for the BGK equation. In EQMOM, the velocity
distribution is reconstructed as a sum of multiple continuous kernel functions
instead of the $\delta$-function [4, 27]. A list of kernels that can be used
for the purpose of EQMOM is summarized in [20]. If the kernel is the Gaussian
distribution (denoted Gaussian-EQMOM), it was found that the method is well-
defined (namely, the unclosed terms can be uniquely determined), and the
resultant moment system respects the structural stability condition [13].
However, little is known for other types of kernels, which may be suitable for
different scenarios. For instance, the space plasmas follow kappa
distributions with high energy tails deviated from a Maxwellian [19]. We also
remark that the EQMOM with the beta function as the kernel is used to treat a
simplified radiative transfer equation [1].

This paper deals with the EQMOM induced by a univariate kernel function (see
(2.4)) and considers the spatial one-dimensional (1-D) BGK equation.
Surprisingly, it is found that not all kernels share the nice properties of
the Gaussian kernel. As the main result of this paper, we reveal the exact
constraints on the kernel such that the two-node EQMOM ($n=2$) is well-
defined, that the resultant moment system is strictly hyperbolic, and that the
system properly respects the dissipation property (as judged by the structural
stability condition). These constraints are inequalities of the moments of the
kernel function, and are thus easy to check for specific kernels. Moreover, a
sufficient condition is estalished for the $n$-node EQMOM-derived moment
system to be strictly hyperbolic. It is worth mentioning that, in our
argument, the realizable moment set can be characterized for general $n\geq
1$, which essentially contains the results in [6] (Proposition 3.1 for the
two-node Gaussian-EQMOM) as a special case. This leads to a polynomial which
is an answer to the question raised in [16] (footnote 16 on p.88).
Furthermore, we present abundant examples of the kernel functions to which our
theory can be applied in a straightforward manner.

The remainder of the paper is organized as follows. Our main results are
presented in Section 2. Section 3 is devoted to studying the realizability of
the extended quadrature method of moments. Hyperbolicity of the resultant
moment systems is analyzed in Section 4. For two-node systems, the structural
stability condition is verified in Section 5. Section 6 presents a number of
specific kernel functions and investigates their behaviors numerically. The
conclusions are given in Section 7.

##  2\. Preliminaries and main results

Consider the hypothetical 1-D BGK equation for the distribution $f=f(t,x,\xi)$
with time $t\in\mathbb{R}_{+}$, position $x\in\mathbb{R}$ and velocity
$\xi\in\mathbb{R}$:

$\left\\{\begin{aligned}
\partial_{t}f+\xi\partial_{x}f&=\frac{1}{\tau}(f^{eq}-f),\\\
f^{eq}=f^{eq}(\rho,U,\theta;\xi)&=\frac{\rho}{\sqrt{2\pi\theta}}\exp\left(-\frac{(\xi-U)^{2}}{2\theta}\right).\end{aligned}\right.$
(2.1)

Here $\tau$ is a relaxation time. In the local Maxwellian equilibrium
$f^{eq}$, the density $\rho$, mean velocity $U$ and temperature $\theta$ are
determined as the velocity moments of $f$:

$\rho=\int_{\mathbb{R}}fd\xi,\quad\rho U=\int_{\mathbb{R}}\xi
fd\xi,\quad\rho\theta+\rho U^{2}=\int_{\mathbb{R}}\xi^{2}fd\xi.$

Define the $j$th velocity moments of $f$ as

$M_{j}=M_{j}(t,x)=\int_{\mathbb{R}}\xi^{j}fd\xi$

for $j\in\mathbb{N}$. The evolution equation for $M_{j}$ can be immediately
derived from (2.1) as

$\partial_{t}M_{j}+\partial_{x}M_{j+1}=\frac{1}{\tau}\left(\rho\Delta_{j}^{eq}(U,\sqrt{\theta})-M_{j}\right)$
(2.2)

with

$\Delta_{j}^{eq}(u,\sigma)=\int_{\mathbb{R}}\xi^{j}\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(\xi-u)^{2}}{2\sigma^{2}}\right)d\xi.$

Notice that the first $n$ equations for $M_{0},...,M_{n-1}$ contain the term
$\partial_{x}M_{n}$. Therefore, any finite truncation of the above equations
leads to an unclosed system, and a closure procedure is required. In this
paper, we are concerned with an extended quadrature method of moments (EQMOM)
[4, 27].

###  2.1. Extended quadrature method of moments (EQMOM)

Let $\mathcal{K}=\mathcal{K}(\xi)\geq 0$ satisfy

$\mathfrak{m}_{j}:=\int_{\mathbb{R}}\xi^{j}\mathcal{K}(\xi)d\xi<\infty,\
\forall j\in\mathbb{N},\quad\text{and }\mathfrak{m}_{0}=1.$ (2.3)

In EQMOM, the distribution $f$ is approximated with the following ansatz

$f(\xi)=\sum_{i=1}^{n}\frac{w_{i}}{\sigma}\mathcal{K}\left(\frac{\xi-
u_{i}}{\sigma}\right)$ (2.4)

with the weights $w_{i}$, nodes $u_{i}$ and ‘width’ $\sigma>0$ to be
determined. To do this, the first ($2n+1$) lower-order moments are employed:

$M_{j}=\sum_{i=1}^{n}w_{i}\Delta_{j}(u_{i},\sigma),\quad\textrm{with
}\Delta_{j}(u,\sigma):=\int_{\mathbb{R}}\xi^{j}\frac{1}{\sigma}\mathcal{K}\left(\frac{\xi-u}{\sigma}\right)d\xi$
(2.5)

for $j=0,\dots,2n$. This defines a map
$M:=(M_{0},...,M_{2n})^{T}=\mathcal{M}(W)$ for
$W=(w_{1},u_{1},...,w_{n},u_{n},\sigma)^{T}\in\mathbb{R}^{2n+1}$ with
$\sigma>0$. Here the superscript ‘$T$’ denotes the transpose of a vector or
matrix.

Suppose $\mathcal{M}$ is injective on a certain domain
$W\in\Omega\subset\mathbb{R}^{2n+1}$. Then for any
$M\in\mathbb{G}:=\mathcal{M}(\Omega)$, there exists a unique
$W=\mathcal{M}^{-1}(M)$ solving (2.5). In this way, the EQMOM is well-defined
and the next moment $M_{2n+1}$ can be evaluated as a function of $M$:

$M_{2n+1}=\mathcal{M}_{2n+1}(M):=\sum_{i=1}^{n}w_{i}\Delta_{2n+1}(u_{i},\sigma).$
(2.6)

Consequently, the first ($2n+1$) equations in (2.2) are closed as a system of
first-order PDEs:

$\partial_{t}M+\partial_{x}(M_{1},...,M_{2n+1})^{T}=\frac{1}{\tau}\left(\rho\Delta^{eq}(U,\sqrt{\theta})-M\right).$
(2.7)

Here
$\Delta^{eq}(U,\sqrt{\theta})=\left(\Delta_{0}^{eq}(U,\sqrt{\theta}),\dots,\Delta_{2n}^{eq}(U,\sqrt{\theta})\right)^{T}\in\mathbb{R}^{2n+1}$
with $\rho=M_{0}$, $\rho U=M_{1}$ and $\rho(\theta+U^{2})=M_{2}$.

The main goal of this paper is to investigate the injectivity of the map
$\mathcal{M}$ in (2.5) for the general kernel $\mathcal{K}(\xi)$, where the
injectivity is closely related to the realizability of moments. Such analyses
are useful to design efficient and robust algorithms to solve $W$ from (2.5).
Moreover, we analyze hyperbolicity of the moment closure system (2.7) and its
dissipation property inherited from the the $H$-theorem of the kinetic
equation.

######  Remark 2.1.

For the sake of simplicity, throughout this paper we assume that the kernel
function $\mathcal{K}(\xi)$ is normalized, in the sense that
$\mathfrak{m}_{0}=1$, $\mathfrak{m}_{1}=0$ and $\mathfrak{m}_{2}=1$. In fact,
a general $\mathcal{K}(\xi)$ can be normalized as
$\mathcal{K}^{+}(\xi)=h\mathcal{K}(h\xi+\xi_{0})$ with
$\xi_{0}=\mathfrak{m}_{1}$ and
$h=\sqrt{\mathfrak{m}_{2}-\mathfrak{m}_{1}^{2}}$ (note that
$\mathfrak{m}_{2}\geq\mathfrak{m}_{1}^{2}$ due to the Cauchy-Schwartz
inequality). Moreover, the map $\mathcal{M}_{2n+1}(M)$ in (2.6) is the same as
that derived from $\mathcal{K}(\xi)$, so the moment closure system (2.7) is
unchanged with the normalization.

###  2.2. Structural stability condition

For smooth solutions, the balance laws (2.7) can be written as

$\partial_{t}M+A(M)\partial_{x}M=S(M):=\frac{1}{\tau}\left(\rho\Delta^{eq}(U,\sqrt{\theta})-M\right)$
(2.8)

with coefficient matrix

$A(M)=\begin{bmatrix}0&1&&&\\\ &0&1&&\\\ &&\ddots&\ddots&\\\ &&&0&1\\\
a_{0}&a_{1}&\cdots&a_{2n-1}&a_{2n}\end{bmatrix},$ (2.9)

where $a_{j}(M)=\partial\mathcal{M}_{2n+1}/\partial M_{j}$ for $j=0,\dots,2n$.
It is called hyperbolic if $A(M)\in\mathbb{R}^{(2n+1)\times(2n+1)}$ has
$(2n+1)$ linearly-independent real eigenvectors [22]. If $A(M)$ has $(2n+1)$
distinct real eigenvalues, it is called strictly hyperbolic. Obviously, strict
hyperbolicity implies hyperbolicity. The dissipativeness of the moment system
can be characterized with the structural stability condition proposed in [25]
for hyperbolic relaxation systems.

Assume that the equilibrium manifold $\mathcal{E}=\\{M\in\mathbb{G}\mid
S(M)=0\\}$ is not empty. Denote by $S_{M}(M)$ the Jacobian matrix of $S(M)$.
The structural stability condition reads as

  * (I)

For any $M\in\mathcal{E}$, there exist invertible matrices
$P=P(M)\in\mathbb{R}^{(2n+1)\times(2n+1)}$ and
$\hat{T}=\hat{T}(M)\in\mathbb{R}^{r\times r}$ ($0<r\leq 2n+1$) such that

$PS_{M}(M)P^{-1}=\text{diag}(\textbf{0}_{(2n+1-r)\times(2n+1-r)},\hat{T}).$

  * (II)

For any $M\in\mathbb{G}$, there exists a positive definite symmetric matrix
$A_{0}=A_{0}(M)$ such that $A_{0}A(M)=A^{T}(M)A_{0}$.

  * (III)

For any $M\in\mathcal{E}$, the coefficient matrix and the source are coupled
as

$A_{0}S_{M}(M)+S_{M}^{T}(M)A_{0}\leq-P^{T}\begin{bmatrix}0&0\\\
0&I_{r}\end{bmatrix}P.$

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

######  Remark 2.2.

Recently, it has been demonstrated that several moment models from the kinetic
equations respect the structural stability condition, including the Gaussian-
EQMOM [13] and the hyperbolic regularization models [8, 28]. For the 1-D
system (2.8), Condition (II) is satisfied if and only if the system is
hyperbolic [13]. Condition (III) can be regarded as a proper manifestation of
the dissipation property inherited from the kinetic model. See detailed
discussions in [26].

###  2.3. Main results

Our main results are collected in this subsection. As mentioned in Remark 2.1,
we assume that the kernel $\mathcal{K}(\xi)$ is normalized, that is,
$\mathfrak{m}_{0}=\mathfrak{m}_{2}=1$ and $\mathfrak{m}_{1}=0$.

To state the results, we recursively define a sequence of numbers associated
with $\\{\mathfrak{m}_{j}\\}:$

$b_{0}=1,\quad
b_{j}=-\sum_{k=1}^{j}\frac{\mathfrak{m}_{k}}{k!}b_{j-k}\quad\text{for
}j=1,2,\dots$ (2.10)

and a number of auxiliary moments associated with
$M=(M_{0},...,M_{2n})^{T}\in\mathbb{R}^{2n+1}$:

$M^{*}=\\{M_{j}^{*}\\}_{j=0}^{2n},\quad
M_{j}^{*}=M_{j}^{*}(\sigma)=\sum_{k=0}^{j}b_{k}\sigma^{k}\frac{j!}{(j-k)!}M_{j-k}.$
(2.11)

Moreover, for $M=(M_{0},...,M_{2n})^{T}$ we introduce the Hankel matrix [21]
as

$H_{k}(M)=\begin{bmatrix}M_{0}&M_{1}&\cdots&M_{k}\\\
M_{1}&M_{2}&\cdots&M_{k+1}\\\ \vdots&\vdots&\ddots&\vdots\\\
M_{k}&M_{k+1}&\cdots&M_{2k}\end{bmatrix}\in\mathbb{R}^{(k+1)\times(k+1)}\quad\text{for
}k\leq n.$ (2.12)

This is a real symmetric matrix.

Our first result is

######  Theorem 2.3.

Given $M=(M_{0},\dots,M_{2n})^{T}\in\mathbb{R}^{2n+1}$, set

$P_{n}(\sigma;M)=\det H_{n}(M^{*}(\sigma)).$ (2.13)

The following statements are equivalent.

(i). There exists a unique $W\in\Omega=\Omega^{\prime}\times\\{\sigma>0\\}$
with

$\Omega^{\prime}=\\{(w_{1},u_{1},...,w_{n},u_{n})\in\mathbb{R}^{2n}\mid\
w_{i}>0,\ \forall i;\ u_{1}<\dots<u_{n}\\},$ (2.14)

such that $\mathcal{M}(W)=M$.

(ii). $P_{n}(\sigma;M)=0$ has a unique positive root $\sigma_{0}$ such that
the Hankel matrix $H_{n-1}(M^{*}(\sigma_{0}))$ is positive definite.

######  Remark 2.4.

Clearly, a similar conclusion can be formulated when some of the weights
$w_{i}$ are zero or the centers $u_{i}$ coincide. In that case, one can find a
unique index $k\leq n$ such that $w_{i}>0$ for any $i=1,\dots,k$ and all the
$u_{i}$’s are distinct.

######  Remark 2.5.

Statement (ii) serves as an implicit realizable condition for $M$. For $n=2$,
this condition is exactly that in [6] for the Gaussian kernel. Moreover,
Theorem 2.3 suggests a key step in inverting the map $M=\mathcal{M}(W)$,
namely, finding $\sigma$ as a root of the polynomial $P_{n}(\sigma)=\det
H_{n}(M^{*}(\sigma))$ of degree $n(n+1)$. For $n=2$, the polynomial for even
and normalized kernels is

$P_{2}(\sigma)=(5-\mathfrak{m}_{4})\sigma_{1}^{3}+2\theta(3-\mathfrak{m}_{4})\sigma_{1}^{2}+(M_{4}^{\prime}-\mathfrak{m}_{4}\theta^{2})\sigma_{1}+M_{3}^{\prime
2}$

with $\sigma_{1}=\sigma^{2}-\theta$,

$M_{3}^{\prime}=\frac{M_{3}}{M_{0}}-3U\theta-U^{3},\quad
M_{4}^{\prime}=\frac{M_{4}}{M_{0}}-4UM_{3}^{\prime}-6U^{2}\theta-U^{4},$

$U=M_{1}/M_{0}$ and $\theta=M_{2}/M_{0}-U^{2}$. Once $\sigma$ is found, other
components of $W$ can be determined by the existing algorithms [16] (see
Remark 3.1).

As a corollary of this theorem, we have

######  Corollary 2.6 (Injectivity).

For $n=2$, the map $M=\mathcal{M}(W)$ in (2.5) is injective for $W\in\Omega$
if and only if the inequality $\mathfrak{m}_{4}\geq
3+\frac{9}{8}\mathfrak{m}_{3}^{2}$ holds for $\mathcal{K}(\xi)$.

######  Remark 2.7.

For $n=2$, it is not difficult to see that $M_{5}$ can be expressed in terms
of $M_{0},...,M_{4}$ when $u_{1}=u_{2}$. As a consequence, the condition of
Corollary 2.6 ensures that the EQMOM is well defined on
$\mathcal{M}(\Omega^{tot})$ with

$\Omega^{tot}=\\{W\in\mathbb{R}^{5}|w_{1}>0,\ w_{2}>0,\ u_{1}\leq u_{2},\
\sigma>0\\}.$

Suppose the map $M=\mathcal{M}(W)$ is injective on $\Omega$ for general $n$.
Our second result is

######  Theorem 2.8 (Hyperbolicity).

If the $b$-polynomial

$p(t):=\sum_{j=0}^{2n+1}b_{2n+1-j}t^{j}$

has $(2n+1)$ real roots (counting multiplicity) and at least two roots are
nonzero, then the $n$-node EQMOM moment system (2.7) is strictly hyperbolic on
$\mathcal{M}(\Omega)$.

Furthermore, for even kernels we have

######  Theorem 2.9 (Hyperbolicity).

Let the kernel $\mathcal{K}(\xi)$ be an even function. The two-node EQMOM
moment system (2.7) is strictly hyperbolic for $M\in\mathcal{M}(\Omega^{tot})$
if and only if $3\leq\mathfrak{m}_{4}<6$.

######  Theorem 2.10 (Dissipativeness).

Let the kernel $\mathcal{K}(\xi)$ be an even function and
$3\leq\mathfrak{m}_{4}<6$. Then the two-node EQMOM moment system (2.7)
satisfies the structural stability condition if and only if
$3\leq\mathfrak{m}_{4}<5$.

In the next section, Theorem 2.3 and Corollary 2.6 will be proved. Section 4.1
is devoted to a proof of Theorem 2.8, while Theorems 2.9 & 2.10 are proved in
Sections 4.2 and 5, respectively.

##  3\. Injectivity

In this section, we prove Theorem 2.3 and Corollary 2.6. To start with, we
recall the definition of the map $M=\mathcal{M}(W)$ in (2.5):

$M_{j}=\sum_{i=1}^{n}w_{i}\Delta_{j}(u_{i},\sigma),\quad\textrm{with
}\Delta_{j}(u,\sigma):=\int_{\mathbb{R}}\xi^{j}\frac{1}{\sigma}\mathcal{K}\left(\frac{\xi-u}{\sigma}\right)d\xi,$

for $j=0,\dots,2n$.

By performing the change of variables $\frac{\xi-u}{\sigma}\mapsto\xi$, we can
easily see that

$\Delta_{j}(u,\sigma)=\sum_{k=0}^{j}\binom{j}{k}\mathfrak{m}_{k}\sigma^{k}u^{j-k},$
(3.1)

indicating that $\Delta_{j}(u,\sigma)$ is a homogeneous bivariate polynomial
of $u$ and $\sigma$. It is not difficult to verify

$\displaystyle\partial_{u}\Delta_{j}(u,\sigma)$
$\displaystyle=j\Delta_{j-1}(u,\sigma),$ (3.2a)
$\displaystyle\partial_{\sigma}\Delta_{j}(u,\sigma)$
$\displaystyle=j\sum_{k=0}^{j-1}\binom{j-1}{k}\mathfrak{m}_{k+1}\sigma^{k}u^{j-1-k}.$
(3.2b)

Notice that $M=(M_{0},...,M_{2n})^{T}$ can be conversely expressed in terms of
the auxiliary moments $M_{j}^{*}=M_{j}^{*}(\sigma)$ defined in (2.11) as

$M_{j}=\sum_{k=0}^{j}\binom{j}{k}\mathfrak{m}_{k}\sigma^{k}M_{j-k}^{*}\quad\text{for
}j=0,\dots,2n.$ (3.3)

Indeed, a straightforward calculation of the right-hand side, incorporating
(2.11), yields

r.h.s.
$\displaystyle=\sum_{k=0}^{j}\binom{j}{k}\mathfrak{m}_{k}\sigma^{k}\sum_{l=0}^{j-k}b_{l}\sigma^{l}\frac{(j-k)!}{(j-k-l)!}M_{j-k-l}$
$\displaystyle=\sum_{s=0}^{j}\frac{j!}{(j-s)!}\sigma^{s}M_{j-s}\sum_{k=0}^{s}\frac{\mathfrak{m}_{k}}{k!}b_{s-k}.$

The second equality is derived after the change of variables $s=k+l$.
Rewriting (2.10) as

$\sum_{k=0}^{s}\frac{\mathfrak{m}_{k}}{k!}b_{s-k}=0,\quad s\geq 1,$

we obtain $\text{r.h.s}=M_{j}$. Similarly, with (3.1) involved, a direct
calculation of the right-hand side of (2.11) results in

$M_{j}^{*}(\sigma)=\sum_{i=1}^{n}w_{i}u_{i}^{j},\quad j=0,...,2n.$ (3.4)

######  Remark 3.1.

Together with Remark 2.5, the last formula suggests a practical method to
solve $W\in\mathbb{R}^{2n+1}$ from the EQMOM map $M=\mathcal{M}(W)$. Once
$\sigma$ has been determined as shown in Remark 2.5, the $(w_{i},u_{i})$’s can
be determined by solving the first $2n$ equations in (3.4) with existing
algorithms [16].

About the Hankel matrix, we quote the following lemma.

######  Lemma 3.2 ([21], Theorem 9.7).

Given $M^{\prime}=(M_{0},\dots,M_{2n-1})\in\mathbb{R}^{2n}$, if the nonlinear
equations

$\sum_{i=1}^{n}w_{i}u_{i}^{j}=M_{j}\quad\text{for }j=0,\dots,2n-1$

have a solution in $\Omega^{\prime}$ defined in (2.14), then the Hankel matrix
$H_{n-1}(M^{\prime})$ is positive definite. Conversely, if
$H_{n-1}(M^{\prime})$ is positive definite, then the last equations have a
unique solution in $\Omega^{\prime}$.

To prove Theorem 2.3, we first notice the following fact:

######  Proposition 3.3.

If $M_{j}^{*}=\sum_{i=1}^{n}w_{i}u_{i}^{j}$ for $j=0,\dots,2n$, the Hankel
matrix $H_{n}(\\{M_{j}^{*}\\})$ is singular.

###### Proof.

A direct calculation gives

$\det\begin{bmatrix}M_{0}^{*}&\cdots&M_{n}^{*}\\\ \vdots&&\vdots\\\
M_{n}^{*}&\cdots&M_{2n}^{*}\end{bmatrix}=\sum_{1\leq i_{0},...,i_{n}\leq
n}\det\begin{bmatrix}w_{i_{0}}&\cdots&w_{i_{n}}u_{i_{n}}^{n}\\\
\vdots&&\vdots\\\
w_{i_{0}}u_{i_{0}}^{n}&\cdots&w_{i_{n}}u_{i_{n}}^{2n}\end{bmatrix}.$

For each determinant in the summation, at least two of the $(n+1)$ indices
$1\leq i_{0},\dots,i_{n}\leq n$ are identical, therefore the determinant is
zero. Hence the Hankel matrix is singular. ∎

###### Proof of Theorem 2.3.

(i) $\Rightarrow$ (ii). In this case, it has been shown in (3.4) that
$M^{*}_{j}(\sigma)$ can be expressed as $\sum_{i=1}^{n}w_{i}u_{i}^{j}$ for
$j=0,\dots,2n$ with $w_{i}>0$ and all the $u_{i}$’s distinct. Then it follows
from Proposition 3.3 that $P_{n}(\sigma;M)=\det
H_{n}(\\{M_{j}^{*}(\sigma)\\})=0$. Because all $w_{i}>0$ and the $u_{i}$’s are
distinct, we deduce from Lemma 3.2 that $H_{n-1}(\\{M_{j}^{*}(\sigma)\\})$ is
positive definite.

For the uniqueness, suppose that $P_{n}(\sigma;M)$ has another root
$\sigma_{1}>0$ such that $H_{n-1}(\\{M_{j}^{*}(\sigma_{1})\\})$ is positive
definite. It follows from Lemma 3.2 that there exists a $2n$-tuple
$\\{q_{i}>0,v_{i}\\}_{1\leq i\leq n}$ such that $v_{1}<\cdots<v_{n}$ and
$M_{j}^{*}(\sigma_{1})=\sum_{i=1}^{n}q_{i}v_{i}^{j}$ for $j=0,...,2n-1$. Set
$S_{j}=\sum_{i=1}^{n}q_{i}v_{i}^{j}$. From $P_{n}(\sigma_{1};M)=0$ and
Proposition 3.3 we see that $\det
H_{n}(\\{M_{j}^{*}(\sigma_{1})\\})=P_{n}(\sigma_{1};M)=\det
H_{n}(\\{S_{j}\\})=0$. On the other hand, we observe that

$\det H_{n}(\\{M_{j}^{*}(\sigma_{1})\\})=aM_{2n}^{*}(\sigma_{1})+b,\quad\det
H_{n}(\\{S_{j}\\})=aS_{2n}+b$

with

$a=\det H_{n-1}(\\{M_{j}^{*}\\})=\det H_{n-1}(\\{S_{j}\\})>0$

and $b$ depending only on $M^{*}_{j}(\sigma_{1})=S_{j}$ with $j\leq 2n-1$.
Thus, we have $M^{*}_{2n}(\sigma_{1})=S_{2n}$ and thereby get another solution
to the equations $\mathcal{M}(W^{\prime})=M$, violating the uniqueness in (i).
This proves (ii).

(ii) $\Rightarrow$ (i). Assume that $P_{n}(\sigma_{0};M)=0$ and
$H_{n-1}(\\{M_{j}^{*}(\sigma_{0})\\})$ is positive definite. The reasoning
above shows that there exists a unique $2n$-tuple $\\{w_{i}>0,u_{i}\\}$ such
that $u_{1}<\cdots<u_{n}$ and $W=(w_{i},u_{i},\sigma_{0})$ solves
$\mathcal{M}(W)=M$. If $W_{1}=(q_{i},v_{i},\sigma_{1})\neq W$ is another
solution, then $\sigma_{1}\neq\sigma_{0}$ and the reasoning in (i) shows that
$\sigma_{1}$ is another root of $P_{n}(\sigma;M)=0$ which contradicts (ii).
This completes the proof. ∎

Now we are in a position to prove Corollary 2.6.

###### Proof of Corollary 2.6.

Assume $\mathfrak{m}_{4}\geq 3+\frac{9}{8}\mathfrak{m}_{3}^{2}$ for
$\mathcal{K}(\xi)$. It suffices to show that the Jacobian
$\frac{\partial\mathcal{M}}{\partial W}$ is invertible for $W\in\Omega$.
Recall the explicit expressions of $\Delta_{j}(u,\sigma)$ in (3.1) and its
derivatives in (3.2a) & (3.2b). Using

$M_{j}=w_{1}\Delta_{j}(u_{1},\sigma)+w_{2}\Delta_{j}(u_{2},\sigma),\quad
j=0,1,\dots,4,$

we compute the $(j+1)$th row of the Jacobian
$\frac{\partial\mathcal{M}}{\partial W}$ as

$\left(\Delta_{j}(u_{1},\sigma),jw_{1}\Delta_{j-1}(u_{1},\sigma),\Delta_{j}(u_{2},\sigma),jw_{2}\Delta_{j-1}(u_{2},\sigma),\sum_{i=1}^{2}w_{i}\partial_{\sigma}\Delta_{j}(u_{i},\sigma)\right)$

and, by resorting to MATLAB,

$\det\frac{\partial\mathcal{M}}{\partial
W}=w_{1}w_{2}\sigma^{3}(u_{1}-u_{2})^{4}\left[w_{1}q\left(\frac{u_{1}-u_{2}}{\sigma}\right)+w_{2}q\left(\frac{u_{2}-u_{1}}{\sigma}\right)\right]$

with $q(x)=x^{2}+3\mathfrak{m}_{3}x+2(\mathfrak{m}_{4}-3)$. Obviously, we have
$q(x)\geq 0$ for any real $x$ due to $\mathfrak{m}_{4}\geq
3+\frac{9}{8}\mathfrak{m}_{3}^{2}$, and $q(x)$ does not have two distinct
zeros. Therefore, we have $\det\frac{\partial\mathcal{M}}{\partial W}>0$ for
$W\in\Omega$.

Conversely, we assume $\mathfrak{m}_{4}<3+\frac{9}{8}\mathfrak{m}_{3}^{2}$ for
$\mathcal{K}(\xi)$ and show that $\mathcal{M}$ is not injective on $\Omega$.
To do this, we set
$M_{b}=\mathcal{M}(\frac{1}{2},0,\frac{1}{2},0,1)=\\{\mathfrak{m}_{j}\\}_{j=0}^{4}$
and compute

$\displaystyle M_{0}^{*}=1,\quad M_{1}^{*}=0,\quad
M_{2}^{*}(\sigma)=1-\sigma^{2},$ $\displaystyle
M_{3}^{*}(\sigma)=\mathfrak{m}_{3}(1-\sigma^{3}),\quad
M_{4}^{*}(\sigma)=\mathfrak{m}_{4}(1-\sigma^{4})-6\sigma^{2}(1-\sigma^{2})$

and

$\displaystyle\frac{P_{2}(\sigma;M_{b})}{(1-\sigma)^{2}}=$
$\displaystyle(\mathfrak{m}_{4}-\mathfrak{m}_{3}^{2}-5)\sigma^{4}+2(\mathfrak{m}_{4}-\mathfrak{m}_{3}^{2}-5)\sigma^{3}+(2\mathfrak{m}_{4}-3\mathfrak{m}_{3}^{2}-6)\sigma^{2}$
$\displaystyle+2(\mathfrak{m}_{4}-\mathfrak{m}_{3}^{2}-1)\sigma+(\mathfrak{m}_{4}-\mathfrak{m}_{3}^{2}-1).$

The corresponding $2\times 2$ Hankel matrix, which we rewrite as
$H_{1}(\sigma;M_{b})$, is positive definite if and only if the polynomial

$\det
H_{1}(\sigma;M_{b})=M^{*}_{0}(\sigma)M^{*}_{2}(\sigma)-M^{*2}_{1}(\sigma)=1-\sigma^{2}>0,$

namely, $0<\sigma<1$. On the other hand, we denote
$\tilde{P}_{2}(\sigma;M)=\frac{P_{2}(\sigma;M)}{(1-\sigma)^{2}}$ and notice

$\tilde{P}_{2}(1;M_{b})=8\mathfrak{m}_{4}-9\mathfrak{m}_{3}^{2}-24<0\quad\text{and}\quad\tilde{P}_{2}(0;M_{b})=\mathfrak{m}_{4}-\mathfrak{m}_{3}^{2}-1=\det
H_{2}(M_{b})>0$

where the second inequality follows from the positivity of $H_{2}(M_{b})$. By
continuity, we may choose $\epsilon>0$ and
$M\in\mathcal{M}(\Omega^{\prime}\times\\{\sigma=1\\})$ such that

$P_{2}(\epsilon;M)>0>P_{2}(1-\epsilon;M)\quad\text{and}\quad\det
H_{1}(\sigma;M)>0$

for $\sigma\in[\epsilon,1-\epsilon]$. Therefore, $P_{2}(\sigma;M)$ has a root
$\sigma_{0}$ in $(\epsilon,1-\epsilon)$. Clearly, $\sigma=1$ is another root
of $P_{2}(\sigma;M)$ such that $H_{1}(1;M)>0$. By Theorem 2.3, $\mathcal{M}$
is not injective and hence the proof is complete. ∎

##  4\. Hyperbolicity

In this section we prove Theorems 2.8 & 2.9. To this end, it suffices to show
that the characteristic polynomial

$c=c(u;M)=u^{2n+1}-\sum_{j=0}^{2n}a_{j}u^{j}$ (4.1)

of the coefficient matrix $A(M)$ in (2.9) has distinct real roots.

###  4.1. A proof of Theorem 2.8

First of all, we follow [13] and introduce an auxiliary polynomial associated
with the characteristic polynomial:

$g=g(u;W)=\Delta_{2n+1}(u,\sigma)-\sum_{j=0}^{2n}a_{j}\Delta_{j}(u,\sigma).$
(4.2)

With $g$ thus defined, we claim that

$c(u;M)=\sum_{k=0}^{2n+1}b_{k}\sigma^{k}\partial_{u}^{k}g(u;W).$ (4.3)

Indeed, we recall 2.11 & 3.4 that

$\sum_{i=1}^{n}w_{i}u_{i}^{j}=M_{j}^{*}=\sum_{k=0}^{j}b_{k}\sigma^{k}\frac{j!}{(j-k)!}\sum_{i=1}^{n}w_{i}\Delta_{j-k}(u_{i},\sigma).$

In this identity, we set $w_{1}=1,w_{2}=\cdots=w_{n}=0$ and derive from (3.2a)
that

$u_{1}^{j}=\sum_{k=0}^{j}b_{k}\sigma^{k}\frac{j!}{(j-k)!}\Delta_{j-k}(u_{1},\sigma)=\sum_{k=0}^{j}b_{k}\sigma^{k}\partial_{u}^{k}\Delta_{j}(u_{1},\sigma).$

Thus the claim becomes clear.

Furthermore, $g(u;W)$ has the following elegant property for general kernel
$\mathcal{K}(\xi)$.

######  Proposition 4.1.

$g(u;W)=(u-u_{1})^{2}\cdots(u-u_{n})^{2}(u-\tilde{u}(W)).$

###### Proof.

With the expressions for $M_{j}$ in (2.5) and $\Delta_{j}(u,\sigma)$ in (3.1)
($j=0,\dots,2n$), we calculate the $(j+1)$th-row of the Jacobian matrix
$\frac{\partial\mathcal{M}}{\partial W}\in\mathbb{R}^{(2n+1)\times(2n+1)}$ as

$\left(\Delta_{j}(u_{1},\sigma),w_{1}\partial_{u}\Delta_{j}(u_{1},\sigma),\dots,\Delta_{j}(u_{n}),w_{n}\partial_{u}\Delta_{j}(u_{n},\sigma),\partial_{\sigma}M_{j}\right).$

Moreover, from (2.6) we compute $\frac{\partial\mathcal{M}_{2n+1}}{\partial
W}$ as

$\left(\Delta_{2n+1}(u_{1},\sigma),w_{1}\partial_{u}\Delta_{2n+1}(u_{1},\sigma),\dots,\Delta_{2n+1}(u_{n},\sigma),w_{n}\partial_{u}\Delta_{2n+1}(u_{n},\sigma),\partial_{\sigma}M_{2n+1}\right).$

Then we use the relation

$(a_{0},\dots,a_{2n})\frac{\partial\mathcal{M}}{\partial
W}=\frac{\partial\mathcal{M}_{2n+1}}{\partial
M}\frac{\partial\mathcal{M}}{\partial
W}=\frac{\partial\mathcal{M}_{2n+1}}{\partial W}$ (4.4)

to obtain, for $i=1,\dots,n$,

$\displaystyle
a_{0}\Delta_{0}(u_{i},\sigma)+\cdots+a_{2n}\Delta_{2n}(u_{i},\sigma)$
$\displaystyle=\Delta_{2n+1}(u_{i},\sigma),$ (4.5a) $\displaystyle
a_{0}\partial_{u}\Delta_{0}(u_{i},\sigma)+\cdots+a_{2n}\partial_{u}\Delta_{2n}(u_{i},\sigma)$
$\displaystyle=\partial_{u}\Delta_{2n+1}(u_{i},\sigma),$ (4.5b)

implying $g(u_{i})=\partial_{u}g(u_{i})=0$ for $i=1,\dots,n$. This immediately
leads to the expected expression. ∎

We also need the following elementary fact.

######  Proposition 4.2.

If a polynomial $f(u)$ of degree $N$ has $N$ real roots (counting
multiplicity) and the maximum multiplicity of the roots is $m(f)$, then

$g_{s}(u)=f(u)+sf^{\prime}(u)$

also has $N$ real roots (counting multiplicity) and
$m(g_{s})=\max\\{m(f)-1,1\\}$ if $s\neq 0$.

###### Proof.

According to the condition, we can write
$f(u)=C(u-u_{1})^{m_{1}}\cdots(u-u_{r})^{m_{r}}$ with $\sum_{i=1}^{r}m_{i}=N$.
Denote by $u_{i}^{\prime}\in(u_{i},u_{i+1})$ the root of $f^{\prime}(u)$ for
$i=1,\dots,r-1$. Set $u_{0}^{\prime}=-\infty$ and $u_{r}^{\prime}=\infty$. It
suffices to show that, if $s\neq 0$,

$g_{s}(u)=C\prod_{i=1}^{r}(u-u_{i})^{m_{i}-1}\cdot\prod_{i=1}^{r}(u-v_{i})$

with $v_{i}\in(u_{i-1}^{\prime},u_{i}^{\prime})$ and $v_{i}\neq u_{i}$ for
$i=1,\dots,r$.

For this purpose, we first notice that $u_{i}$ is a root of $g_{s}(u)$ with
multiplicity $(m_{i}-1)$. Then we look into the interval
$(u_{i-1}^{\prime},u_{i}^{\prime})$ which contains $u_{i}$ for $i=1,\dots,r$.
Note $g_{s}(u_{i}^{\prime})=f(u_{i}^{\prime})$ for $i=0,...,r$ (this equality
holds for $i=0,r$ where $f(u_{i}^{\prime})=\pm\infty$). If $m_{i}$ is even, we
have $f(u_{i-1}^{\prime})f(u_{i}^{\prime})>0$ while
$g_{s}(u_{i}-\epsilon)g_{s}(u_{i}+\epsilon)<0$ because the multiplicity of
$u_{i}$ is odd ($=m_{i}-1$) for $g_{s}(u)$. Therefore, there exists one root
of $g_{s}(u)$, denoted by $v_{i}$, in either
$(u_{i-1}^{\prime},u_{i}-\epsilon)$ or $(u_{i}+\epsilon,u_{i}^{\prime})$. It
is hence distinct from $u_{i}$. Similarly, such a $v_{i}$ also exists for odd
$m_{i}$. In this way, we get $r$ additional roots $v_{i}$ ($i=1,\dots,r$) of
$g_{s}(u)$. These roots are all simple because the degree of $g_{s}(u)$ is
$\sum_{i=1}^{r}(m_{i}-1)+r=N$. Hence, $g_{s}(u)$ can be factorized as above. ∎

Now we are in a position to prove Theorem 2.8.

###### Proof of Theorem 2.8.

Referring to the condition of Theorem 2.8, we can write the $b$-polynomial as
$p(t)=(t+d_{1})\cdots(t+d_{2n+1})$ with at least two of the $d_{i}$’s being
nonzero. It is straightforward to verify that

$\sum_{j=0}^{2n+1}b_{j}t^{j}=t^{2n+1}p\left(\frac{1}{t}\right)=(1+d_{1}t)\cdots(1+d_{2n+1}t).$

Substituting $t$ with $\sigma\partial_{u}$, we obtain

$\sum_{j=0}^{2n+1}b_{j}\sigma^{j}\partial_{u}^{j}g(u;W)=(1+d_{1}\sigma\partial_{u})\cdots(1+d_{2n+1}\sigma\partial_{u})g(u;W),$

which is the characteristic polynomial $c(u;M)$ according to (4.3).

As shown in Proposition 4.1, $g(u;W)$ has $(2n+1)$ real roots (of $u$) and the
maximum multiplicity $m(g)\leq 3$ if $W\in\Omega$. Then by repeatedly using
Proposition 4.2 with $s=\sigma d_{i}$ ($i=2n+1,2n,\dots,1$), we see that the
characteristic polynomial $c(u;M)$ has $(2n+1)$ real roots. Moreover, since at
least two of the $d_{i}$’s are nonzero, the maximum multiplicity of each root
is reduced to 1. Hence, all the roots are distinct. ∎

###  4.2. A proof of Theorem 2.9

First of all, we recall (4.3) for $n=2$ that the characteristic polynomial is

$c(u;M)=\sum_{k=0}^{5}b_{k}\sigma^{k}\partial_{u}^{k}g(u;W)$

and Proposition 4.1 reads as

$g(u;W)=(u-u_{1})^{2}(u-u_{2})^{2}(u-\tilde{u}(W)).$

Moreover, for even kernels we have
$\mathfrak{m}_{1}=\mathfrak{m}_{3}=\mathfrak{m}_{5}=0$. Then it follows from
the definition (2.10) that

$b_{2}=-\frac{1}{2},\quad b_{4}=\frac{1}{4}-\frac{1}{24}\mathfrak{m}_{4},\quad
b_{1}=b_{3}=b_{5}=0.$ (4.6)

Notice that $\tilde{u}=a_{2n}-2\sum_{i=1}^{n}u_{i}$. We can obtain $a_{4}$ by
solving (4.4) via Crammer’s rule and thereby

$\tilde{u}(W)=\frac{w_{1}u_{1}+w_{2}u_{2}}{w_{1}+w_{2}}+\frac{4(\mathfrak{m}_{4}-3)\sigma^{2}(w_{1}-w_{2})(u_{1}-u_{2})}{(w_{1}+w_{2})[(u_{1}-u_{2})^{2}+2(\mathfrak{m}_{4}-3)\sigma^{2}]}.$
(4.7)

Clearly, $\tilde{u}(W)$ and thereby $g(u;W)$ are well-defined for
$u_{1}=u_{2}$.

With these preparations, we turn to

The proof of Theorem 2.9.  By Corollary 2.6 and Remark 2.7, the condition
$\mathfrak{m}_{4}\geq 3$ ensures that the EQMOM is well defined. Assume that
the two-node EQMOM is strictly hyperbolic on $\mathcal{M}(\Omega^{tot})$.
Taking $W=(\frac{1}{2},0,\frac{1}{2},0,1)\in\Omega^{eq}$, we have
$g(u;W)=u^{5}$ and $c(u;W)=u(u^{4}-10u^{2}+120b_{4})$. The strict
hyperbolicity that $c(u;W)$ has 5 distinct real roots implies $b_{4}>0$ and
hence $\mathfrak{m}_{4}<6$ due to (4.6).

Conversely, if $3\leq\mathfrak{m}_{4}<6$, we see from (4.6) that
$b_{2}=-\frac{1}{2}<0<b_{4}\leq\frac{1}{8}<\frac{5}{6}b_{2}^{2}$. Thus,
$c(u;W)=g(u;W)+b_{2}\sigma^{2}\partial_{u}^{2}g(u;W)+b_{4}\sigma^{4}\partial_{u}^{4}g(u;W)$
has 5 distinct real roots due to Proposition 4.3 to be proved below.

######  Proposition 4.3.

Set $g(u)=g(u;W)$. For any $s>0$, the polynomial

$h(u;s)=g(u)-c_{1}sg^{(2)}(u)+c_{2}s^{2}g^{(4)}(u)$

has 5 distinct roots if the constants satisfy $c_{1}>0$ and
$0<c_{2}<\frac{5}{6}c_{1}^{2}$.

###### Proof.

Without loss of generality, we assume $u_{1}\leq u_{2}$ and consider four
cases: (I) $u_{1}=u_{2}=\tilde{u}$, (II) $u_{1}<\tilde{u}<u_{2}$, (III)
$u_{1}<u_{2}\leq\tilde{u}$ and (IV) $\tilde{u}\leq u_{1}<u_{2}$.

Case I. In this case, we have $g(u)=(u-u_{1})^{5}$ and

$h(u;s)=(u-u_{1})^{5}-20c_{1}s(u-u_{1})^{3}+120c_{2}s^{2}(u-u_{1})=(u-u_{1})(t^{2}-20c_{1}st+120c_{2}s^{2})$

with $t=(u-u_{1})^{2}$. Thanks to $c_{1}>0$ and
$0<c_{2}<\frac{5}{6}c_{1}^{2}$, the quadratic function in $t$ has two distinct
positive roots. Therefore $h(u;s)$ has 5 distinct real roots.

For other cases, notice that $h(u;s)$ is a monic polynomial of order 5. Thus,
it suffices to find $\hat{u}_{1}<\dots<\hat{u}_{4}$ such that

$h(\hat{u}_{1};s)>0\quad\text{and}\quad
h(\hat{u}_{i};s)h(\hat{u}_{i+1};s)<0,\quad i=1,2,3.$

For this purpose, we set $P_{u}^{+}:=\\{s>0|h(u;s)>0\\}$ and
$P_{u}^{-}:=\\{s>0|h(u;s)<0\\}$ for each $u$. Then our main tast is to show

$\bigcup\limits_{u\in I}P_{u}^{+}=(0,\infty)=\bigcup\limits_{u\in
I^{\prime}}P_{u}^{-}$ (4.8)

for some intervals $I,I^{\prime}$. To do this, the key idea is to view
$h(u;s)$ as a quadratic function of $s$ and check its symmetric axis and
discriminant

$A(u)=\frac{c_{1}g^{(2)}(u)}{2c_{2}g^{(4)}(u)},\quad d(u)=g^{(2)}(u)^{2}-\beta
g^{(4)}(u)g(u)$

with $\beta=4c_{2}/c_{1}^{2}\in(0,\frac{10}{3})$. If the leading coefficient
$c_{2}g^{(4)}(u)$ is negative (or positive) and the discriminant $d(u)$ is
positive in $I$ (resp. $I^{\prime}$), then $P_{u}^{+}$ (resp. $P_{u}^{-}$) is
an open interval centered at $s=A(u)$.

In addition, we recall the form of $g(u)$ given in Proposition 4.1 and denote
by $u^{(j)}_{1}\leq\dots\leq u^{(j)}_{5-j}$ the $(5-j)$ roots of $g^{(j)}(u)$
(counting multiplicity). It is not difficult to see that $u^{(j)}_{i}\leq
u^{(j+1)}_{i}\leq u^{(j)}_{i+1}$ for $j=1,2,3$ and $i=1,\dots,4-j$. The
equalities occur only when $j=1$ and $\tilde{u}=u_{1}$ or $u_{2}$.

Case II. In this case, we have
$u_{1}=u_{1}^{(1)}<u^{(2)}_{1}<u^{(1)}_{2}<\tilde{u}$ and
$u^{(2)}_{1}<u^{(3)}_{1}<u^{(4)}_{1}$, resulting in $g(u^{(2)}_{1})<0$ and
$g^{(4)}(u^{(2)}_{1})<0$. Thus, we have $h(u^{(2)}_{1};s)<0$ for any $s>0$. A
similar argument shows $h(u^{(2)}_{3};s)>0$. As a consequence, it suffices to
prove that there exist $\hat{u}_{1}\leq u_{1}$ and $\hat{u}_{4}\geq u_{2}$
such that $h(\hat{u}_{1};s)>0$ and $h(\hat{u}_{4};s)<0$; in other words, we
need to verify (4.8) for $I=(-\infty,u_{1}]$ and $I^{\prime}=[u_{2},\infty)$.

To this end, we notice $A(u)\to+\infty$ as $u\to-\infty$ and
$P_{u_{1}}^{+}=(0,2A(u_{1}))$. Owing to the continuity of $A(u)$, it suffices
to show that $d(u)>0$ for $u\leq u_{1}$. Notice that $d(u)$ is a polynomial of
order 6. A straightforward calculation yields (the $u$-dependence is omitted
for clarity)

$\displaystyle d^{(1)}$
$\displaystyle=2g^{(2)}g^{(3)}-\beta(g^{(5)}g+g^{(4)}g^{(1)}),$ (4.9)
$\displaystyle d^{(2)}$
$\displaystyle=2g^{(3)2}+(2-\beta)g^{(2)}g^{(4)}-2\beta g^{(1)}g^{(5)},$
$\displaystyle d^{(3)}$
$\displaystyle=(6-\beta)g^{(3)}g^{(4)}+(2-3\beta)g^{(2)}g^{(5)},$
$\displaystyle d^{(4)}$
$\displaystyle=(6-\beta)g^{(4)2}+(8-4\beta)g^{(3)}g^{(5)},$ $\displaystyle
d^{(5)}$ $\displaystyle=(20-6\beta)g^{(4)}g^{(5)},$ $\displaystyle d^{(6)}$
$\displaystyle=(20-6\beta)g^{(5)2}.$

Obviously, we have $d(u_{1})>0$, $d^{(1)}(u_{1})<0$, $d^{(5)}(u_{1})<0$ and
$d^{(6)}(u_{1})>0$. Setting $a=u_{2}-u_{1}>0$ and
$t=\frac{\tilde{u}-u_{1}}{u_{2}-u_{1}}\in(0,1)$, we obtain

$\displaystyle d^{(2)}(u_{1})$
$\displaystyle=24a^{4}[(16-2\beta)t^{2}+(20-4\beta)t+3],$ (4.10)
$\displaystyle d^{(3)}(u_{1})$
$\displaystyle=-48a^{3}[6(6-\beta)t^{2}+5(20-6\beta)t+6(6-\beta)],$
$\displaystyle d^{(4)}(u_{1})$
$\displaystyle=576a^{2}[(6-\beta)t^{2}+(44-14\beta)t+(34-9\beta)].$

For $0<\beta<\frac{10}{3}$, the above $t$-parabolae are positive on
$[0,\infty)$. Thus, we have $d^{(2)}(u_{1})>0$, $d^{(3)}(u_{1})<0$ and
$d^{(4)}(u_{1})>0$. The Taylor expansion of $d(u)$ at $u_{1}$ then verifies
the positivity of $d(u)$ for $u\leq u_{1}$.

A similar argument can be used to show $\bigcup\limits_{u\geq
u_{2}}P_{u}^{-}=(0,\infty)$.

Case III. In this case we take $\hat{u}_{2}=u^{(2)}_{1}$ and if $u_{2}\geq
u^{(4)}_{1}$, we take $\hat{u}_{3}=u_{2}$. It is then easy to see that
$h(\hat{u}_{2};s)<0<h(\hat{u}_{3};s)$ for any $s>0$.

For $u_{2}<u^{(4)}_{1}$, we can show the existence of
$\hat{u}_{3}\in(u_{2},u^{(4)}_{1})$ such that $h(\hat{u}_{3},s)>0$ by
verifying the first equality in (4.8) with $I=(u_{2},u_{1}^{(4)})$. Clearly,
for $u\in I$ we have $c_{2}g^{(4)}(u)<0$, $P_{u_{2}}^{+}=(0,2A(u_{2}))$ and
$A(u)\to\infty$ as $u\to u^{(4)}_{1}$. Then it suffice to show $d(u)>0$ for
$u\in(u_{2},u_{1}^{(4)})$. Note that $g^{(j)}(u)<0$ for
$u\in(u_{2},u^{(4)}_{1})$ because $u_{2}=u^{(1)}_{3}>u^{(2)}_{2}>u^{(3)}_{1}$
and $u^{(4)}_{1}<u^{(3)}_{2}<u^{(1)}_{4}<\tilde{u}$. Then we see that
$d(u_{2})>0$ and $d^{(1)}(u_{2})>0$. Moreover, if $\beta\leq 2$, we see from
(4.9) that $d^{(2)}(u)>0$ throughout $(u_{2},u^{(4)}_{1})$. Consequently, we
have $d(u)>0$ for $u\in(u_{2},u^{(4)}_{1})$. On the other hand, if $\beta>2$,
it follows from (4.9) that $d^{(3)}(u)>0$ for $u\in(u_{2},u^{(4)}_{1})$. As
for $d^{(2)}(u_{2})$, it has the same form as in (4.10) except that $u_{1}$
and $u_{2}$ are permuted. Clearly, we have $d^{(2)}(u_{2})>0$. In conclusion,
we have shown that $d(u)>0$ for $u\in(u_{2},u^{(4)}_{1})$ in both cases.

To show the existence of $\hat{u}_{4}$, we verify the second equality in (4.8)
with $I^{\prime}=(\tilde{u},\infty)$. As above, it suffices to show $d(u)>0$
for $u>\tilde{u}$. From (4.9), we see that $d(\tilde{u})>0$,
$d^{(5)}(\tilde{u})>0$ and $d^{(6)}(\tilde{u})>0$. Furthermore, we set
$a=\tilde{u}-u_{1}>0$ and $t=\frac{\tilde{u}-u_{2}}{\tilde{u}-u_{1}}\in[0,1)$.
A straightforward calculation gives

$g^{(1)}(\tilde{u})=a^{4}p_{1}(t),\ g^{(2)}(\tilde{u})=4a^{3}p_{2}(t),\
g^{(3)}(\tilde{u})=6a^{2}p_{3}(t),\ g^{(4)}(\tilde{u})=48ap_{4}(t)$

with $p_{1}(t)=t^{2}$, $p_{2}(t)=t^{2}+t$, $p_{3}(t)=t^{2}+4t+1$ and
$p_{4}(t)=t+1$. Owing to the following inequalities

$\displaystyle p_{2}(t)p_{3}(t)$
$\displaystyle=(t^{2}+t)(t^{2}+4t+1)\geq(t^{2}+t)\cdot 6t=6p_{1}(t)p_{4}(t),$
$\displaystyle p_{3}(t)^{2}$ $\displaystyle=((t+1)^{2}+2t)^{2}\geq
8t(t+1)^{2}=8p_{2}(t)p_{4}(t),$ $\displaystyle p_{3}(t)^{2}$
$\displaystyle=(t^{2}+4t+1)^{2}\geq 36t^{2}=36p_{1}(t),$ $\displaystyle
p_{3}(t)p_{4}(t)$ $\displaystyle=(t^{2}+4t+1)(t+1)\geq 6t(t+1)=6p_{2}(t),$
$\displaystyle 3p_{4}(t)^{2}$ $\displaystyle=3(t+1)^{2}\geq
4t+2(t+1)^{2}=2p_{3}(t),$

we derive, by a direct calculation, that

$\displaystyle d^{(1)}(\tilde{u})$
$\displaystyle=48a^{5}\left(p_{2}(t)p_{3}(t)-\beta p_{4}(t)p_{1}(t)\right)\geq
48a^{5}(6-\beta)p_{1}(t)p_{4}(t)>0,$ $\displaystyle d^{(2)}(\tilde{u})$
$\displaystyle=72a^{4}p_{3}(t)^{2}-192(\beta-2)a^{4}p_{2}(t)p_{4}(t)-240\beta
a^{4}p_{1}(t)$
$\displaystyle\geq(256-192(\beta-2))a^{4}p_{2}(t)p_{4}(t)+240(6-\beta)a^{4}p_{1}(t)>0,$
$\displaystyle d^{(3)}(\tilde{u})$
$\displaystyle=288(6-\beta)a^{3}p_{3}(t)p_{4}(t)-480(3\beta-2)a^{3}p_{2}(t)\geq
96(118-33\beta)a^{3}p_{2}(t)>0,$ $\displaystyle d^{(4)}(\tilde{u})$
$\displaystyle=48^{2}(6-\beta)a^{2}p_{4}(t)^{2}-720(4\beta-8)a^{2}p_{3}(t)\geq
192(78-23\beta)a^{2}p_{3}(t)>0.$

These ensure that $d(u)>0$ and thereby the existence of $\hat{u}_{4}$.

A similar argument as in Case (II) can be used to prove the existence of
$\hat{u}_{1}\leq u_{1}$ such that $h(\hat{u}_{1};s)>0$ for any $s>0$.

Case IV. This case can be converted to Case (III) with
$-u_{2}<-u_{1}\leq-\tilde{u}$ by introducing
$\tilde{g}(u)=-g(-u)=(u+u_{1})^{2}(u+u_{2})^{2}(u+\tilde{u})$. We see that
$-h(-u;s)=\tilde{g}(u)-c_{1}s\tilde{g}^{(2)}(u)+c_{2}s^{2}\tilde{g}^{(4)}(u)$
has 5 distinct roots, so is $h(u;s)$. Hence, the proof is complete. ∎

##  5\. Structural stability

In this section we prove Theorem 2.10, namely, checking the structural
stability condition (I)–(III) in Subsection 2.2.

For Condition (I), we calculate the Jacobian of
$S(M)=\frac{1}{\tau}\left(\rho\Delta^{eq}(U,\sqrt{\theta})-M\right)$ as

$S_{M}(M)=\frac{1}{\tau}\begin{bmatrix}0&&&&\\\ &0&&&\\\ &&0&&\\\
s_{1}&s_{2}&s_{3}&-1&\\\ s_{4}&s_{5}&s_{6}&&-1\end{bmatrix}$ (5.1)

with

$\displaystyle s_{1}=U^{3}-3U\theta,$ $\displaystyle
s_{2}=3(\theta-U^{2}),\quad s_{3}=3U,$ $\displaystyle
s_{4}=3(U^{4}-2U^{2}\theta-\theta^{2}),$ $\displaystyle s_{5}=-8U^{3},\quad
s_{6}=6(U^{2}+\theta).$

Take

$P(M)=\begin{bmatrix}1&&&&\\\ &1&&&\\\ &&1&&\\\ -s_{1}&-s_{2}&-s_{3}&1&\\\
-s_{4}&-s_{5}&-s_{6}&&1\end{bmatrix},$ (5.2)

It is obvious that $P(M)S_{M}(M)=\tau^{-1}\text{diag}(0,0,0,-1,-1)P(M)$. This
justifies Condition (I). Note that the choice of $P$ is unique up to a block-
diagonal matrix.

As to Condition (II), we know from Theorem 2.9 that the two-node moment system
(2.7) with even kernels is strictly hyperbolic if $3\leq\mathfrak{m}_{4}<6$.
Namely, the coefficient matrix $A(M)$ has 5 distinct real eigenvalues
$\lambda_{i}=\lambda_{i}(M)(i=1,...,5)$. Corresponding to these eigenvalues,
the left eigenvectors form the following matrix

$L=L(M)=\begin{bmatrix}\lambda_{1}^{4}&\lambda_{1}^{3}&\lambda_{1}^{2}&\lambda_{1}&1\\\
\lambda_{2}^{4}&\lambda_{2}^{3}&\lambda_{2}^{2}&\lambda_{2}&1\\\
\lambda_{3}^{4}&\lambda_{3}^{3}&\lambda_{3}^{2}&\lambda_{3}&1\\\
\lambda_{4}^{4}&\lambda_{4}^{3}&\lambda_{4}^{2}&\lambda_{4}&1\\\
\lambda_{5}^{4}&\lambda_{5}^{3}&\lambda_{5}^{2}&\lambda_{5}&1\end{bmatrix}\begin{bmatrix}-1&&&&\\\
a_{4}&-1&&&\\\ a_{3}&a_{4}&-1&&\\\ a_{2}&a_{3}&a_{4}&-1&\\\
a_{1}&a_{2}&a_{3}&a_{4}&-1\end{bmatrix},$ (5.3)

which can be easily verified. With this matrix, the symmetrizer
$A_{0}=A_{0}(M)$ in Condition (II) must be chosen as $A_{0}=L^{T}\Lambda L$
with $\Lambda$ an arbitrary positive definite diagonal matrix to be determined
[25].

The rest of this section is to choose the diagonal matrix $\Lambda=\Lambda(M)$
such that Condition (III) is satisfied for $M$ in the equilibrium manifold.
Since $P(M)S_{M}(M)=\tau^{-1}\text{diag}(\mathbf{0}_{3},-I_{2})P(M)$, it is
equivalent to find $\Lambda$ such that the matrix
$P^{-T}A_{0}P^{-1}=P^{-T}L^{T}\Lambda LP^{-1}$ is block diagonal with the same
partition as $\text{diag}(\mathbf{0}_{3},-I_{2})$, meaning that the first
three columns of $\sqrt{\Lambda}LP^{-1}$ are orthogonal to the last two
columns. Note that the existence of such a $\Lambda$ is independent of the
choice of $P$. Denote by $r_{i}\in\mathbb{R}^{5}$ the $i$th column of
$\sqrt{\Lambda}L$. Since

$P^{-1}=\begin{bmatrix}1&&&&\\\ &1&&&\\\ &&1&&\\\ s_{1}&s_{2}&s_{3}&1&\\\
s_{4}&s_{5}&s_{6}&&1\end{bmatrix},$

the orthogonality gives six equations

$(r_{i},r_{j})+s_{i}(r_{4},r_{j})+s_{i+3}(r_{5},r_{j})=0\quad\text{for
}i=1,2,3\text{ and }j=4,5.$ (5.4)

Here $(\cdot,\cdot)$ represents the dot product of vectors.

To show that (5.4) can be used to determine $\Lambda$, we write
$\sqrt{\Lambda}=\text{diag}(x_{i})_{i=1}^{5}$. Then it follows from (5.3) that

$r_{i}=\left(x_{1}\sum_{j=i}^{5}a_{j}\lambda_{1}^{j-i},\ \dots,\
x_{5}\sum_{j=i}^{5}a_{j}\lambda_{5}^{j-i}\right)^{T}$ (5.5)

with $a_{5}=-1$. With this, the dot products can be written as

$(r_{k},r_{5})=-\sum_{i=1}^{5}x_{i}^{2}\sum_{j=k}^{5}a_{j}\lambda_{i}^{j-k}.$

Moreover, we introduce

$r_{4}^{\prime}=r_{4}+a_{4}r_{5}=-(x_{1}\lambda_{1},\dots,x_{5}\lambda_{5})^{T}.$

It is clear that the equations in (5.4) with $j=4$ can be replaced with

$(r_{i},r_{4}^{\prime})+s_{i}(r_{4},r_{4}^{\prime})+s_{i+3}(r_{5},r_{4}^{\prime})=0,\quad
i=1,2,3,$

and

$\begin{split}(r_{k},r_{4}^{\prime})&=\sum_{i=1}^{5}x_{i}^{2}\left[\lambda_{i}\left(\lambda_{i}^{5-k}-a_{4}\lambda_{i}^{4-k}-\cdots-
a_{k}\right)-a_{k-1}+a_{k-1}\right]\\\
&=(r_{k-1},r_{5})+a_{k-1}\sum_{i=1}^{5}x_{i}^{2}.\end{split}$ (5.6)

These indicate that (5.4) is a system of six linear equations for the five
unknowns $x_{i}^{2}$. Note that the coefficients of this system all depend on
$M$.

Since Condition (III) is posed only for $M$ in the equilibrium, we only need
to calculate the coefficients for $M$ in the equilibrium manifold for the
moment system (2.7):

$\mathcal{E}=\\{M\in\mathcal{M}(\Omega^{tot})|M_{j}=\rho\Delta_{j}^{eq}(U,\sqrt{\theta})\text{
for }j=0,\dots,4\\}$

with $\Delta_{j}^{eq}(U,\sqrt{\theta})$ defined in (2.2) and

$\rho=M_{0},\quad
U=M_{1}/M_{0}\quad\text{and}\quad\theta=(M_{0}M_{2}-M_{1}^{2})/M_{0}^{2}.$

About this $\mathcal{E}$, we have

######  Proposition 5.1.

For any $\rho,\theta>0$ and $U\in\mathbb{R}$, consider equations

$M_{j}:=\sum_{i=1}^{2}w_{i}\Delta_{j}(u_{i},\sigma)=\rho\Delta_{j}^{eq}(U,\sqrt{\theta}),\quad
j=0,1,...,4,$

for $W=(w_{1},u_{1},w_{2},u_{2},\sigma)\in\Omega^{tot}$. When
$\mathfrak{m}_{4}<3$, the equations have no solution; when
$\mathfrak{m}_{4}=3$, the solutions satisfy
$u_{1}=u_{2}=U,\sigma=\sqrt{\theta}$ and $w_{1}+w_{2}=\rho$; and when
$\mathfrak{m}_{4}>3$, there is a unique solution given as

$w_{1,2}=\frac{\rho}{2},\quad
u_{1,2}=U\mp\nu\sigma,\quad\sigma=\left(\frac{\theta}{\nu^{2}+1}\right)^{\frac{1}{2}},\quad\nu=\left(\frac{\mathfrak{m}_{4}-3}{2}\right)^{\frac{1}{4}}.$

In particular, the equilibrium manifold is nonempty if and only if
$\mathfrak{m}_{4}\geq 3$.

###### Proof.

Notice that

$\displaystyle\sum_{i=1}^{2}w_{i}\Delta_{j}\left(\frac{u_{i}-U}{\sigma},1\right)=$
$\displaystyle\int\xi^{j}\sum_{i=1}^{2}w_{i}\mathcal{K}\left(\xi-\frac{u_{i}-U}{\sigma}\right)d\xi$
$\displaystyle=$
$\displaystyle\int\left(\frac{\eta-U}{\sigma}\right)^{j}\sum_{i=1}^{2}\frac{w_{i}}{\sigma}\mathcal{K}\left(\frac{\eta-
u_{i}}{\sigma}\right)d\eta=\sum_{k=0}^{j}\binom{j}{k}(-U)^{k}\frac{M_{j-k}}{\sigma^{j}}$

and similarly,

$\rho\Delta_{j}^{eq}\left(0,\frac{\sqrt{\theta}}{\sigma}\right)=\sum_{k=0}^{j}\binom{j}{k}(-U)^{k}\frac{\rho\Delta_{j-k}^{eq}(U,\sqrt{\theta})}{\sigma^{j}}.$

Then the given equations $M_{j}=\rho\Delta_{j}^{eq}(U,\sqrt{\theta})$ are
equivalent to

$\sum_{i=1}^{2}\frac{w_{i}}{\rho}\Delta_{j}\left(\frac{u_{i}-U}{\sigma},1\right)=\Delta_{j}^{eq}\left(0,\frac{\sqrt{\theta}}{\sigma}\right),\quad
j=0,...,4.$ (5.7)

Denote $w_{i}^{\prime}=\frac{w_{i}}{\rho}$,
$u_{i}^{\prime}=\frac{u_{i}-U}{\sigma}$ and
$\theta^{\prime}=\frac{\theta}{\sigma^{2}}$. Recall from (3.1) that

$\Delta_{j}(u,\sigma)=\sum_{k=0}^{j}\binom{j}{k}\mathfrak{m}_{k}\sigma^{k}u^{j-k}$

with $\mathfrak{m}_{3}=0$ for even kernels. Then equations (5.7) can be
rewritten as

$\displaystyle w_{1}^{\prime}+w_{2}^{\prime}$ $\displaystyle=1,$ (5.8a)
$\displaystyle
w_{1}^{\prime}u_{1}^{\prime}+w_{2}^{\prime}u_{2}^{\prime}=w_{1}^{\prime}u_{1}^{\prime
3}+w_{2}^{\prime}u_{2}^{\prime 3}$ $\displaystyle=0,$ (5.8b) $\displaystyle
w_{1}^{\prime}u_{1}^{\prime 2}+w_{2}^{\prime}u_{2}^{\prime 2}$
$\displaystyle=\theta^{\prime}-1,$ (5.8c) $\displaystyle
w_{1}^{\prime}u_{1}^{\prime 4}+w_{2}^{\prime}u_{2}^{\prime 4}$
$\displaystyle=3\theta^{\prime 2}-6(\theta^{\prime}-1)-\mathfrak{m}_{4}.$
(5.8d)

Multiplying both sides of (5.8c) with $(u_{1}^{\prime}+u_{2}^{\prime})$ and
using (5.8b), we obtain

$(\theta^{\prime}-1)(u_{1}^{\prime}+u_{2}^{\prime})=0.$

This gives $u_{1}^{\prime}=-u_{2}^{\prime}:=\nu\neq 0$ if $\theta^{\prime}\neq
1$. Then we see from (5.8a) and (5.8b) that
$w_{1}^{\prime}=w_{2}^{\prime}=\frac{1}{2}$. Thus, (5.8d) becomes
$2\nu^{4}=\mathfrak{m}_{4}-3$ with $\theta^{\prime}=1+\nu^{2}$. When
$\theta^{\prime}=1$, it follows from (5.8c) that
$u_{1}^{\prime}=u_{2}^{\prime}=0$ due to $w_{i}^{\prime}>0$ and
$\mathfrak{m}_{4}=3$ from (5.8d). Hence, the given equations have a solution
if and only if $\mathfrak{m}_{4}\geq 3$. ∎

At equilibrium, it is seen from Proposition 5.1 and (4.7) that $\tilde{u}=U$
and

$g(u;W)=(u-U-\sigma\nu)^{2}(u-U+\sigma\nu)^{2}(u-U).$

We then rewrite $g(u;W)$ and

$c(u;W)=g(u;W)+b_{2}\sigma^{2}g^{(2)}(u;W)+b_{4}\sigma^{4}g^{(4)}(u;W),$

with $b_{j}$ in (4.6) as $g(u;U,\sigma)$ and $c(u;U,\sigma)$ to manifest the
dependence on $U$ and $\sigma$. Clearly, we have
$g(u;U,\sigma)=\sigma^{5}g\left(\frac{u-U}{\sigma};0,1\right)$ and
$c(u;U,\sigma)=\sigma^{5}c\left(\frac{u-U}{\sigma};0,1\right)$.

A direct calculation yields $c(u;0,1)=u^{5}-B_{1}u^{3}+B_{2}u$ with

$B_{1}=2\nu^{2}+10\quad\text{and }B_{2}=-9\nu^{4}+6\nu^{2}+15.$ (5.9)

Recalling $c(u;W)=u^{5}-a_{4}u^{4}-\cdots-a_{0}$, we derive

$\begin{split}a_{0}&=U^{5}-B_{1}U^{3}\sigma^{2}+B_{2}U\sigma^{4},\
a_{1}=-5U^{4}+3B_{1}U^{2}\sigma^{2}-B_{2}\sigma^{4},\\\
a_{2}&=10U^{3}-3B_{1}U\sigma^{2},\quad a_{3}=-10U^{2}+B_{1}\sigma^{2},\quad
a_{4}=5U.\end{split}$ (5.10)

This determines the matrix $L=L(M)$ in (5.3) on the equilibrium manifold.

Denote by $\pm\mu_{1}$, $\pm\mu_{2}$ and 0 the five distinct roots of the
polynomial $c(u;0,1)$. Then the roots $\lambda_{i}$ of $c(u;U,\sigma)$ can be
written as

$\lambda_{1,2}=U\pm\mu_{1}\sigma,\quad\lambda_{3,4}=U\pm\mu_{2}\sigma,\quad\lambda_{5}=U.$
(5.11)

By introducing $V=U/\sigma$, it is seen that $r_{i}$ is a homogeneous
polynomial of $\sigma$ (of degree $5-i$), so the relations in (5.4) are all
homogeneous with $\sigma$. Thus, we may set $\sigma=1$ for the following
calculations.

Set

$Y_{0}=\sum_{i=1}^{5}x_{i}^{2},\quad
Y_{j}=\left(x_{1}^{2}+(-1)^{j}x_{2}^{2}\right)\mu_{1}^{j}+\left(x_{3}^{2}+(-1)^{j}x_{4}^{2}\right)\mu_{2}^{j},\quad\text{for
}j=1,\dots,4.$

Having (5.10) & (5.11), the dot products $(r_{i},r_{5})$ can be calculated as

$\displaystyle(r_{1},r_{5})$
$\displaystyle=(U^{4}-B_{1}U^{2}+B_{2})Y_{0}-(U^{3}-B_{1}U)Y_{1}+(U^{2}-B_{1})Y_{2}-UY_{3}+Y_{4},$
$\displaystyle(r_{2},r_{5})$
$\displaystyle=(-4U^{3}+2B_{1}U)Y_{0}+(3U^{2}-B_{1})Y_{1}-2UY_{2}+Y_{3},$
$\displaystyle(r_{3},r_{5})$
$\displaystyle=(6U^{2}-B_{1})Y_{0}-3UY_{1}+Y_{2},$
$\displaystyle(r_{4},r_{5})$ $\displaystyle=-4UY_{0}+Y_{1},$
$\displaystyle(r_{5},r_{5})$ $\displaystyle=Y_{0}.$

Furthermore, we use (5.6) to get

$\displaystyle(r_{1},r_{4}^{\prime})$
$\displaystyle=a_{0}Y_{0}=(U^{5}-B_{1}U^{3}+B_{2}U)Y_{0},$
$\displaystyle(r_{2},r_{4}^{\prime})$
$\displaystyle=(2B_{1}U^{2}-4U^{4})Y_{0}-(U^{3}-B_{1}U)Y_{1}+(U^{2}-B_{1})Y_{2}-UY_{3}+Y_{4},$
$\displaystyle(r_{3},r_{4}^{\prime})$
$\displaystyle=(6U^{3}-B_{1}U)Y_{0}+(3U^{2}-B_{1})Y_{1}-2UY_{2}+Y_{3},$
$\displaystyle(r_{4},r_{4}^{\prime})$
$\displaystyle=-4U^{2}Y_{0}-3UY_{1}+Y_{2},$
$\displaystyle(r_{5},r_{4}^{\prime})$ $\displaystyle=UY_{0}+Y_{1}.$

With these and those in (5.1), the equations in (5.4) can be written as

$\displaystyle 0=$ $\displaystyle\left[(6\theta-
B_{1})U^{2}+B_{2}-3\theta^{2}\right]Y_{0}+(B_{1}-3\theta)UY_{1}$ (5.12a)
$\displaystyle+(U^{2}-B_{1})Y_{2}-UY_{3}+Y_{4},$ $\displaystyle 0=$
$\displaystyle(2B_{1}-12\theta)UY_{0}+(3\theta-B_{1})Y_{1}-2UY_{2}+Y_{3},$
(5.12b) $\displaystyle 0=$ $\displaystyle(6\theta-B_{1})Y_{0}+Y_{2},$ (5.12c)
$\displaystyle 0=$ $\displaystyle\left[(6\theta-
B_{1})U^{3}+(B_{2}-3\theta^{2})U\right]Y_{0}+(3U^{2}\theta-3\theta^{2})Y_{1}+(U^{3}-3U\theta)Y_{2},$
(5.12d) $\displaystyle 0=$
$\displaystyle(2B_{1}-12\theta)U^{2}Y_{0}+(B_{1}-9\theta)UY_{1}-(B_{1}-3\theta+2U^{2})Y_{2}-UY_{3}+Y_{4},$
(5.12e) $\displaystyle 0=$ $\displaystyle(6\theta-B_{1})UY_{0}+(6\theta-
B_{1})Y_{1}+UY_{2}+Y_{3}.$ (5.12f)

We shall show that there exists $x_{i}^{2}$ solving these equations if and
only if $\mathfrak{m}_{4}<5$. To do this, we use (5.12c) and deduce from
(5.12b) & (5.12f) that $Y_{3}=(B_{1}-3\theta)Y_{1}=(B_{1}-6\theta)Y_{1}$ and
hence $Y_{1}=Y_{3}=0$. By the definitions of $Y_{1}$ and $Y_{3}$, it follows
that $x_{1}^{2}=x_{2}^{2}$ and $x_{3}^{2}=x_{4}^{2}$ due to
$\mu_{1}\neq\mu_{2}$. On the other hand, notice that $\theta=\nu^{2}+1$ due to
Proposition 5.1 with $\sigma=1$. We deduce from (5.9) that

$B_{2}-3\theta^{2}=3\theta(B_{1}-6\theta)=-12(\nu^{4}-1)=-6(\mathfrak{m}_{4}-5).$
(5.13)

Thus, it is not difficult to see that all the equations in (5.12) are linear
combinations of (5.12a) and (5.12c).

Furthermore, since $\mu_{1}$ and $\mu_{2}$ are the nonzero roots of
$c(u;0,1)=u^{5}-B_{1}u^{3}+B_{2}u$, we have $\mu_{1}^{2}+\mu_{2}^{2}=B_{1}>0$
and $\mu_{1}^{2}\mu_{2}^{2}=B_{2}>0$. Thus, by the definitions of $Y_{2}$ and
$Y_{4}$, we have $B_{1}Y_{2}=Y_{4}+B_{2}(x_{1}^{2}+\cdots+x_{4}^{2})$.
Consequently, (5.12a) and (5.12c) are equivalent to (using
$x_{1}^{2}=x_{2}^{2}$ and $x_{3}^{2}=x_{4}^{2}$)

$\displaystyle(B_{2}-3\theta^{2})Y_{0}$ $\displaystyle=B_{2}\cdot
2(x_{1}^{2}+x_{3}^{2}),$ (5.14a) $\displaystyle Y_{2}$
$\displaystyle=(B_{1}-6\theta)Y_{0}.$ (5.14b)

Therefore, if $\mathfrak{m}_{4}\geq 5$, from (5.13) and (5.14b) we see that
$Y_{0}$ and $Y_{2}$ cannot be both positive. This together with the
definitions of $Y_{0}$ and $Y_{2}$ indicates the nonexistence of the
$x_{i}^{2}$’s and thereby the diagonal positive matrix $\Lambda$.

Finally, we show the existence of the $x_{i}^{2}$’s if $\mathfrak{m}_{4}<5$.
From (5.13) we see that $\mathfrak{m}_{4}<5$ if and only if $0\leq\nu<1$.
Substituting (5.13) into (5.14a) and using $Y_{0}=\sum_{i=1}^{5}x_{i}^{2}$, we
easily obtain

$x_{1}^{2}+x_{3}^{2}=2\frac{1-\nu^{2}}{\theta}x_{5}^{2}.$

Substituting this into (5.14b) and putting them into matrix form, we arrive at

$\begin{bmatrix}1&1\\\
\mu_{1}^{2}&\mu_{2}^{2}\end{bmatrix}\begin{bmatrix}x_{1}^{2}\\\
x_{3}^{2}\end{bmatrix}=\frac{2(1-\nu^{2})}{\theta}x_{5}^{2}\begin{bmatrix}1\\\
5-3\nu^{2}\end{bmatrix}.$

Assume $\mu_{1}<\mu_{2}$. It is not difficult to verify that the two
components $x_{1}^{2}$ and $x_{3}^{2}$ of the solution to this system are
positive if and only if $\mu_{1}^{2}<5-3\nu^{2}<\mu_{2}^{2}$. Notice that

$\mu_{2}^{2}>\frac{\mu_{1}^{2}+\mu_{2}^{2}}{2}=\frac{B_{1}}{2}=\nu^{2}+5\geq
5-3\nu^{2}\quad\text{and}$
$\mu_{1}^{2}=\frac{B_{1}}{2}-\sqrt{\frac{1}{4}B_{1}^{2}-B_{2}}=5+\nu^{2}-\sqrt{10\nu^{4}+4\nu^{2}+10}<5-3\nu^{2}$

for $0\leq\nu<1$. Hence, the existence of positive $x_{1}^{2}$ and $x_{3}^{2}$
is demonstrated and the proof is complete.

##  6\. Specific kernel functions

In this section we present a number of specific kernel functions which satisfy
the conditions required by our previous theoretical results, and then examine
their performance with numerical tests.

###  6.1. Examples of kernel functions

Here are examples of the kernel functions, which may not be normalized. To
check the conditions in our main results, we notice that, for even normalized
kernels,

$\mathfrak{m}_{4}=\int\sqrt{\tilde{\mathfrak{m}}_{2}}\mathcal{K}(\sqrt{\tilde{\mathfrak{m}}_{2}}\xi)\xi^{4}d\xi=\tilde{\mathfrak{m}}_{2}^{-2}\int\mathcal{K}(\eta)\eta^{4}d\eta=\tilde{\mathfrak{m}}_{4}/\tilde{\mathfrak{m}}_{2}^{2},$

where $\tilde{\mathfrak{m}}_{4},\tilde{\mathfrak{m}}_{2}$ are the $4^{th}$ and
$2^{nd}$ moments of the unnormalized kernel, respectively.

######  Example 6.1 (Gaussian distribution).

Our first example is the most widely-used Gaussian distribution

$\mathcal{K}(\xi)=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{\xi^{2}}{2}\right)$

in the EQMOM approach for the BGK equation [6]. It is even and normalized with
$\mathfrak{m}_{4}=3$. According to Theorem 2.10, the two-node Gaussian-EQMOM
satisfies the structural stability condition. Indeed, for this kernel, it has
been even shown in [13] that the structural stability condition is fulfilled
for the $n$-node EQMOM.

######  Example 6.2 (Piecewise polynomials).

Our next example is the piecewise polynomials

$\mathcal{K}_{j}(\xi)=\frac{j+1}{2}(1-|\xi|)^{j}\mathbf{1}_{|\xi|\leq 1},\quad
j=0,1,\dots.$

These are even functions and $\mathcal{K}_{1}$ was used in [7]. A
straightforward calculation of the unnormalized $k$th-moment
$\tilde{\mathfrak{m}}_{j,k}$ gives

$\frac{\tilde{\mathfrak{m}}_{j,4}}{\tilde{\mathfrak{m}}_{j,2}^{2}}=\frac{6(j+2)(j+3)}{(j+4)(j+5)}$

which is in $[3,6)$ for $j\geq 3$ and in $[3,5)$ for $3\leq j\leq 18$.
According to Theorems 2.9 & 2.10, the corresponding two-node EQMOM moment
system is well-defined and strictly hyperbolic if $j\geq 3$, and satisfies the
structural stability condition if $3\leq j\leq 18$.

On the other hand, even polynomials

$\tilde{\mathcal{K}}_{j}(\xi)=\frac{j+1}{2j}(1-|\xi|^{j})\mathbf{1}_{|\xi|\leq
1},\quad j=1,2,\dots$

cannot be used as EQMOM kernels since
$\tilde{\mathfrak{m}}_{j,4}/\tilde{\mathfrak{m}}_{j,2}^{2}=\frac{9(j+3)^{2}}{5(j+1)(j+5)}<3$,
violating the conditions of Corollary 2.6.

######  Example 6.3 (Kappa distribution).

For kernels to be the kappa distribution [19]

$\mathcal{K}(\xi)=\frac{\Gamma(\kappa+1)}{\kappa\sqrt{\pi\left(\kappa-\frac{3}{2}\right)}\Gamma\left(\kappa-\frac{1}{2}\right)}\left(\frac{\xi^{2}}{\kappa-\frac{3}{2}}+1\right)^{-\kappa}$

with $\kappa>3$, where $\Gamma(z)=\int_{0}^{\infty}t^{z-1}e^{-t}dt$ is the
gamma function, we claim that the two-node EQMOM map $\mathcal{M}$ in (2.5) is
injective. The resultant moment system is strictly hyperbolic if
$\kappa>\frac{7}{2}$ and satisfies the structural stability condition if
$\kappa>4$.

To see this, we only need to consider the rescaled kernel
$\mathcal{K}(\xi)=I_{0,\kappa}^{-1}(\xi^{2}+1)^{-\kappa}$ with
$I_{j,\kappa}=\int_{\mathbb{R}}\xi^{2j}(\xi^{2}+1)^{-\kappa}d\xi$. Obviously
we have

$I_{j,\kappa}=\int_{\mathbb{R}}(\xi^{2j}-\xi^{2j-2}+\xi^{2j-2})(\xi^{2}+1)^{-\kappa}d\xi=I_{j-1,\kappa-1}+I_{j-1,\kappa}.$

Moreover, integrating by parts gives

$I_{j,\kappa}=\frac{-1}{2(\kappa-1)}\int_{\mathbb{R}}\xi^{2j-1}d(\xi^{2}+1)^{1-\kappa}=\frac{2j-1}{2(\kappa-1)}I_{j-1,\kappa-1}$

for $\kappa>j+1/2$. Thus we obtain
$I_{j,\kappa}=\frac{2j-1}{2j-2\kappa+1}I_{j-1,\kappa}$ and

$\frac{\tilde{\mathfrak{m}}_{4}}{\tilde{\mathfrak{m}}_{2}^{2}}=\frac{I_{2,\kappa}/I_{0,\kappa}}{\left(I_{1,\kappa}/I_{0,\kappa}\right)^{2}}=\frac{I_{2,\kappa}/I_{1,\kappa}}{I_{1,\kappa}/I_{0,\kappa}}=\frac{3(2\kappa-3)}{2\kappa-5}>3.$

Therefore, by Corollary 2.6 the two-node EQMOM is well defined. Moreover, we
have $\tilde{\mathfrak{m}}_{4}/\tilde{\mathfrak{m}}_{2}^{2}<6$ and
$\tilde{\mathfrak{m}}_{4}/\tilde{\mathfrak{m}}_{2}^{2}<5$ if
$\kappa>\frac{7}{2}$ and $\kappa>4$, respectively. Consequently, the
conditions of Theorems 2.9 & 2.10 are satisfied if $\kappa>\frac{7}{2}$ or
$\kappa>4$.

######  Example 6.4.

Let

$\mathcal{K}_{j}(\xi)=\frac{1}{2\Gamma(j+1)}|\xi|^{j}e^{-|\xi|},\quad
j=1,2,\dots.$

A straightforward calculation gives
$\tilde{\mathfrak{m}}_{4,j}/\tilde{\mathfrak{m}}_{2,j}^{2}=\frac{(j+3)(j+4)}{(j+1)(j+2)}$,
which is greater than 3 if and only if $j=1$. For $\mathcal{K}_{1}(\xi)$, the
ratio is $\frac{10}{3}<5$. According to Theorem 2.10, the resultant two-node
EQMOM moment system satisfies the structural stability condition.

######  Example 6.5.

Let

$\mathcal{K}(\xi)=\frac{1}{3}(1+\cos\xi)e^{-|\xi|}.$

A direct calculation leads to
$\tilde{\mathfrak{m}}_{4}/\tilde{\mathfrak{m}}_{2}^{2}=14$. According to
Corollary 2.6 and Theorem 2.10, the two-node EQMOM is well-defined, but the
resultant moment system is not hyperbolic.

######  Example 6.6 (Uneven kernels).

Let us give an example of uneven kernels which allows a well-defined and
strictly hyperbolic two-node EQMOM moment system. This kernel function reads
as

$\mathcal{K}(\xi)=\frac{1}{c}\left(\frac{1-\xi}{\xi^{\alpha}}\mathbf{1}_{0<\xi\leq
1}+\frac{1+\xi}{|\xi|^{\beta}}\mathbf{1}_{-1\leq\xi<0}\right)$

with $\alpha=0.6060$, $\beta=0.5340$ and
$c=\frac{1}{(1-\alpha)(2-\alpha)}+\frac{1}{(1-\beta)(2-\beta)}=3.2845$. The
normalized moments $\mathfrak{m}_{j}$ of $\mathcal{K}(\xi)$ are
$\mathfrak{m}_{3}=-0.04006$ and $\mathfrak{m}_{4}=4.7455$. Thus, we have
$\mathfrak{m}_{4}>3+\frac{9}{8}\mathfrak{m}_{3}^{2}$ and the corresponding
two-node EQMOM is well-defined. Its strict hyperbolicity is ensured by Theorem
2.8 as the polynomial

$p(t)=t^{5}-0.5t^{3}+6.6766\times 10^{-3}t^{2}+0.052271t-2.7789\times 10^{-3}$

can be shown to have 5 distinct real roots (this can also be verified
numerically).

######  Example 6.7.

Our last kernel function is

$\mathcal{K}(\xi)=\xi e^{-\xi}\mathbf{1}_{\xi>0}$

and its unnormalized moments are obviously $\tilde{\mathfrak{m}}_{j}=(j+1)!$.
It is not difficult to verify the conditions of Corollary 2.6 and, thereby,
the two-node EQMOM is well-defined. Moreover, a direct calculation via (2.10)
gives $b_{1}=-2$, $b_{2}=1$ and $b_{j}=0$ for $j\geq 3$. Thus, the polynomial
in Theorem 2.8 for the two-node case is $p(t)=t^{3}(t-1)^{2}$ and has two
nonzero roots. By Theorem 2.8, the two-node EQMOM moment system is strictly
hyperbolic.

###  6.2. Numerical validation

In this subsetion, we use a Riemann problem of the Euler equations to show
that the kernel functions given in the previous subsection produce
satisfactory results if they satisfy the conditions required by our theory;
otherwise they lead to spurious results. The initial data of the Riemann
problem are

$\rho(0,x)=\left\\{\begin{aligned} 3.093,\quad&x<0,\\\
1,\quad&x>0,\end{aligned}\right.\quad U(0,x)=0,\quad\theta(0,x)=1,$

for the Euler equation. These data are used to determine the equilibrium
distribution $f^{eq}(\rho,U,\theta;\xi)$ in the kinetic equation (2.1) and
thereby initial moments.

The computational domain $-1\leq x\leq 1$ is discretized into 1000 uniform
cells. The Neumann boundary condition $\partial_{x}f=0$ is applied on the
endpoints $x=\pm 1$. The spatial fluxes are treated as in [6] and the time
step is chosen so that the CFL number is less than 0.5. Two limiting cases of
$\tau$ in (2.7) are considered here. The continuum limit has infinitely fast
collisions with $\tau=0$, while the free-molecular limit assumes no collision
($\tau=\infty$). The analytical solutions for the two cases can be found in
[24] and [6], respectively. We test six different kernels for the both cases.

Figure 1. 1-D Riemann problem with no collision: Profiles of density $\rho$,
velocity $U$ and temperature $\theta$ at $t=0.1$. Kernels: 1\. Analytical
solution; 2. Gaussian; 3. Kappa distribution with $\kappa=6$; 4. $K_{3}(\xi)$
in Example 6.2; 5. $\mathcal{K}(\xi)$ in Example 6.6; 6. $\mathcal{K}(\xi)$ in
Example 6.5; 7. $K_{25}(\xi)$ in Example 6.2. Figure 2. 1-D Riemann problem in
the continuum limit: Profiles of density $\rho$, velocity $U$ and temperature
$\theta$ at $t=0.1$. Kernels: 1\. Analytical solution; 2. Gaussian; 3. Kappa
distribution with $\kappa=6$; 4. $K_{3}(\xi)$ in Example 6.2; 5.
$\mathcal{K}(\xi)$ in Example 6.6; 6. $\mathcal{K}(\xi)$ in Example 6.5; 7.
$K_{25}(\xi)$ in Example 6.2.

Fig. 1 shows the spatial profiles of the macroscopic quantities
$(\rho,u,\theta)$ at $t=0.1$ for the free-molecular case. Both the simulated
results and analytical solutions are plotted. It is seen that the kernel
greatly affects the simulation results of this highly non-equilibrium flow.
The Gaussian kernel, the kappa distribution in Example 6.3 with $\kappa=6$ and
the even polynomial $\mathcal{K}_{3}(\xi)$ in Example 6.2 all satisfy the
structural stability condition, and exhibit similar precision in the
simulation. The results with the uneven kernel in Example 6.6 can still
roughly reproduce the profiles but shows larger errors than the former
kernels. Notice that the corresponding moment system is strictly hyperbolic,
while whether the structural stability condition (III) is respected remains
unclear.

On the other hand, Figs. 1 (b1)-(b3) include the results from two improper
kernels. The kernel in Example 6.5 yields a non-hyperbolic system, leading to
huge unphysical peaks (termed ‘$\delta$-shocks’) in the region where flow
quantities change drastically. This is a common phenomenon for non-hyperbolic
moment systems [10]. By contrast, the kernel $\mathcal{K}_{25}(\xi)$ in
Example 6.2 is hyperbolic but violates the stability condition (III). The
errors are larger than those from the kernel $\mathcal{K}_{3}(\xi)$.

Fig. 2 gives the numerical results of the continuum case. The flow can be
reasonably predicted, including a right-moving shock wave, a left-moving
rarefaction wave and a discontinuity at $x=0$. The result indicates that the
continuum flow regime may be less sensitive to the choice of the kernels.

##  7\. Conclusions

This paper is concerned with a class of moment closure systems derived with an
extended quadrature method of moments (EQMOM) for the one-dimensional BGK
equation. The class is characterized with a kernel function and the unknown
distribution is approximated with Ansatz Eq.(2.4). We investigate the
realizability of the extended method of moments (see Theorem 2.3 & Corollary
2.6). A sufficient condition (Theorem 2.8) on the kernel is identified for the
EQMOM-derived moment systems to be strictly hyperbolic.

Furthermore, sufficient and necessary conditions are established for the two-
node systems to be well-defined and strictly hyperbolic, and to preserve the
dissipation property of the kinetic equation. For normalized kernels, the
condition is $\mathfrak{m}_{4}\geq 3+\frac{9}{8}\mathfrak{m}_{3}^{2}$ for the
well-definedness, where $\mathfrak{m}_{3}$ and $\mathfrak{m}_{4}$ are the 3rd
and 4th moments of the kernel function. When the kernel is even, the
conditions are $3\leq\mathfrak{m}_{4}<6$ for hyperbolicity and
$3\leq\mathfrak{m}_{4}<5$ for the dissipativeness corresponding to the
$H$-theorem of the kinetic equation.

In addition, we present a number of examples of the kernel functions and
examine their performance numerically. Precisely, we use a Riemann problem of
the Euler equations to show that the kernel functions produce satisfactory
results if they satisfy the conditions required by our theoretical results;
otherwise they lead to spurious results.

## Acknowledgments

This work is supported by the National Key Research and Development Program of
China (Grant no. 2021YFA0719200) and the National Natural Science Foundation
of China (Grant no. 12071246). The authors are grateful to Prof. Shuiqing Li
and Mr. Yihong Chen at Tsinghua University for insightful discussions.

## References

  * [1] G. W. Alldredge, R. Li, W. Li, Approximating the m-2 method by the extended quadrature method of moments for radiative transfer in slab geometry, Kinet. Relat. Models, 9 (2016) 237–249.
  * [2] P. L. Bhatnagar, E. P. Gross, M. Krook, A model for collision processes in gases. I. small amplitude processes in charged and neutral one-component systems, Phys. Rev., 94 (1954) 511–525.
  * [3] Z. Cai, Y. Fan, R. Li, Globally hyperbolic regularization of Grad’s moment system, Commun. Pur. Appl. Math., 67 (2013) 464–518.
  * [4] C. Chalons, R. O. Fox, M. Massot, A multi-Gaussian quadrature method of moments for gas-particle flows in a LES framework, in Proceedings of the Summer Program 2010, Center for Turbulence Research, Stanford University, Stanford, CA, (2010) 347–358.
  * [5] C. Chalons, D. Kah, M. Massot, Beyond pressureless gas dynamics: Quadrature-based velocity moment models, Commun. Math. Sci., 10 (2012) 1241–1272.
  * [6] C. Chalons, R. Fox, F. Laurent, M. Massot, A. Vi$\acute{e}$, Multivariate Gaussian extended quadrature method of moments for turbulent disperse multiphase flow, Multiscale Model. Simul. 15 (2017) 1553–1583.
  * [7] Y. Cheng, J.A. Rossmanith, A class of quadrature-based moment-closure methods with application to the Vlasov-Poisson-Fokker-Planck system in the high-field limit, J. Comput. Appl. Math., 262 (2014) 384–398.
  * [8] Y. Di, Y. Fan, R. Li, L. Zheng, Linear stability of hyperbolic moment models for Boltzmann equation, Numer. Math. Theor. Meth. Appl., 10 (2017) 255–277.
  * [9] G. Dimarco, L. Pareschi, Numerical methods for kinetic equations, Acta Numerica, 23 (2014) 369–520.
  * [10] R. Fox, A quadrature-based third-order moment method for dilute gas-particle flows, J. Comput. Phys. 227 (2008) 6313–6350.
  * [11] H. Grad, On the kinetic theory of rarefied gases, Comm. Pure Appl. Math., 2 (1949), 331–407.
  * [12] S. Harris, An Introduction to the Theory of the Boltzmann Equation, Dover Publications, New York, 2004.
  * [13] Q. Huang, S.Q. Li, W.-A. Yong, Stability analysis of quadrature-based moment methods for kinetic equations, SIAM J. Appl. Math. 80 (2020) 206–231.
  * [14] Q. Huang, J. Koellermeier, W.-A. Yong, Equilibrium stability analysis of hyperbolic shallow water moment equations, Math. Meth. Appl. Sci. 45 (2022) 6459–6480.
  * [15] J. Koellermeier and M. Rominger, Analysis and numerical simulation of hyperbolic shallow water moment equations, Commun. Comp. Phys., 28 (2020), 1038–1084.
  * [16] D. L. Marchisio, R. O. Fox, Computational Models for Polydisperse Particulate and Multiphase Systems, Cambridge University Press, Cambridge, 2013.
  * [17] R. McGraw, Description of aerosol dynamics by the quadrature method of moments, Aerosol Sci. Technol., 27 (1997), 255–265.
  * [18] L. Mieussens, Discrete velocity model and implicit scheme for the BGK equation of rarefield gas dynamics, Math. Mod. Meth. Appl. Sci., 10 (2000) 1121–1149.
  * [19] V. Pierrard and M. Lazar, Kappa Distributions: Theory and Applications in Space Plasma, Sol. Phys., 267 (2010) 153–174.
  * [20] M. Pigou, J. Morchain, P. Fede, M. I. Penet, G. Laronze, New developments of the Extended Quadrature Method of Moments to solve Population Balance Equations, J. Comput. Phys., 365 (2018) 243–268.
  * [21] K. Schmüdgen, The Moment Problem. Graduate Texts in Mathematics, vol 277, Springer, Cham, 2017.
  * [22] D. Serre, Systems of Conservation Laws 1: Hyperbolicity, Entropies, Shock Waves, Cambridge University Press, Cambridge, 1999.
  * [23] H. Struchtrup and M. Torrilhon, Regularization of Grad’s 13 moment equaitons: Derivation and linear analysis, Phys. Fluids, 15 (2003) 2668–2680.
  * [24] E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd edition, Springer, New York, 2009.
  * [25] W.-A. Yong, Singular perturbations of fisrt-order hyperbolic systems with stiff source terms, J. Differ. Equations, 155 (1999) 89–132.
  * [26] W.-A. Yong, Basic Aspects of Hyperbolic Relaxation Systems, in Advances in the Theory of Shock Waves, Birkhäuser Boston, Boston, MA, (2001) 259–305.
  * [27] C. Yuan, F. Laurent, R.O. Fox, An extended quadrature method of moments for population balance equations, J. Aerosol Sci., 51 (2012) 1–23.
  * [28] W. Zhao, W.-A. Yong, L.-S. Luo, Stability analysis of a class of globally hyperbolic moment systems, Commun. Math. Sci., 15 (2017) 609–633.