# Discrete-velocity-direction models of BGK-type with minimum entropy: I.
Basic idea

Qian Huang  Department of Energy and Power Engineering, Tsinghua University
Beijing, 100084, China<EMAIL_ADDRESS>,  Yihong Chen  Department
of Mathematical Sciences, Tsinghua University
Beijing, 100084, China<EMAIL_ADDRESS>and  Wen-An Yong
Department of Mathematical Sciences, Tsinghua University
Beijing, 100084, China<EMAIL_ADDRESS>
###### Abstract.

In this series of works, we develop a discrete-velocity-direction model (DVDM)
with collisions of BGK-type for simulating rarefied flows. Unlike the
conventional kinetic models (both BGK and discrete-velocity models), the new
model restricts the transport to finite fixed directions but leaves the
transport speed to be a 1-D continuous variable. Analogous to the BGK
equation, the discrete equilibriums of the model are determined by minimizing
a discrete entropy. In this first paper, we introduce the DVDM and investigate
its basic properties, including the existence of the discrete equilibriums and
the $H$-theorem. We also show that the discrete equilibriums can be
efficiently obtained by solving a convex optimization problem. The proposed
model provides a new way in choosing discrete velocities for the computational
practice of the conventional discrete-velocity methodology. It also
facilitates a convenient multidimensional extension of the extended quadrature
method of moments. We validate the model with numerical experiments for two
benchmark problems at moderate computational costs.

###### Key words and phrases:

Kinetic equations; discrete-velocity-direction model (DVDM); minimum entropy
principle; discrete-velocity model (DVM); extended quadrature method of
moments (EQMOM)

* Corresponding author

##  1\. Introduction

Kinetic theories and the Boltzmann equation lay the foundation for
investigating non-equilibrium many-body interacting systems. Besides the
classical rarefied gas dynamics [30], kinetic theories have also found
substantial applications in multiphase flow problems [13, 19] and the emerging
field of active matter [3, 31]. However, solving the Boltzmann-like equation
can be challenging. The first obstacle is due to the collision mechanisms of
gas molecules [16] or aggregation-breakage of aerosol/colloid particles [13,
23]. To overcome this obstacle, the BGK model [5] and its improvements (for
instance, Shakhov model [29] and ES model [17]) have been proposed to simplify
the collision term while keeping fundamental properties of the original
equation and hence become prevalent [16].

The second obstacle is the velocity dependence of the unknowns which makes
solving the Boltzmann equation computationally costly. Many numerical
approaches have thus been developed, including the linearization method [27,
30], the spectral method [11], the discrete-velocity method (DVM) [1, 14, 24,
25], the semi-continuous method [20, 28], and a variety of moment methods [6,
21, 7, 23, 19]. All these methods aim at good mathematical properties
(realizability, model stability, convergence to the original equation, etc.)
and numerical performance (computational complexity, numerical stability,
etc.), and have their advantages and weaknesses [4].

Among the methods above, the Gaussian-extended quadrature method of moments
(Gaussian-EQMOM) [7, 23, 34] serves as the original motivation of this
project. In 1-D case, this method assumes the unknown or distribution to be a
linear combination of $N$ ($N\geq 1$) Gaussian functions with unknown centers
and a common (unknown) variance to be inversely solved from the velocity
moments of the distribution [34]. The resulting moment closure system was
proved to satisfy the structural stability condition widely respected by
physical systems [18, 33]. The quadrature-based method of moments stands out
also because of its potential to characterize systems far from equilibrium
[22, 23]. However, a satisfactory extension of EQMOM to 2/3-D velocity
space(s) seems not available; see Refs. [7, 23].

On the other hand, we notice that a discrete-velocity-direction model (DVDM)
was proposed physically and investigated in Refs. [35, 36, 37]. As a semi-
continuous model different from the existing ones [20, 28], the DVDM
constrains particle transport in finite fixed directions while keeping the
speed continuous. The resulting model has complicated terms for collisions and
therefore does not offer a convenient way to get multidimensional EQMOM [35].
Additionally, the semi-continuous system can be numerically solved by
discretizing the transport speed in each orientation [35], and thus provides a
new way in selecting discrete velocity nodes, which is distinct from the
common DVM practice in a uniform cubic lattice [24].

In this project, we propose a discrete-velocity-direction model (DVDM) with
collisions of BGK-type. For the new model, the equilibriums are determined by
minimizing a discrete entropy. We investigate the existence and uniqueness of
the discrete equilibriums and establish an $H$-theorem characterizing the
dissipation property of the original kinetic equation. Because the DVDM is
different from the discrete-velocity BGK models [24], the analysis is quite
involved in comparison with that in Ref. [24]. We also show that the discrete
equilibriums can be obtained efficiently by solving a convex optimization
problem. Moreover, the continuous transport speed can be treated with simple
discretizations or combinations with 1-D Gaussian-EQMOM. The latter provides a
new multidimensional extension of EQMOM and yields a hyperbolic moment system,
which seems satisfactory. The proposed methods are numerically validated with
two benchmark tests. Further numerical and analytical results will be reported
in forthcoming papers [8].

The remainder of the paper is organized as follows. Section 2 presents the
discrete-velocity-direction model (DVDM) with BGK-collisions. The existence of
local equilibriums is studied in Section 3. Section 4 provides two approaches
to treat the continuous velocity-modulus in DVDM, including DVD-DVM in Section
4.1 and DVD-EQMOM in Section 4.2. Section 5 discusses numerical issues, with
the space-time discretization schemes in Section 5.1 and the computation of
discrete equilibrium in Section 5.2. In Section 6, two benchmark flows are
simulated. Finally, we conclude our paper in Section 7.

##  2\. Discrete-velocity-direction models of BGK-type

Let $\bm{\xi}\in\mathbb{R}^{D}$ be the molecule velocity ($D=2,3$) and
$\bm{x}\in\mathbb{R}^{D_{x}}$ the spatial position ($1\leq D_{x}\leq D$). We
consider the BGK equation [5] for the distribution $f=f(t,\bm{x},\bm{\xi})$:

$\left\\{\begin{aligned}
\partial_{t}f+\bm{\xi}\cdot\nabla_{\bm{x}}f&=\frac{1}{\tau}\left(\mathcal{E}[f]-f\right),\quad
t\geq 0,\\\
\mathcal{E}[f]&=\exp\left(\bm{\alpha}_{eq}\cdot\bm{m}(\bm{\xi})\right),\\\
\bm{\alpha}_{eq}&=\left(\ln\frac{\rho}{(2\pi\theta)^{D/2}}-\frac{U^{2}}{2\theta},\
\frac{\bm{U}}{\theta},\ -\frac{1}{\theta}\right)^{T}\in\mathbb{R}^{D+2},\\\
\rho\left(1,\bm{U},E\right)^{T}&=\langle\bm{m}(\bm{\xi})f\rangle\in\mathbb{R}^{D+2},\quad\theta=\frac{2E-U^{2}}{D}.\end{aligned}\right.$
(2.1)

Here $\tau$ is the relaxation time, $\bm{m}(\bm{\xi})=\left(1,\ \bm{\xi},\
\frac{|\bm{\xi}|^{2}}{2}\right)^{T}\in\mathbb{R}^{D+2}$, $U=|\bm{U}|$ is the
Euclidean length of the vector $\bm{U}$, and the bracket $\langle\cdot\rangle$
is defined as $\langle
g(\bm{\xi})\rangle=\int_{\mathbb{R}^{D}}g(\bm{\xi})d\bm{\xi}$ for any
measurable function $g(\bm{\xi})$. The local Maxwellian equilibrium
$\mathcal{E}[f]$ is implicitly defined by $f$ through the macroscopic fluid
density $\rho$, velocity $\bm{U}\in\mathbb{R}^{D}$, energy $E$ (or temperature
$\theta$) which are the velocity moments of $f$.

It can be easily verified that $\mathcal{E}[f]$ reproduces the local
macroscopic quantities [30]:

$\langle\bm{m}(\bm{\xi})\mathcal{E}[f]\rangle=\bm{\rho}:=\rho\left(1,\bm{U},E\right)^{T}\in\mathbb{R}^{D+2}.$
(2.2)

Moreover, it was shown [24] that given any $\bm{\rho}$ with $\rho>0$ and
$\theta>0$, $\mathcal{E}[f]$ is the unique solution that minimizes the
following kinetic entropy $H[f]$ :

$H[f]=\langle f\ln f-f\rangle$ (2.3)

subject to the constraint Eq. (2.2).

Solving the multidimensional BGK Eq. (2.1) can be computationally costly. In
this work, we propose a class of discrete-velocity-direction models (DVDM)
with a minimum entropy.

In the DVDM, the particle transport is limited to $N$ prescribed directions
$\bm{l}_{m}\in\mathbb{R}^{D}$ with $|\bm{l}_{m}|=1$ ($m=1,\dots,N$). Denote by
$\bm{l}_{\mathcal{L}}=\left(\bm{l}_{m}\right)_{m=1}^{N}$. Usually, we take
$N>D$ and choose the directions such that the matrix
$L^{T}=(\bm{l}_{1},\dots,\bm{l}_{N})\in\mathbb{R}^{D\times N}$ has rank $D$.
The velocity distribution $f(t,\bm{x},\bm{\xi})$ is replaced by $N$ 1-D
velocity distributions $f_{\mathcal{L}}=\\{f_{m}(t,\bm{x},\xi)\\}_{m=1}^{N}$
with $\xi\in\mathbb{R}$. The governing equation for each
$f_{m}=f_{m}(t,\bm{x},\xi)$ has the following form:

$\partial_{t}f_{m}+\xi\bm{l}_{m}\cdot\nabla_{\bm{x}}f_{m}=\frac{1}{\tau}(\mathcal{E}_{m}-f_{m})$
(2.4)

with the local equilibrium
$\mathcal{E_{L}}=\\{\mathcal{E}_{m}(t,\bm{x},\xi)\\}_{m=1}^{N}$ yet to be
determined by $f_{\mathcal{L}}$.

For further references, we introduce a bracket
$\langle(\cdot)_{\mathcal{L}}\rangle_{\mathcal{L}}$ for any $N$-tuple
$g_{\mathcal{L}}=\left(g_{m}(\xi)\right)_{m=1}^{N}$ as

$\langle
g_{\mathcal{L}}\rangle_{\mathcal{L}}=\sum_{m=1}^{N}\int_{\mathbb{R}}g_{m}(\xi)d\xi.$
(2.5)

Define

$\rho=\langle
f_{\mathcal{L}}\rangle_{\mathcal{L}},\quad\rho\bm{U}=\langle\xi\bm{l}_{\mathcal{L}}f_{\mathcal{L}}\rangle_{\mathcal{L}},\quad\rho
E=\langle\frac{\xi^{2}}{2}f_{\mathcal{L}}\rangle_{\mathcal{L}},$

or equivalently

$\bm{\rho}=\langle\bm{m}_{\mathcal{L}}f_{\mathcal{L}}\rangle_{\mathcal{L}},$
(2.6)

with $\bm{m}_{\mathcal{L}}=\\{\bm{m}_{m}\\}_{m=1}^{N}$ and
$\bm{m}_{m}=\left(1,\xi\bm{l}_{m},\frac{1}{2}\xi^{2}\right)^{T}\in\mathbb{R}^{D+2}$
for $m=1,\dots,N$.

The local equilibrium $\mathcal{E_{L}}$ is determined so that it minimizes the
discrete analogue of the kinetic entropy

$H[f_{\mathcal{L}}]=\langle f_{\mathcal{L}}\ln
f_{\mathcal{L}}-f_{\mathcal{L}}\rangle_{\mathcal{L}}.$ (2.7)

among all possible $N$-tuples $f_{\mathcal{L}}\geq 0$ satisfying
$\langle\bm{m}_{\mathcal{L}}f_{\mathcal{L}}\rangle_{\mathcal{L}}=\bm{\rho}$
for a given $\bm{\rho}$. The same idea was utilized to develop a conservative
and entropy-decreasing discrete-velocity model (DVM) of the BGK equation [24].

As to the existence and uniqueness of the minimizer, we have the following
result.

######  Theorem 2.1.

Suppose $f_{\mathcal{L}}=\\{f_{m}(\xi)\\}_{m=1}^{N}$ and
$\bm{\rho}=\langle\bm{m}_{\mathcal{L}}f_{\mathcal{L}}\rangle_{\mathcal{L}}$
satisfy $f_{m}\geq 0$ for all $m$ and $0<|\bm{\rho}|<\infty$. Then the
discrete kinetic entropy Eq. (2.7) has a unique minimizer $\mathcal{E_{L}}$
satisfying the constraint
$\langle\bm{m}_{\mathcal{L}}\mathcal{E_{L}}\rangle_{\mathcal{L}}=\bm{\rho}$.
Moreover, the minimizer has the exponential form

$\mathcal{E}_{m}=\exp(\bm{\alpha}\cdot\bm{m}_{m})$ (2.8)

with a certain $\bm{\alpha}\in\mathbb{R}^{D+2}$.

The proof of this theorem will be completed in the next section, where we show
that there exists
$\bm{\alpha}=(\alpha_{0},\hat{\bm{\alpha}}^{T},\alpha_{D+1})^{T}\in\mathbb{R}^{D+2}$
with $\hat{\bm{\alpha}}\in\mathbb{R}^{D}$ such that

$\langle\bm{m}_{\mathcal{L}}\exp\left(\bm{\alpha}\cdot\bm{m}_{\mathcal{L}}\right)\rangle_{\mathcal{L}}=\bm{\rho}.$
(2.9)

This implies $\alpha_{D+1}<0$. Having such an $\bm{\alpha}$, we can show that
the local equilibrium $\mathcal{E_{L}}$ with
$\mathcal{E}_{m}=\exp(\bm{\alpha}\cdot\bm{m}_{m})$ is the unique minimizer.

To do this, we use the fact that the function $(x\ln x-x)$ is strictly convex.
Then, for any $f_{\mathcal{L}}=\\{f_{m}\\}_{m=1}^{N}$ satisfying
$\langle\bm{m}_{\mathcal{L}}f_{\mathcal{L}}\rangle_{\mathcal{L}}=\bm{\rho}$,
we have

$f_{m}\ln
f_{m}-f_{m}\geq\mathcal{E}_{m}\ln\mathcal{E}_{m}-\mathcal{E}_{m}+(\ln\mathcal{E}_{m})(f_{m}-\mathcal{E}_{m})=\mathcal{E}_{m}\ln\mathcal{E}_{m}-\mathcal{E}_{m}+(\bm{\alpha}\cdot\bm{m}_{m})(f_{m}-\mathcal{E}_{m})$

for all $m$, and thereby,

$\displaystyle H[f_{\mathcal{L}}]$ $\displaystyle\geq
H[\mathcal{E}_{\mathcal{L}}]+\langle(\bm{\alpha}\cdot\bm{m}_{\mathcal{L}})(f_{\mathcal{L}}-\mathcal{E_{L}})\rangle_{\mathcal{L}}$
$\displaystyle=H[\mathcal{E}_{\mathcal{L}}]+\alpha_{0}\langle
f_{\mathcal{L}}-\mathcal{E_{L}}\rangle_{\mathcal{L}}+\hat{\bm{\alpha}}\cdot\langle\xi\bm{l}_{\mathcal{L}}(f_{\mathcal{L}}-\mathcal{E_{L}})\rangle_{\mathcal{L}}+\alpha_{D+1}\langle\frac{\xi^{2}}{2}(f_{\mathcal{L}}-\mathcal{E_{L}})\rangle_{\mathcal{L}}$
$\displaystyle=H[\mathcal{E}_{\mathcal{L}}].$

The equality holds if and only if $f_{m}=\mathcal{E}_{m}$ for all $m=1,...,N$.
Thus, if there is another minimizer of the exponential form, say,
$\exp(\bm{\alpha}^{\prime}\cdot\bm{m}_{m})$, then
$0=(\bm{\alpha}^{\prime}-\bm{\alpha})\cdot\bm{m}_{m}=(\alpha_{0}^{\prime}-\alpha_{0})+\xi\bm{l}_{m}\cdot(\hat{\bm{\alpha}}^{\prime}-\hat{\bm{\alpha}})+\frac{\xi^{2}}{2}(\alpha_{D+1}^{\prime}-\alpha_{D+1})$,
which means $\bm{\alpha}^{\prime}=\bm{\alpha}$ by choosing different $\xi$ and
using the assumption that $L^{T}=(\bm{l}_{1},\dots,\bm{l}_{N})$ is of full-
rank. Hence, the uniqueness of the minimizer is proved.

With the local equilibrium $\mathcal{E_{L}}$ determined as above, the DVDM Eq.
(2.4) is well defined as

$\left\\{\begin{aligned}
\partial_{t}f_{m}+\xi\bm{l}_{m}\cdot\nabla_{\bm{x}}f_{m}&=\frac{1}{\tau}(\exp(\bm{\alpha}\cdot\bm{m}_{m})-f_{m}),\quad
m=1,\dots,N,\\\
\langle\bm{m}_{\mathcal{L}}\exp(\bm{\alpha}\cdot\bm{m}_{\mathcal{L}})\rangle_{\mathcal{L}}&=\bm{\rho}:=\langle\bm{m}_{\mathcal{L}}f_{\mathcal{L}}\rangle_{\mathcal{L}}.\end{aligned}\right.$
(2.10)

The $N$ differential equations are coupled through $\bm{\alpha}$ determined by
the last $(D+2)$ nonlinear algebraic equations in terms of
$f_{\mathcal{L}}=\\{f_{m}\\}_{m=1}^{N}$. Note that the local equilibrium in
Eq. (2.10) can be rewritten as

$\mathcal{E}_{m}=\exp(\bm{\alpha}\cdot\bm{m}_{m})=\frac{\rho_{m}}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(\xi-
u_{m})^{2}}{2\sigma^{2}}\right)$ (2.11)

for $\alpha_{D+1}<0$. The parameter $(\rho_{m},u_{m},\sigma^{2})$ is related
to $\bm{\alpha}$ as follows:

$\sigma^{2}=-\frac{1}{\alpha_{D+1}},\quad
u_{m}=\bm{l}_{m}\cdot\left(\sigma^{2}\hat{\bm{\alpha}}\right),\quad\rho_{m}=\sqrt{2\pi}\sigma\exp\left(\alpha_{0}+\frac{u_{m}^{2}}{2\sigma^{2}}\right).$
(2.12)

######  Remark 2.2.

While $\mathcal{E_{L}}$ preserves the local fluid quantity
$\bm{\rho}=\langle\bm{m}_{\mathcal{L}}f_{\mathcal{L}}\rangle_{\mathcal{L}}$,
the conservation property generally does not hold in any specific direction.
That is to say, $\rho_{m}\neq\int_{\mathbb{R}}f_{m}d\xi$,
$\rho_{m}u_{m}\neq\int_{\mathbb{R}}\xi f_{m}d\xi$, and
$\rho_{m}(u_{m}^{2}+\sigma^{2})\neq\int_{\mathbb{R}}\xi^{2}f_{m}d\xi$ for
$m=1,\dots,N$. Indeed, this is the key mechanism reallocating molecules among
the prescribed directions in the DVDM-BGK model Eq. (2.10).

The DVDM-BGK model Eq. (2.10) has also the following properties including an
$H$-theorem.

######  Theorem 2.3.

Suppose the DVDM-BGK Eq. (2.10) with positive initial data
$f_{m}(0,\bm{x},\xi)>0$ has a solution
$f_{\mathcal{L}}=\\{f_{m}(t,\bm{x},\xi)\\}_{m=1}^{N}$. Then we have

$\displaystyle f_{m}(t,\bm{x},\xi)>0,\quad\forall\ m,t,$
$\displaystyle\bm{x},\xi,$ (2.13)
$\displaystyle\partial_{t}\langle\bm{m}_{\mathcal{L}}f_{\mathcal{L}}\rangle_{\mathcal{L}}+\langle\bm{m}_{\mathcal{L}}\xi\bm{l}_{\mathcal{L}}\cdot\nabla_{\bm{x}}f_{\mathcal{L}}\rangle_{\mathcal{L}}$
$\displaystyle=0,$ (2.14)
$\displaystyle\partial_{t}H[f_{\mathcal{L}}]+\langle\xi\bm{l}_{\mathcal{L}}\cdot\nabla_{\bm{x}}(f_{\mathcal{L}}\ln
f_{\mathcal{L}}-f_{\mathcal{L}})\rangle_{\mathcal{L}}$ $\displaystyle\leq 0.$
(2.15)

###### Proof.

For Eq. (2.13), we define $g_{m}(t,\bm{x},\xi)=f_{m}(t,\bm{x}+\xi
t\bm{l}_{m},\xi)$ and deduce from Eq. (2.10) that

$\displaystyle\frac{d}{dt}g_{m}(t,\bm{x},\xi)$
$\displaystyle=\partial_{t}f_{m}(t,\bm{x}+\xi
t\bm{l}_{m},\xi)+\xi\bm{l}_{m}\cdot\nabla_{\bm{x}}f_{m}(t,\bm{x}+\xi
t\bm{l}_{m},\xi)$
$\displaystyle=\frac{1}{\tau}\left(\mathcal{E}_{m}(t,\bm{x}+\xi
t\bm{l}_{m},\xi)-g_{m}(t,\bm{x},\xi)\right),$

yielding

$\frac{d}{dt}\left(e^{\frac{t}{\tau}}g_{m}(t,\bm{x},\xi)\right)=\frac{1}{\tau}e^{\frac{t}{\tau}}\mathcal{E}_{m}(t,\bm{x}+\xi
t\bm{l}_{m},\xi).$

Solving $g_{m}(t,\bm{x},\xi)$ from the above equation, we obtain

$f_{m}(t,\bm{x},\xi)=\frac{1}{\tau}\int_{0}^{t}e^{\frac{s-t}{\tau}}\mathcal{E}_{m}(s,\bm{x}+\xi(s-t)\bm{l}_{m},\xi)\mathrm{d}s+e^{\frac{-t}{\tau}}f_{m}(0,\bm{x}-\xi
t\bm{l}_{m},\xi).$

Thus, Eq. (2.13) follows immediately from the positivity of
$f_{m}(0,\bm{x},\xi)$ and $\mathcal{E}_{m}=\exp(\bm{\alpha}\cdot\bm{m}_{m})$.

To derive Eq. (2.14), we just exert the operation
$\langle\bm{m}_{\mathcal{L}}(\cdot)_{\mathcal{L}}\rangle_{\mathcal{L}}$ on
both sides of the differential equations of Eq. (2.10) and notice that the RHS
vanishes as
$\langle\bm{m}_{\mathcal{L}}\left(\exp(\bm{\alpha}\cdot\bm{m}_{\mathcal{L}})-f_{\mathcal{L}}\right)\rangle_{\mathcal{L}}=0$.

For Eq. (2.15), we multiply the two sides of Eq. (2.4) with $\ln f_{m}$ to
obtain

$\displaystyle\partial_{t}(f_{m}\ln
f_{m}-f_{m})+\xi\bm{l}_{m}\cdot\nabla_{\bm{x}}(f_{m}\ln f_{m}-f_{m})$
$\displaystyle=\frac{1}{\tau}(\mathcal{E}_{m}-f_{m})(\ln f_{m})$
$\displaystyle\leq\frac{1}{\tau}(\mathcal{E}_{m}-f_{m})(\ln\mathcal{E}_{m}),$

where the inequality $(x-y)(\ln x-\ln y)\geq 0$ has been used. Then, Eq.
(2.15) follows by applying $\langle(\cdot)\rangle_{\mathcal{L}}$ to both sides
of the last inequality and noticing that
$\langle(\mathcal{E}_{\mathcal{L}}-f_{\mathcal{L}})(\ln\mathcal{E}_{\mathcal{L}})\rangle_{\mathcal{L}}=0$.
This completes the proof. ∎

Note that Eq. (2.14) is just the classical conservation laws for the
macroscopic fluid quantities:

$\partial_{t}\bm{\rho}+\nabla_{\bm{x}}\cdot\begin{pmatrix}\rho\bm{U}\\\
\rho\bm{U}\otimes\bm{U}+\mathbf{P}\\\
\rho\bm{U}E+\mathbf{P}\cdot\bm{U}+\bm{q}\end{pmatrix}=\bm{0}.$

Here $\mathbf{P}=\left\langle(\xi\bm{l}_{\mathcal{L}}-\bm{U})^{\otimes
2}f_{\mathcal{L}}\right\rangle_{\mathcal{L}}$ and
$\bm{q}=\frac{1}{2}\langle(\xi\bm{l}_{\mathcal{L}}-\bm{U})|\xi\bm{l}_{\mathcal{L}}-\bm{U}|^{2}f_{\mathcal{L}}\rangle_{\mathcal{L}}$.

######  Remark 2.4 (Planar flows).

For real-world planar rarefied flows, we have $D=3$,
$\bm{\xi}=(\xi_{x},\xi_{y},\xi_{z})^{T}$ and $D_{x}=2$. In DVDM, it is
straightforward to select velocity orientations $\bm{l}_{m}$ in
$\mathbb{R}^{3}$ to obtain the governing Eq. (2.4) [35]. On the other hand,
before discretizing velocity directions, we can also resort to the technique
of reduced distribution functions to eliminate the dependence of $\xi_{z}$.
There are two approaches for this purpose: (i) Define
$(g,h)^{T}=\int(1,\xi_{z}^{2})^{T}fd\xi_{z}$ [9]; (ii) Define
$(g,h)^{T}=\frac{1}{\sqrt{\pi}}\int\left(e^{-\xi_{z}^{2}},\left(\xi_{z}^{2}-\frac{1}{2}\right)e^{-\xi_{z}^{2}}\right)^{T}fd\xi_{z}$
[30]. In both cases we have $g=g(t,\bm{x},\xi_{x},\xi_{y})$ and
$h=h(t,\bm{x},\xi_{x},\xi_{y})$. The governing equations for $g$ and $h$ can
be derived from the BGK equation and the orientations $\bm{l}_{m}$ are hence
chosen on the $(\xi_{x},\xi_{y})$-plane. The discrete equilibriums can be
modeled and solved out of the minimum entropy principle in a similar manner.

We conclude this section with the following proposition for an alternative way
to solve $\bm{\alpha}$ from the nonlinear algebraic equations Eq. (2.9).

######  Proposition 2.5.

$\bm{\alpha}$ solves Eq. (2.9) if and only if $\bm{\alpha}$ minimizes

$J(\bm{\alpha}):=\langle\exp\left(\bm{\alpha}\cdot\bm{m}_{\mathcal{L}}\right)\rangle_{\mathcal{L}}-\bm{\rho}\cdot\bm{\alpha}$
(2.16)

for $\bm{\alpha}=(\alpha_{0},\hat{\bm{\alpha}},\alpha_{D+1})$ with
$\alpha_{D+1}<0$.

###### Proof.

If $\bm{\alpha}$ is a minimizer, we have $\nabla J(\bm{\alpha})=0$, which is
just Eq. (2.9). Conversely, if $\bm{\alpha}$ solves Eq. (2.9), then we have
$\nabla J(\bm{\alpha})=0$. Thus, it suffices to show that $J=J(\bm{\alpha})$
is strictly convex. To do so, we compute the Hessian

$\nabla_{\bm{\alpha}\bm{\alpha}}J=\langle\bm{m}_{\mathcal{L}}\otimes\bm{m}_{\mathcal{L}}\exp\left(\bm{\alpha}\cdot\bm{m}_{\mathcal{L}}\right)\rangle_{\mathcal{L}}$

and consider the quadratic form

$\bm{\eta}^{T}\nabla_{\bm{\alpha}\bm{\alpha}}J\bm{\eta}=\langle(\bm{m}_{\mathcal{L}}\cdot\bm{\eta})^{2}\exp\left(\bm{\alpha}\cdot\bm{m}_{\mathcal{L}}\right)\rangle_{\mathcal{L}}\geq
0$

for any $\bm{\eta}\in\mathbb{R}^{D+2}$. The equality holds if and only if
$\bm{m}_{m}\cdot\bm{\eta}=0$ for all $m=1,\dots,N$ and thus $\bm{\eta}=\bm{0}$
due to the assumption that $L^{T}=(\bm{l}_{1},...,\bm{l}_{N})$ is of full-
rank. This indicates that the Hessian $\nabla_{\bm{\alpha}\bm{\alpha}}J$ is
positive definite and therefore $J$ is strictly convex. This completes the
proof. ∎

######  Remark 2.6.

This proposition does not claim the existence of the solution to Eq. (2.9).
But the strict convexity of $J=J(\bm{\alpha})$ implies the uniqueness.

##  3\. Existence of $\alpha$

This section is devoted to completing the proof of Theorem 2.1. Namely, we
will prove that $J(\bm{\alpha})$ defined in Eq. (2.16) attains its minimum in
the open set
$\mathcal{D}:=\mathbb{R}\times\mathbb{R}^{D}\times\mathbb{R}_{-}$. Throughout
this section, we assume that the given $f_{\mathcal{L}}$ and
$\bm{\rho}=\langle\bm{m}_{\mathcal{L}}f_{\mathcal{L}}\rangle_{\mathcal{L}}$
satisfy the constraints in Theorem 2.1: $f_{m}\geq 0$ for all $m$ and
$0<|\bm{\rho}|<\infty$.

First of all, we give a further constraint on $\bm{\rho}$.

######  Lemma 3.1.

Let $\bm{\rho}=(\rho,\rho\bm{U},\rho E)^{T}\in\mathbb{R}^{D+2}$ be computed
from $f_{\mathcal{L}}=\\{f_{m}\\}_{m=1}^{N}$ in terms of
$\bm{\rho}=\langle\bm{m}_{\mathcal{L}}f_{\mathcal{L}}\rangle_{\mathcal{L}}$.
Then there is a constant $s_{min}\geq 1$ such that

$\rho>0,\quad E>\frac{1}{2}U^{2}s_{min}^{2}.$

###### Proof.

Set $\hat{\rho}_{m}=\int_{\mathbb{R}}f_{m}d\xi\geq 0$ and
$\hat{\rho}_{m}\hat{u}_{m}=\int_{\mathbb{R}}\xi f_{m}d\xi$ for $m=1,\dots,N$.
We deduce from
$\bm{\rho}=\langle\bm{m}_{\mathcal{L}}f_{\mathcal{L}}\rangle_{\mathcal{L}}$
and the Cauchy-Schwartz inequality that

$\rho=\sum_{m=1}^{N}\hat{\rho}_{m}>0,\quad\rho\bm{U}=\sum_{m=1}^{N}\hat{\rho}_{m}\hat{u}_{m}\bm{l}_{m},$
(3.1) $2\rho E=\sum_{m=1}^{N}\int_{\mathbb{R}}\xi^{2}f_{m}d\xi>\sum_{m:\
\hat{\rho}_{m}>0}\frac{\left(\int_{\mathbb{R}}\xi
f_{m}d\xi\right)^{2}}{\int_{\mathbb{R}}f_{m}d\xi}=\sum_{m=1}^{N}\hat{\rho}_{m}\hat{u}_{m}^{2}.$

The inequality is strict because the equality holds only when $f_{m}=0$ for
all $m$, which contradicts $\rho>0$. Since

$\displaystyle\sum_{m=1}^{N}\hat{\rho}_{m}\hat{u}_{m}^{2}$
$\displaystyle=\frac{1}{\rho}\sum_{m,n}\hat{\rho}_{n}\hat{\rho}_{m}\hat{u}_{m}^{2}=\frac{1}{\rho}\left[\sum_{m}(\hat{\rho}_{m}\hat{u}_{m})^{2}+\sum_{m<n}\hat{\rho}_{m}\hat{\rho}_{n}(\hat{u}_{m}^{2}+\hat{u}_{n}^{2})\right]$
$\displaystyle\geq\frac{1}{\rho}\left(\sum_{m=1}^{N}|\hat{\rho}_{m}\hat{u}_{m}|\right)^{2},$

we have

$2\rho
E>\frac{1}{\rho}\left(\sum_{m=1}^{N}|\hat{\rho}_{m}\hat{u}_{m}|\right)^{2}.$
(3.2)

On the other hand, the second equation in Eq. (3.1) can be rewritten as

$\sum_{m=1}^{N}\left(\frac{\hat{\rho}_{m}\hat{u}_{m}}{\rho
U}\right)\bm{l}_{m}=\bm{e}:=\frac{\bm{U}}{U}.$

Here we have assumed that $U\neq 0$, otherwise the proof is already complete.
Denote by $s_{min}=s_{min}(\bm{e})$ the minimum of $\sum_{m=1}^{N}|a_{m}|$
subject to the constraints
$\sum_{m=1}^{N}a_{m}\bm{l}_{m}=\bm{e}\in\mathbb{R}^{D}$. Obviously we have
$s_{min}(\bm{e})\geq 1$ because

$\sum_{m=1}^{N}|a_{m}|=\sum_{m=1}^{N}|a_{m}\bm{l}_{m}|\geq\left|\sum_{m=1}^{N}a_{m}\bm{l}_{m}\right|=|\bm{e}|=1.$

Hence we see that $\sum_{m=1}^{N}|\hat{\rho}_{m}\hat{u}_{m}|\geq\rho
Us_{min}(\bm{e})$ and the proof is complete. ∎

######  Remark 3.2.

Lemma 3.1 indicates the capacity of DVDM to realize macroscopic flow states.
In the Boltzmann equation, the constraint for $\bm{\rho}$ is $\rho>0$ and
$E>\frac{1}{2}U^{2}$ (because the temperature $\theta=(2E-U^{2})/D>0$). By
contrast, Lemma 3.1 gives a bigger lower bound for $E$ because $s_{min}\geq
1$, exhibiting the price we pay in the DVDM approximation. Generally speaking,
more states can be realized with more discrete velocity directions (see
Appendix A for the case $D=2$). On the other hand, no upper bound is required
for $E$ or $\theta$ in DVDM, which is in contrast to DVM [24].

To show the existence of a minimizer $\bm{\alpha}$ of $J(\bm{\alpha})$, we
take $\bm{\alpha}^{*}\in\mathcal{D}$ and set
$B_{M}:=\\{\bm{\alpha}\in\mathcal{D}:J(\bm{\alpha})\leq
J(\bm{\alpha}^{*})\\}$. Thus it suffices to prove that $B_{M}$ is compact,
namely, $B_{M}$ is both closed in $\mathbb{R}^{D+2}$ and bounded.

The closedness is due to the following lemma, which is similar to Proposition
2 (P2) in Ref. [26] but its proof needs more efforts.

######  Lemma 3.3.

If a sequence $\\{\bm{\alpha}^{(p)}\\}_{p=1}^{\infty}\subset\mathcal{D}$
satisfies

$\lim_{p\to\infty}\alpha_{D+1}^{(p)}=0,$

then there is a subsequence $\\{\bm{\alpha}^{(p_{k})}\\}_{k=1}^{\infty}$ such
that $J(\bm{\alpha}^{(p_{k})})$ goes to $+\infty$.

###### Proof.

Note that
$\|\bm{b}\|_{\mathcal{L}}=\sqrt{\sum_{m=1}^{N}(\bm{l}_{m}\cdot\bm{b})^{2}}$ is
a norm on $\mathbb{R}^{D}$. The reason for this is that
$L^{T}=(\bm{l}_{1},...,\bm{l}_{N})\in\mathbb{R}^{D\times N}$ is full-rank,
$L^{T}L$ is positive definite, and
$\sum_{m=1}^{N}(\bm{l}_{m}\cdot\bm{b})^{2}=\bm{b}^{T}L^{T}L\bm{b}$. Thus there
exists a constant $C>0$ such that $\|\bm{b}\|_{\mathcal{L}}^{2}\geq
C|\bm{b}|^{2}$ for all $\bm{b}\in\mathbb{R}^{D}$. In what follows, we use $C$
as a generic constant.

A direct computation from Eq. (2.16) yields

$J(\bm{\alpha})=\sum_{m=1}^{N}\sqrt{\frac{2\pi}{-\alpha_{D+1}}}\exp\left(\alpha_{0}+\frac{(\hat{\bm{\alpha}}\cdot\bm{l}_{m})^{2}}{2(-\alpha_{D+1})}\right)-\rho(\alpha_{0}+\hat{\bm{\alpha}}\cdot\bm{U}+\alpha_{D+1}E).$
(3.3)

Since $e^{x}$ is a convex function, we have

$\sum_{m=1}^{N}\exp\left(\frac{(\hat{\bm{\alpha}}\cdot\bm{l}_{m})^{2}}{2(-\alpha_{D+1})}\right)\geq
N\exp\left(\frac{\sum_{m=1}^{N}(\hat{\bm{\alpha}}\cdot\bm{l}_{m})^{2}}{2N(-\alpha_{D+1})}\right)\geq
C\exp\left(C\frac{|\hat{\bm{\alpha}}|^{2}}{|\alpha_{D+1}|}\right)$

and therefore

$J(\bm{\alpha})\geq\frac{C}{\sqrt{|\alpha_{D+1}|}}\exp\left(\alpha_{0}+C\frac{|\hat{\bm{\alpha}}|^{2}}{|\alpha_{D+1}|}\right)-\rho(\alpha_{0}+U|\hat{\bm{\alpha}}|).$
(3.4)

With this inequality, the rest of this proof is to show the right-hand side of
Eq. (3.4) goes to $+\infty$ as $\alpha_{0}^{(p)}\to 0^{-}$, which can be found
in Ref. [26] and we omit it here. ∎

It remains to show that $B_{M}$ is bounded. Otherwise, there is an unbounded
sequence $\\{\bm{\alpha}^{(p)}\\}_{p=1}^{\infty}\subset B_{M}$. Thus we only
need to show that $J(\bm{\alpha}^{(p)})$ has no upper bound. To do this, we
formulate the following lemma.

######  Lemma 3.4.

Set

$\mathbb{S}^{D+1}_{-}:=\\{\bm{w}=(w_{0},\hat{\bm{w}},w_{D+1})\in\mathbb{R}^{D+2}:|\bm{w}|=1,w_{D+1}<0\\}.$

For any $\bm{w}\in\mathbb{S}^{D+1}_{-}$, there is an open neighborhood
$\Omega(\bm{w})$ of $\bm{w}$ in $\mathbb{S}^{D+1}_{-}$ such that
$J(t\bm{w}^{\prime})\to+\infty$ uniformly on $\Omega(\bm{w})$, namely,

$\lim_{t\to+\infty}\inf_{\bm{w}^{\prime}\in\Omega(\bm{w})}J(t\bm{w}^{\prime})=+\infty.$

###### Proof.

From Eq. (3.3) we have

$J(t\bm{w})=\sum_{m=1}^{N}\sqrt{\frac{2\pi}{t|w_{D+1}|}}\exp\left(tw_{0}+t\frac{(\hat{\bm{w}}\cdot\bm{l}_{m})^{2}}{2|w_{D+1}|}\right)-t\bm{\rho}\cdot\bm{w}.$
(3.5)

If there exists $m$ such that
$(\hat{\bm{w}}\cdot\bm{l}_{m})^{2}>2w_{0}w_{D+1}$, then $\bm{w}$ has a
neighborhood $\Omega(\bm{w})$ in $\mathbb{S}^{D+1}_{-}$ such that for all
$\bm{w}^{\prime}\in\Omega(\bm{w})$,

$w_{0}^{\prime}+\frac{(\hat{\bm{w}}^{\prime}\cdot\bm{l}_{m})^{2}}{2|w_{D+1}^{\prime}|}\geq
d>0\quad\text{and}\quad\bm{\rho}\cdot\bm{w}^{\prime}\leq C(\bm{w})<+\infty.$

Then it follows from Eq. (3.5) that
$J(t\bm{w}^{\prime})\geq\frac{C}{\sqrt{t}}\exp(td)-tC(\bm{w})\to+\infty$
uniformly for $\bm{w}^{\prime}\in\Omega(\bm{w})$.

Otherwise, we have $(\hat{\bm{w}}\cdot\bm{l}_{m})^{2}\leq 2w_{0}w_{D+1}$ for
all $m$. If $w_{0}=0$, $\hat{\bm{w}}$ must be $\bm{0}$ and we have
$\bm{\rho}\cdot\bm{w}=\rho(w_{0}+\bm{U}\cdot\hat{\bm{w}}+Ew_{D+1})<0$.
Otherwise, $w_{0}$ must be negative. Due to Lemma 3.1, we can find
$\\{a_{m}\\}_{m=1}^{N}\subset\mathbb{R}$ satisfying
$U\sum_{m=1}^{N}a_{m}\bm{l}_{m}=\bm{U}$ and
$\sqrt{2E}>U\sum_{m=1}^{N}|a_{m}|$. Then we deduce that

$\displaystyle-w_{0}-Ew_{D+1}$ $\displaystyle\geq
2\sqrt{Ew_{0}w_{D+1}}>U\sum_{m=1}^{N}|a_{m}|\sqrt{2w_{0}w_{D+1}}$
$\displaystyle\geq U\sum_{m=1}^{N}|a_{m}||\hat{\bm{w}}\cdot\bm{l}_{m}|\geq
U|\sum_{m=1}^{N}a_{m}\hat{\bm{w}}\cdot\bm{l}_{m}|=|\hat{\bm{w}}\cdot\bm{U}|.$

Therefore, we also have $\bm{\rho}\cdot\bm{w}<0$. Now we can find a
neighborhood $\Omega(\bm{w})$ of $\bm{w}$ in $\mathbb{S}^{D+1}_{-}$ such that
$\bm{\rho}\cdot\bm{w}^{\prime}\leq-d<0$ for all
$\bm{w}^{\prime}\in\Omega(\bm{w})$. Thus we have
$J(t\bm{w}^{\prime})>td\to+\infty$ uniformly for
$\bm{w}^{\prime}\in\Omega(\bm{w})$. This completes the proof. ∎

######  Remark 3.5.

Since $\mathbb{S}_{-}^{D+1}$ is not compact, this lemma cannot directly imply
that $J(t\bm{w})$ goes to infinity uniformly for all directions $\bm{w}$. In
contrast, in DVM [24], the compactness of $\mathbb{S}^{D+1}$ helps to find a
finite open cover and prove the coercivity.

Having Lemma 3.4, we first suppose $\\{\bm{\alpha}^{(p)}\\}$ is contained in a
cone $C_{r}$ defined by

$C_{r}:=\left\\{\bm{\alpha}\in\mathcal{D}:\sqrt{\alpha_{0}^{2}+|\hat{\bm{\alpha}}|^{2}}\leq
r|\alpha_{D+1}|\right\\}$

for some $r>0$. Then $\mathbb{S}_{r}^{D+1}:=\mathbb{S}_{-}^{D+1}\cap C_{r}$ is
a compact set. Lemma 3.4 implies that for any $\bm{w}\in\mathbb{S}_{r}^{D+1}$,
there exist a neighborhood $\Omega(\bm{w})$ in $\mathbb{S}^{D+1}_{r}$ and
$T_{\bm{w}}>0$, such that

$J(t\bm{w}^{\prime})>J(\bm{\alpha}^{*})+1\ \text{for all }t>T_{\bm{w}},\
\bm{w}^{\prime}\in\Omega(\bm{w}).$

Since $\cup_{\bm{w}\in\mathbb{S}^{D+1}_{r}}\Omega(\bm{w})$ is an open cover of
$\mathbb{S}^{D+1}_{r}$, we can find a finite set
$\\{\bm{w}^{(q)}\\}_{q=1}^{Q}\subset\mathbb{S}^{D+1}_{r}$ such that
$\mathbb{S}^{D+1}_{r}\subset\cup_{q=1}^{Q}\Omega(\bm{w}^{(q)})$. We take
$T=\max_{1\leq q\leq Q}T_{\bm{w}^{q}}$ and get
$J(t\bm{w})>J(\bm{\alpha}^{*})+1$ for all $t>T$ and
$\bm{w}\in\mathbb{S}^{D+1}_{r}$. Therefore, $\\{\bm{\alpha}^{(p)}\\}$ is
bounded in $C_{r}\cap\\{|\bm{\alpha}|\leq T\\}$, which is a contradiction.

The last argument indicates that
$\sqrt{|\alpha_{0}^{(p)}|^{2}+|\hat{\bm{\alpha}}^{(p)}|^{2}}/|\alpha_{D+1}^{(p)}|\to\infty$,
which implies that either $|\alpha_{0}^{(p)}|/|\alpha_{D+1}^{(p)}|$ or
$|\hat{\bm{\alpha}}^{(p)}|/|\alpha_{D+1}^{(p)}|$ is unbounded. Note that we
can assume $|\alpha_{D+1}|\geq\delta>0$ due to Lemma 3.3.

To show that $J(\bm{\alpha}^{(p)})$ is unbounded, we recall Eq. (3.4):

$J(\bm{\alpha})\geq\frac{C}{\sqrt{|\alpha_{D+1}|}}\exp\left(\alpha_{0}+C\frac{|\hat{\bm{\alpha}}|^{2}}{|\alpha_{D+1}|}\right)-\rho(\alpha_{0}+U|\hat{\bm{\alpha}}|).$

Then we follow Ref. [26] and divide the argument into five cases, where we use
$M^{(p)}$ to denote terms that go to $+\infty$ as $p\to\infty$.

  1. (1)

$\\{\alpha_{0}^{(p)}\\}$ is bounded. Then
$|\alpha_{0}^{(p)}|/|\alpha_{D+1}^{(p)}|$ is bounded and we have

$\frac{|\hat{\bm{\alpha}}^{(p)}|}{|\alpha_{D+1}^{(p)}|}\to\infty\quad\text{and}\quad|\hat{\bm{\alpha}}^{(p)}|\geq\frac{\delta|\hat{\bm{\alpha}}^{(p)}|}{|\alpha_{D+1}^{(p)}|}\to\infty.$

Thus it follows that

$\displaystyle J(\bm{\alpha}^{(p)})$
$\displaystyle\geq\frac{C}{\sqrt{|\alpha_{D+1}^{(p)}|}}\exp\left(C\frac{|\hat{\bm{\alpha}}^{(p)}|^{2}}{|\alpha_{D+1}^{(p)}|}\right)-C|\hat{\bm{\alpha}}^{(p)}|-C$
$\displaystyle\geq
C\frac{|\hat{\bm{\alpha}}^{(p)}|^{4}}{|\alpha_{D+1}^{(p)}|^{\frac{5}{2}}}-C|\hat{\bm{\alpha}}^{(p)}|-C\to+\infty.$

  2. (2)

$\alpha_{0}^{(p)}\to+\infty$ and
$|\alpha_{0}^{(p)}|/|\alpha_{D+1}^{(p)}|\to+\infty$. In this case, we have

$\displaystyle J(\bm{\alpha}^{(p)})$
$\displaystyle\geq\frac{C}{\sqrt{|\alpha_{D+1}^{(p)}|}}\left(1+\frac{|\hat{\bm{\alpha}}^{(p)}|^{2}}{|\alpha_{D+1}^{(p)}|}\right)\exp(\alpha_{0}^{(p)})-C\alpha_{0}^{(p)}-C|\hat{\bm{\alpha}}^{(p)}|$
$\displaystyle\geq
M^{(p)}\alpha_{0}^{(p)}+M^{(p)}|\hat{\bm{\alpha}}^{(p)}|^{2}-C\alpha_{0}^{(p)}-C|\hat{\bm{\alpha}}^{(p)}|\to+\infty.$

  3. (3)

$\alpha_{0}^{(p)}\to+\infty$ and $|\alpha_{0}^{(p)}|/|\alpha_{D+1}^{(p)}|$ is
bounded. As in the first case, we have

$\frac{|\hat{\bm{\alpha}}^{(p)}|}{|\alpha_{D+1}^{(p)}|}\to\infty\quad\text{and}\quad|\hat{\bm{\alpha}}^{(p)}|\geq\frac{\delta|\hat{\bm{\alpha}}^{(p)}|}{|\alpha_{D+1}^{(p)}|}\to\infty.$

Thus it follows that

$\displaystyle J(\bm{\alpha}^{(p)})$
$\displaystyle\geq\frac{C}{\sqrt{|\alpha_{D+1}^{(p)}|}}(1+\alpha_{0}^{(p)})\exp\left(C\frac{|\hat{\bm{\alpha}}^{(p)}|^{2}}{|\alpha_{D+1}^{(p)}|}\right)-C\alpha_{0}^{(p)}-C|\hat{\bm{\alpha}}^{(p)}|$
$\displaystyle\geq
M^{(p)}|\hat{\bm{\alpha}}^{(p)}|+M^{(p)}\alpha_{0}^{(p)}-C\alpha_{0}^{(p)}-C|\hat{\bm{\alpha}}^{(p)}|\to+\infty.$

  4. (4)

$\alpha_{0}^{(p)}\to-\infty$ and
$\alpha_{0}^{(p)}+U|\hat{\bm{\alpha}}^{(p)}|\to-\infty$. It is easy to see
that $J(\bm{\alpha}^{(p)})\to+\infty$ in this case.

  5. (5)

$\alpha_{0}^{(p)}\to-\infty$ and
$\alpha_{0}^{(p)}+U|\hat{\bm{\alpha}}^{(p)}|\geq-M_{L}$. Then we have
$|\hat{\bm{\alpha}}^{(p)}|\to+\infty$. Furthermore, we claim
$|\hat{\bm{\alpha}}^{(p)}|/|\alpha_{D+1}^{(p)}|\to\infty$. Otherwise,
$|\alpha_{0}^{(p)}|/|\alpha_{D+1}^{(p)}|$ is unbounded. We can assume
$M_{L}>0$ and get

$\frac{\alpha_{0}^{(p)}}{|\alpha_{D+1}^{(p)}|}+\frac{\hat{\bm{\alpha}}^{(p)}\cdot\bm{U}}{|\alpha_{D+1}^{(p)}|}\geq\frac{-M_{L}}{|\alpha_{D+1}^{(p)}|}\geq\frac{-M_{L}}{\delta}>-\infty.$

Thus $|\alpha_{0}^{(p)}|/|\alpha_{D+1}^{(p)}|\to\infty$ yields
$|\hat{\bm{\alpha}}^{(p)}|/|\alpha_{D+1}^{(p)}|\to\infty$. Now we have

$\displaystyle J(\bm{\alpha}^{(p)})$
$\displaystyle\geq\frac{C}{\sqrt{|\alpha_{D+1}^{(p)}|}}\exp\left(C\frac{|\hat{\bm{\alpha}}^{(p)}|^{2}}{|\alpha_{D+1}^{(p)}|}-C|\hat{\bm{\alpha}}^{(p)}|-M_{L}\right)-C|\hat{\bm{\alpha}}^{(p)}|-C$
$\displaystyle\geq
M^{(p)}\exp\left(\frac{C}{2}\frac{|\hat{\bm{\alpha}}^{(p)}|^{2}}{|\alpha_{D+1}^{(p)}|}-C|\hat{\bm{\alpha}}^{(p)}|\right)-C|\hat{\bm{\alpha}}^{(p)}|-C\to+\infty.$

Consequently, we have shown that $J(\bm{\alpha}^{(p)})$ has no upper bound.
Hence $B_{M}$ is bounded and the proof is completed.

We close this section with the following remark.

######  Remark 3.6.

Recall Eq. (2.12). At the discrete equilibrium $\mathcal{E_{L}}$, we have
$\sigma^{2}<D\theta$. To see this, we deduce from
$\rho\bm{U}=\sum_{m=1}^{N}\rho_{m}u_{m}\bm{l}_{m}$ that

$\rho
U^{2}=\frac{1}{\rho}\left(\rho\bm{U}\cdot\rho\bm{U}\right)=\frac{1}{\rho}\sum_{m,n}\rho_{m}\rho_{n}\left(u_{m}\bm{l}_{m}\cdot
u_{n}\bm{l}_{n}\right).$

Then we have

$\displaystyle\sum_{m}\rho_{m}u_{m}^{2}-\rho U^{2}$
$\displaystyle=\frac{1}{2\rho}\sum_{m,n}\rho_{m}\rho_{n}\left(u_{m}^{2}+u_{n}^{2}-2u_{m}\bm{l}_{m}\cdot
u_{n}\bm{l}_{n}\right)$
$\displaystyle=\frac{1}{2\rho}\sum_{m,n}\rho_{m}\rho_{n}\|u_{m}\bm{l}_{m}-u_{n}\bm{l}_{n}\|^{2}>0.$

Thus, the conclusion follows from the relation $\rho
U^{2}+D\rho\theta=\sum_{m}\rho_{m}u_{m}^{2}+\rho\sigma^{2}$.

##  4\. Spatial-time models

Besides the spatial-time variables $\bm{x}$ and $t$, $\xi\in\mathbb{R}$ is
another continuous variable in the DVDM model Eq. (2.10). In this section, we
derive models only with $\bm{x}$ and $t$ as continuous variables by treating
$\xi$ in two ways.

###  4.1. Discretizing $\xi$

In this first way, the velocity variable $\xi$ in Eq. (2.10) is replaced with
a set of fixed nodes $\xi_{mk}=k\Delta\xi+\xi_{0}$ for
$k=1,\dots,N_{m}^{\prime}$ and $\Delta\xi,\xi_{0}\in\mathbb{R}$. Namely, each
$f_{m}(t,\bm{x},\xi)$ is represented by an $N_{m}^{\prime}$-vector
$\left(f_{mk}(t,\bm{x})\right)_{k=1}^{N_{m}^{\prime}}$. The resultant
governing equation for each $f_{mk}=f_{mk}(t,\bm{x})$ reads as

$\partial_{t}f_{mk}+\xi_{mk}\bm{l}_{m}\cdot\nabla_{\bm{x}}f_{mk}=\frac{1}{\tau}\left(\mathcal{E}_{mk}-f_{mk}\right)$
(4.1)

for $1\leq k\leq N_{m}^{\prime}$ and $1\leq m\leq N$. Note that $\Delta\xi$
may depend on the direction $m$ and the index $k$. For the following reason,
we denote this kind of models as DVD-DVM.

The DVD-DVM in Eq. (4.1) is similar to DVM [24] but allows a new way in
selecting the discrete velocities. Indeed, common DVM practices use discrete-
velocity nodes in a uniform cubic lattice of $\mathbb{R}^{D}$ [14, 24]. By
contrast, the proposed DVD-DVM creates discrete velocities radially
distributed in $\mathbb{R}^{D}$ and centered at the origin of the velocity
space.

To determine the equilibrium $\mathcal{E}_{mk}$ in Eq. (4.1), we define the
macroscopic fluid quantity $\bm{\rho}=(\rho,\rho\bm{U},\rho E)^{T}$ as

$\rho=\sum_{m,k}f_{mk}\Delta\xi,\quad\rho\bm{U}=\sum_{m,k}\xi_{mk}f_{mk}\bm{l}_{m}\Delta\xi,\quad\rho
E=\sum_{m,k}\frac{1}{2}\xi_{mk}^{2}f_{mk}\Delta\xi.$ (4.2)

For a given $\bm{\rho}$, we follow the minimum entropy principle and determine
the discrete equilibrium $\left\\{\mathcal{E}_{mk}\geq 0\right\\}_{1\leq m\leq
N,1\leq k\leq N_{m}^{\prime}}$ so that it minimizes the discrete entropy

$H[\\{g_{mk}\\}]=\sum_{m,k}(g_{mk}\ln g_{mk}-g_{mk})\Delta\xi$ (4.3)

among all possible $\\{g_{mk}\geq 0\\}_{1\leq m\leq N,1\leq k\leq
N_{m}^{\prime}}$ satisfying Eq. (4.2). For the existence of such a minimizer,
we have the following analogue of Theorem 2.1, where
$\bm{m}_{mk}=\left(1,\xi_{mk}\bm{l}_{m},\frac{1}{2}\xi_{mk}^{2}\right)^{T}\in\mathbb{R}^{D+2}$.

######  Theorem 4.1.

If $\bm{\rho}\in\mathbb{R}^{D+2}$ is strictly realizable, i.e., there exists a
set $\\{g_{mk}>0\\}_{1\leq m\leq N,1\leq k\leq N_{m}^{\prime}}$ satisfying Eq.
(4.2), then the minimization problem above has a unique solution
$\\{\mathcal{E}_{mk}\\}$ which can be expressed as

$\mathcal{E}_{mk}=\exp(\bm{\alpha}^{\prime}\cdot\bm{m}_{mk})$ (4.4)

with a vector $\bm{\alpha}^{\prime}\in\mathbb{R}^{D+2}$. Moreover,
$\bm{\alpha}^{\prime}$ minimizes the function

$J(\bm{\alpha}):=\sum_{m,k}\exp(\bm{\alpha}\cdot\bm{m}_{mk})\Delta\xi-\bm{\rho}\cdot\bm{\alpha}.$
(4.5)

This theorem can be proved in a similar fashion as the conventional DVM with
the minimum entropy principle [24]. It provides a simple way to compute
$\\{\mathcal{E}_{mk}\\}\in\mathbb{R}^{\mathcal{N}^{\prime}}$ with
$\mathcal{N}^{\prime}=\sum_{m=1}^{N}N_{m}^{\prime}$ by solving the
$(D+2)$-dimensional minimization problem Eq. (4.5) without constraints.

Besides, our DVDM in Eq. (2.10) offers another possibility to determine
$\\{\mathcal{E}_{mk}\\}$ by conforming ‘more loosely’ to the minimum entropy
principle. Recall Eq. (2.11) that the DVDM equilibrium in each direction
$\bm{l}_{m}$ is a Gaussian distribution parameterized with
$(\rho_{m},u_{m},\sigma^{2})$. We may choose $\\{\mathcal{E}_{mk}\\}$ so that
for each $m$, $\\{\mathcal{E}_{mk}\\}_{1\leq k\leq N_{m}^{\prime}}$ minimizes
the discrete entropy

$H_{m}[\\{g_{mk}\\}]=\sum_{k=1}^{N_{m}^{\prime}}(g_{mk}\ln
g_{mk}-g_{mk})\Delta\xi$ (4.6)

among all possible $\\{g_{mk}\geq 0\\}_{1\leq k\leq N_{m}^{\prime}}$
satisfying

$\rho_{m}=\sum_{k=1}^{N_{m}^{\prime}}g_{mk}\Delta\xi,\quad\rho_{m}u_{m}=\sum_{k=1}^{N_{m}^{\prime}}\xi_{mk}g_{mk}\Delta\xi,\quad\rho_{m}(u_{m}^{2}+\sigma^{2})=\sum_{k=1}^{N_{m}^{\prime}}\xi_{mk}^{2}g_{mk}\Delta\xi.$

For these $N$ minization problems, the analogue of Theorem 4.1 holds. Namely,
each of them has a unique solution $\\{\mathcal{E}_{mk}\\}_{1\leq k\leq
N_{m}^{\prime}}$ which can be expressed as

$\mathcal{E}_{mk}=\exp\left(\bm{\alpha}^{\prime}_{m}\cdot\bm{m}^{\prime}_{mk}\right)$
(4.7)

with a vector $\bm{\alpha}_{m}^{\prime}\in\mathbb{R}^{3}$ and
$\bm{m}^{\prime}_{mk}=\left(1,\xi_{mk},\frac{1}{2}\xi_{mk}^{2}\right)^{T}$.
Moreover, $\bm{\alpha}_{m}^{\prime}$ minimizes the function

$J_{m}(\bm{\alpha}):=\sum_{k=1}^{N_{m}^{\prime}}\exp\left(\bm{\alpha}\cdot\bm{m}^{\prime}_{mk}\right)\Delta\xi-\bm{\rho}_{m}^{\prime}\cdot\bm{\alpha}$

with
$\bm{\rho}_{m}^{\prime}=\rho_{m}\left(1,u_{m},\frac{u_{m}^{2}+\sigma^{2}}{2}\right)^{T}\in\mathbb{R}^{3}$.
We note that this approach does not necessarily find a minimum of the ‘total’
discrete entropy Eq. (4.3).

In summary, we have derived two kinds of DVD-DVM models Eq. (4.1):

  1. (1)

(DVD-DVM-I). $\mathcal{E}_{mk}$ is given with Eq. (4.4) by minimizing the
discrete entropy Eq. (4.3).

  2. (2)

(DVD-DVM-II). $\mathcal{E}_{mk}$ is given with Eq. (4.7) by minimizing the
discrete entropy Eq. (4.6) in each direction. This approach requires
$(\rho_{m},u_{m},\sigma^{2})_{m=1}^{N}$ defined in Eqs. (2.11 & 2.12) to be
computed by DVDM as a first step.

###  4.2. Gaussian-EQMOM

Besides the DVD-DVM, the method of moment can be applied to the DVDM Eq.
(2.10). For this purpose, we define the $k$th velocity moment of
$f_{m}(t,\bm{x},\xi)$ as

$M_{m,k}=M_{m,k}(t,\bm{x})=\int_{\mathbb{R}}\xi^{k}f_{m}(t,\bm{x},\xi)d\xi$
(4.8)

for $k\in\mathbb{N}$. The evolution of $M_{m,k}$ can be derived from Eq.
(2.10) as

$\partial_{t}M_{m,k}+\bm{l}_{m}\cdot\nabla_{\bm{x}}M_{m,k+1}=\frac{1}{\tau}\left(\rho_{m}\Delta_{k}(u_{m},\sigma^{2})-M_{m,k}\right),$
(4.9)

where $\Delta_{k}(u,\sigma^{2})$ denotes the $k$th moment of the normalized
Gaussian function centered at $u$ with a variance $\sigma^{2}$. Meanwhile, the
macroscopic quantity $\bm{\rho}$ in Eq. (2.6) is computed with

$\rho=\sum_{m=1}^{N}M_{m,0},\quad\rho\bm{U}=\sum_{m=1}^{N}\bm{l}_{m}M_{m,1},\quad\rho
E=\sum_{m=1}^{N}\frac{M_{m,2}}{2}.$ (4.10)

There are infinitely many equations in Eq. (4.9). For each $m$, the first
$N^{\prime}$ equations for moments $M_{m,0},\dots,M_{m,N^{\prime}-1}$ are not
closed because the $M_{m,N^{\prime}-1}$-equation contains $M_{m,N^{\prime}}$
in the convection term. Hence a closure method is needed.

The Gaussian-EQMOM assumes that the 1-D distribution $f_{m}(\xi)$ is a sum of
$N_{m}^{\prime}$ Gaussian functions [23]:

$f_{m}(\xi)=\sum_{\alpha=1}^{N_{m}^{\prime}}\frac{w_{m,\alpha}}{\sqrt{2\pi\vartheta_{m}}}\exp\left(-\frac{(\xi-
v_{m,\alpha})^{2}}{2\vartheta_{m}}\right).$ (4.11)

Note that the variance $\vartheta_{m}>0$ is shared among $N_{m}^{\prime}$
kernels (and is thus not dependent on $\alpha$). The weights $w_{m,\alpha}$,
nodes $v_{m,\alpha}$ ($\alpha=1,\dots,N_{m}^{\prime}$) and variance
$\vartheta_{m}$ are unknown parameters to be solved from the
($2N_{m}^{\prime}+1$) nonlinear equations

$M_{m,k}=\sum_{\alpha=1}^{N_{m}^{\prime}}w_{m,\alpha}\Delta_{k}(v_{m,\alpha},\vartheta_{m})\quad\text{for
}k=0,\dots,2N_{m}^{\prime}.$ (4.12)

An algorithm to solve these equations can be found in the literature [7, 23].
Thus, any unknown moments and/or integrals associated with the distribution
can be computed by Eq. (4.11). For example, the convection term
$M_{m,2N_{m}^{\prime}+1}$ is reconstructed as

$\bar{M}_{m,2N_{m}^{\prime}+1}=\sum_{\alpha=1}^{N_{m}^{\prime}}w_{m,\alpha}\Delta_{2N_{m}^{\prime}+1}(v_{m,\alpha},\vartheta_{m}).$
(4.13)

Consequently, we get

$\partial_{t}M_{m}+A_{m}(M_{m})\bm{l}_{m}\cdot\nabla_{\bm{x}}M_{m}=\frac{1}{\tau}(M_{\mathcal{E}m}-M_{m}),$

where

$\displaystyle M_{m}$
$\displaystyle=(M_{m,0},\dots,M_{m,2N_{m}^{\prime}})^{T}\in\mathbb{R}^{2N^{\prime}_{m}+1},$
$\displaystyle M_{\mathcal{E}m}$
$\displaystyle=\rho_{m}\left(\Delta_{0}(u_{m},\sigma^{2}),\dots,\Delta_{2N_{m}^{\prime}}(u_{m},\sigma^{2})\right)^{T}\in\mathbb{R}^{2N^{\prime}_{m}+1},$

and

$A_{m}(M_{m})=\begin{bmatrix}0&1&&&\\\ &0&1&&\\\ &&\ddots&\ddots&\\\ &&&0&1\\\
a_{m,0}&a_{m,1}&\cdots&a_{m,2N^{\prime}_{m}-1}&a_{m,2N^{\prime}_{m}}\end{bmatrix}\in\mathbb{R}^{(2N_{m}^{\prime}+1)\times(2N_{m}^{\prime}+1)}$
(4.14)

with $a_{m,k}=\diffp{\bar{M}_{m,2N_{m}^{\prime}+1}}{M_{m,k}}$ for
$k=0,\dots,2N_{m}^{\prime}$.

Putting the last equations together, we derive the following spatial-time
model

$\partial_{t}\mathcal{M}+\sum_{i=1}^{D_{x}}A^{(i)}\partial_{x_{i}}\mathcal{M}=\frac{1}{\tau}(M_{\mathcal{E}}-\mathcal{M}).$
(4.15)

Here

$\mathcal{M}=\left(M_{1}^{T},\dots,M_{N}^{T}\right)^{T}\in\mathbb{R}^{\mathcal{N}^{\prime}},\quad
M_{\mathcal{E}}=\left(M_{\mathcal{E}1}^{T},\dots,M_{\mathcal{E}N}^{T}\right)^{T}\in\mathbb{R}^{\mathcal{N}^{\prime}}$

with $\mathcal{N}^{\prime}=\sum_{m=1}^{N}(2N_{m}^{\prime}+1)$ and

$A^{(i)}=\text{diag}\left\\{(\bm{l}_{1}\cdot\bm{e}_{i})A_{1}(M_{1}),\dots,(\bm{l}_{N}\cdot\bm{e}_{i})A_{N}(M_{N})\right\\}\in\mathbb{R}^{\mathcal{N}^{\prime}\times\mathcal{N}^{\prime}},$

where $\bm{e}_{i}$ is the $i$-th column of the identity matrix of order
$D_{x}$.

######  Theorem 4.2.

The moment system (4.15) is hyperbolic.

###### Proof.

According to Ref. [18], each matrix $A_{m}(M_{m})$ defined in Eq.(4.14) is
strictly hyperbolic. Namely, it has $(2N_{m}^{\prime}+1)$ distinct real
eigenvalues. Therefore, for any
$\bm{y}=(y_{1},...,y_{D_{x}})^{T}\in\mathbb{R}^{D_{x}}$ the block-diagonal
matrix $\sum_{i}y_{i}A^{(i)}$ is real diagonalizable. ∎

######  Remark 4.3.

This approach, denoted as DVD-EQMOM, leads to a convenient multidimensional
version of quadrature-based method of moments, which seems better understood
than those in Refs. [7, 23].

##  5\. Numerical schemes

In this section, we show that the DVDM provides new numerical solvers to
simulate rarefied flows. The solvers consist of two parts: spatial-time
discretization and computation of equilibrium. As the first step of this
project, we present the solvers only for spatially one-dimensional models,
while the velocity is still in $\mathbb{R}^{D}$.

###  5.1. Spatial-time discretization

In this part, we discretize the spatial-time models constructed in the
previous section. To do this, we denote by $g_{i}^{n}$ an approximation of
function $g(t,x)$ over the grid-block
$]x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}[\times[t_{n},t_{n+1})$ with
$x_{i}=i\Delta x$, $t_{n}=n\Delta t$, $i=0,\pm 1,\pm 2,\cdots$, and
$n=0,1,2,\cdots$.

####  5.1.1. DVD-DVM

For this kind of models, the governing Eq. (4.1) is approximated by the
following implicit-explicit scheme

$\displaystyle f_{mk,i}^{n+1}=f_{mk,i}^{n}-\frac{\Delta t}{\Delta
x}\left(\mathcal{F}^{n}_{mk,i+\frac{1}{2}}-\mathcal{F}^{n}_{mk,i-\frac{1}{2}}\right)+\frac{\Delta
t}{\tau}\left(\mathcal{E}^{n}_{mk,i}-f^{n+1}_{mk,i}\right).$ (5.1)

The numerical flux is taken as

$\displaystyle\mathcal{F}^{n}_{mk,i+\frac{1}{2}}=\frac{1}{2}\Big{(}$
$\displaystyle\xi_{mk}\bm{l}_{m}\cdot\bm{e}_{1}\left(f^{n}_{mk,i+1}+f^{n}_{mk,i}\right)$
(5.2)
$\displaystyle-|\xi_{mk}\bm{l}_{m}\cdot\bm{e}_{1}|\left(f^{n}_{mk,i+1}-f^{n}_{mk,i}-\Phi^{n}_{mk,i+\frac{1}{2}}\right)\Big{)}.$

The flux limiter $\Phi^{n}_{mk,i+\frac{1}{2}}$ allows to obtain a second-order
scheme, and $\Phi^{n}_{mk,i+\frac{1}{2}}=0$ reduces to a first-order scheme.
The equilibrium state $\mathcal{E}^{n}_{mk,i}$ is solved from the local fluid
quantity $\bm{\rho}^{n}_{i}$ by either Eq. (4.4) (DVD-DVM-I) or Eq. (4.7)
(DVD-DVM-II).

Note that the source in Eq. (5.1) is semi-implicit because $f^{n+1}_{mk,i}$ is
used to make the scheme stable, but $f_{mk,i}^{n+1}$ can be explicitly
obtained.

####  5.1.2. DVD-EQMOM

For the momemt system, the governing Eq. (4.9) is approximated by the
implicit-explicit scheme:

$\displaystyle M_{m,k,i}^{n+1}=$ $\displaystyle M_{m,k,i}^{n}-\frac{\Delta
t}{\Delta
x}\bm{l}_{m}\cdot\bm{e}_{1}\left(\mathcal{F}^{n}_{m,k+1,i+\frac{1}{2}}-\mathcal{F}^{n}_{m,k+1,i-\frac{1}{2}}\right)$
(5.3) $\displaystyle+\frac{\Delta
t}{\tau}\left(M^{n}_{\mathcal{E}m,k,i}-M_{m,k,i}^{n+1}\right).$

The flux is modeled by a ‘kinetic-based’ definition [7, 23]. For example,

$\displaystyle\mathcal{F}^{n}_{m,k+1,i+\frac{1}{2}}=\frac{1}{2}\int_{\mathbb{R}}\xi^{k}\Big{[}$
$\displaystyle\xi\left(f^{n}_{m,i+1}+f^{n}_{m,i}\right)$ (5.4)
$\displaystyle-\text{sgn}(\bm{l}_{m}\cdot\bm{e}_{1})|\xi|\left(f^{n}_{m,i+1}-f^{n}_{m,i}-\Phi^{n}_{m,i+\frac{1}{2}}\right)\Big{]}d\xi.$

The limiter $\Phi^{n}_{m,i+\frac{1}{2}}$ allows to obtain a second-order
scheme, and $\Phi^{n}_{m,i+\frac{1}{2}}=0$ reduces to a first-order scheme.
Then we have

$\mathcal{F}^{n}_{m,k+1,i+\frac{1}{2}}=\left\\{\begin{aligned}
\int_{0}^{\infty}\xi^{k+1}f_{m,i}^{n}d\xi+\int_{-\infty}^{0}\xi^{k+1}f_{m,i+1}^{n}d\xi,\quad&\text{if
}\bm{l}_{m}\cdot\bm{e}_{1}>0,\\\
\int_{0}^{\infty}\xi^{k+1}f_{m,i+1}^{n}d\xi+\int_{-\infty}^{0}\xi^{k+1}f_{m,i}^{n}d\xi,\quad&\text{if
}\bm{l}_{m}\cdot\bm{e}_{1}<0.\end{aligned}\right.$ (5.5)

To evaluate the integrals above, we define

$\langle\xi^{k+1}\rangle^{\pm}(v,\vartheta)=\int_{\xi\gtrless
0}\xi^{k+1}\frac{1}{\sqrt{2\pi\vartheta}}\exp\left(-\frac{(\xi-v)^{2}}{2\vartheta}\right)d\xi,$

which can be computed analytically with recursive relations in terms of $k$
(Eq. (B.1) in Ref. [7]). Thus, due to the ansatz Eq. (4.11), the integrals in
Eq. (5.5) become

$\int_{\xi\gtrless
0}\xi^{k+1}f_{m,i}d\xi=\sum_{\alpha=1}^{N_{m}^{\prime}}w_{m,\alpha,i}\langle\xi^{k+1}\rangle^{\pm}\left(v_{m,\alpha,i},\vartheta_{m,i}\right).$

The moment inversion algorithm in Refs. [7, 23] is used to solve
$w_{m,\alpha,i}$, $v_{m,\alpha,i}$, and $\vartheta_{m,i}$
($\alpha=1,\dots,N_{m}^{\prime}$) from the moments
$M_{m,0,i},\dots,M_{m,2N_{m}^{\prime},i}$.

In the scheme Eq. (5.3), the moments correspond to the equilibrium state

$M^{n}_{\mathcal{E}m,k,i}=\rho^{n}_{m,i}\Delta_{k}\left(u^{n}_{m,i},(\sigma^{2})^{n}_{i}\right).$

The analytical expression of $\Delta_{k}(u,\sigma^{2})$ is given in Ref. [18].
The equilibrium state parameters $\rho^{n}_{m,i},\ u^{n}_{m,i},\
(\sigma^{2})^{n}_{i}$ are solved by the local fluid quantity
$\bm{\rho}^{n}_{i}$, as detailed in Section 5.2.

###  5.2. Computation of equilibrium

This subsection presents details of computing the discrete equilibrium
$\mathcal{E_{L}}$ in Eq. (2.11) for a given macroscopic quantity $\bm{\rho}$.
It can be done by minimizing $J(\bm{\alpha})$ in Eq. (2.16), which is strictly
convex. In this work we use the gradient descent method [10] to solve
$\mathcal{E_{L}}$. Future work is needed to develop stable and efficient
higher-order optimization algorithms for this purpose.

The algorithm to minimize $J(\bm{\alpha})$ is stated below.

  1. (1)

Let the initial value
$\bm{\alpha}^{(0)}\in\mathbb{R}^{D+1}\times\mathbb{R}^{-}$ be given. It is
found that the simple choice of $\bm{\alpha}_{eq}$ of the continuous
equilibrium in Eq. (2.1), even with $U=0$, works reasonably well in our test
case.

  2. (2)

Repeat.

    1. (i)

Update $\bm{\alpha}$ by

$\bm{\alpha}^{(n+1)}=\bm{\alpha}^{(n)}-l_{r}\nabla_{\bm{\alpha}}J(\bm{\alpha}^{(n)}).$

The gradient $\nabla_{\bm{\alpha}}J$ has the explicit form as:

$\nabla_{\bm{\alpha}}J=\left(\sum_{m=1}^{N}\rho_{m},\
\sum_{m=1}^{N}\rho_{m}u_{m}\bm{l}_{m},\
\sum_{m=1}^{N}\frac{\rho_{m}(u_{m}^{2}+\sigma^{2})}{2}\right)^{T}-\bm{\rho}.$

Here $\rho_{m},u_{m},\sigma$ are defined in Eq. (2.12).

    2. (ii)

In each iteration, the learning rate $l_{r}\in(0,1)$ is backtracked based on
the Armijo-Goldstein condition [2] and the constraint that
$\alpha_{D+1}^{(n+1)}<0$.

    3. (iii)

The recovered fluid quantity from $\bm{\alpha}$, i.e.,
$\langle\bm{m}_{\mathcal{L}}\exp\left(\bm{\alpha}\cdot\bm{m}_{\mathcal{L}}\right)\rangle_{\mathcal{L}}$,
is slightly different from $\bm{\rho}$ due to the round-off error. To avoid
error accumulation, we re-assign the results
$\rho_{m}^{(n+1)}\leftarrow\rho_{m}^{(n+1)}\left(\rho/\sum_{m}\rho_{m}^{(n+1)}\right)$,
or equivalently,
$\alpha_{0}^{(n+1)}\leftarrow\alpha_{0}^{(n+1)}+\ln\left(\rho/\sum_{m}\rho_{m}^{(n+1)}\right)$.
Here $\alpha_{0}$ is the first component of $\bm{\alpha}$.

  3. (3)

Until: $\|\nabla_{\bm{\alpha}}J(\bm{\alpha}^{(n)})\|_{2}<\epsilon$.

In this work we take $\epsilon=10^{-8}$. The performance of the algorithm is
illustrated in Fig. 1 with different numbers of velocity direction $N$. It is
shown that the stopping criterion can be achieved with no more than 160 steps
of iteration, almost regardless of $N$ from 5 to 500. Thus, the computation of
discrete equilibrium in DVDM is numerically efficient.

Figure 1. Performance of the $\mathcal{E_{L}}$-algorithm for a given
$\bm{\rho}\in\mathbb{R}^{4}$: $\rho=1,\ \bm{U}=(0.4,0.8)^{T}\text{ and
}\theta=1$ with different numbers of velocity directions $N$ from 5 to 500.
The directions are chosen to be $\bm{l}_{m}=(\cos\gamma_{m},\sin\gamma_{m})$
with $\gamma_{m}=\frac{m-1}{N}\pi$ for $m=1,\dots,N$.

Solving the discrete equilibrium $\mathcal{E_{L}}$ is necessary for both DVD-
DVM-II and DVD-EQMOM. Nevertheless, for DVD-DVM-I, the discrete equilibrium is
solved by minimizing $J(\bm{\alpha})$ in Eq. (4.5). The difference between
DVD-DVM-I and II is studied quantitatively in Fig. 2. We see that for a given
$\bm{\rho}\in\mathbb{R}^{4}$ and a fixed set of directions, DVD-DVM-I
generally results in smaller entropy than II, but the difference almost
disappears when the number of velocity nodes in each direction
$N_{m}^{\prime}>10$ while keeping the node distance $\Delta\xi$ unchanged
(indicating that the nodes cover a greater range). Thus, both DVD-DVM methods
behave similarly in this case.

Figure 2. Discrete entropies of DVD-DVM-I & II (left y-axis) and their
difference (right y-axis) for a given $\bm{\rho}\in\mathbb{R}^{4}$: $\rho=1$,
$\bm{U}=(0.6,1.0)^{T}$, and $\theta=1.5$. The number of directions $N=8$, and
the directions are chosen to be $\bm{l}_{m}=(\cos\gamma_{m},\sin\gamma_{m})$
with $\gamma_{m}=\frac{m-1}{8}\pi$ for $m=1,\dots,8$. The number of velocity
nodes in each direction, $N_{m}^{\prime}$ as abscissa, ranges from 6 to 36.
The velocity nodes are $\xi_{mk}=\pm(k-0.5)$ for $k=1,\dots,N_{m}^{\prime}/2$.

######  Remark 5.1 (Temperature influence).

For $D=2$, Eq. (A.3) specifies a lower bound for the realizable temperature
$\theta$. It is found that temperature has a considerable influence on the
number of iterations. As revealed in Table 1, for given $\rho$ and $\bm{U}$,
lower values of $\theta$ leads to substantially more iteration and longer CPU
time.

Table 1. CPU time for the $J(\bm{\alpha})$-algorithm with $\rho=1$,
$\bm{U}=(0.8,0.4)^{T}$ and different $\theta$

$\theta$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 | 1.1 | 1.2
---|---|---|---|---|---|---|---|---|---|---|---|---
Time (ms) | 14.7 | 6.32 | 4.31 | 2.41 | 2.15 | 1.91 | 1.59 | 1.82 | 1.45 | 1.33 | 1.14 | 1.20

  * •

The CPU time was recorded by running Matlab code on Intel(R) Core(TM)
i7-1065G7. Here, 15 directions are chosen to be
$\bm{l}_{m}=(\cos\gamma_{m},\sin\gamma_{m})$ with
$\gamma_{m}=\frac{m-1}{15}\pi$ ($m=1,\dots,15$) and, according to Eq. (A.3),
the lower bound of $\theta$ is 0.0036.

##  6\. Numerical experiments

Numerical results for two benchmark flow problems, the Couette flow and the
1-D Riemann problem, are reported in this section to show the performance of
the DVDM-BGK models. We assume the velocity space to be two-dimensional,
namely, $D=2$, $\bm{\xi}=(\xi_{x},\xi_{y})^{T}$ and
$f=f(t,\bm{x},\xi_{x},\xi_{y})$. Both DVD-EQMOM (Section 4.2) and DVD-DVM
(Section 4.1) are tested, together with the first-order schemes detailed in
Section 5.

###  6.1. Couette flow

The first problem is the planar Couette flow between two infinite parallel
walls located at $x=\pm H$ with a distance $L=2H$. The left and right walls
move with constant velocities $\pm v_{w}\bm{e}_{y}$. We set $H=0.5$ and
$v_{w}=0.1$. The two walls drive the fluid between them from rest to a final
steady state. This implies that the distribution $f$ is only dependent on $x$.
We take the wall temperature $\theta_{w}=1$, and the Mach number
$Ma=\frac{v_{w}}{\sqrt{\frac{D+2}{D}\theta}}=\frac{1}{10\sqrt{2}}\approx
0.0707$, representing a micro-Couette flow.

For the DVDM simulation, $N$ velocity directions are chosen to be
$\bm{l}_{m}=(\cos\gamma_{m},\sin\gamma_{m})$ with
$\gamma_{m}=\frac{m-1}{N}\pi$ ($m=1,...,N$). The initial temperature and
velocity are taken to be 1 and $\bm{0}$, respectively. The 1-D computational
domain $-0.5<x<0.5$ is discretized into 200 uniform cells with $\Delta
x=0.005$. The time step ensures that the CFL number is less than 0.5. While
the boundary velocity and temperature are kept constant, the density at the
boundaries are adjusted at every moment to ensure zero mass fluxes at the
boundary.

A number of simulations for different flow regimes $\kappa=0.1,\ 1,\text{ and
}10$ are conducted, where $\kappa=(\sqrt{\pi}/2)\rm Kn$ and the Knudsen number
Kn is defined as [15]

${\rm Kn}=\frac{\lambda}{L}=\frac{\tau}{L}\sqrt{\frac{\pi\theta}{2}}.$

Here $\lambda$ denotes the mean free path of the molecules. Thus, varying
$\tau$ in the source term of BGK equation gives the designed values of
$\kappa$.

Fig. 3 shows the dimensionless velocity profiles for various $\kappa$ using
two-node DVD-EQMOM (namely, $N_{m}^{\prime}=2$ in Eq. (4.11)). Only the right
side ($x>0$) is plotted due to symmetry. Increasing $\kappa$ from 0.1 to 10
(namely, increasing $\tau$) shifts the flow from the hydrodynamic regime (with
velocity slip) to a free-molecular regime, and thus the flow is less efficient
in following the moving boundaries. This property is successfully captured by
all involved numbers of directions $N=5,\ 15,\text{ and }45$. For $N=5$, the
DVD-EQMOM velocities for $\kappa=1.0\text{ and }10$ deviate from the DSMC
results to a greater extent than that for $\kappa=0.1$. Increasing $N$ is
effective in reducing the numerical errors especially for $\kappa=1.0\text{
and }10$.

Figure 3. Velocity profiles of the Couette flow for various Knudsen numbers
$\kappa=(\sqrt{\pi}/2){\rm Kn}$ obtained with the two-node DVD-EQMOM. The
number of velocity direction is varied as (a) $N=5$, (b) $N=15$, and (c)
$N=45$. The directions are chosen to be
$\bm{l}_{m}=(\cos\gamma_{m},\sin\gamma_{m})$ with
$\gamma_{m}=\frac{m-1}{N}\pi$ for $m=1,\dots,N$. The DSMC results are
extracted from Ref. [15].

Fig. 4 presents the dimensionless velocity profiles for various $\kappa$ using
both DVD-DVM-I and DVD-DVM-II with 15 prescribed velocity directions. In each
direction 26 discrete velocities are selected as $\xi_{mk}=\pm(0.2k-0.1)$ for
$k=1,\dots,13$. It is seen that both approaches yield very accurate velocity
profiles as compared with the DSMC results. In particular, the nonlinearity of
the velocity profiles near the wall is correctly reproduced.

Figure 4. Velocity profiles of the Couette flow for various Knudsen numbers
$\kappa=(\sqrt{\pi}/2){\rm Kn}$ obtained with (a) DVD-DVM-I and (b) DVD-DVM-
II. The number of directions $N=15$, and the directions are chosen to be
$\bm{l}_{m}=(\cos\gamma_{m},\sin\gamma_{m})$ with
$\gamma_{m}=\frac{m-1}{15}\pi$ for $m=1,\dots,15$. The velocity nodes in each
direction are $\xi_{mk}=\pm(0.2k-0.1)$ for $k=1,\dots,13$. The DSMC results
are extracted from Ref. [15].

Fig. 5 further shows the shear stress

$\tau_{xy}=\langle(\xi\cos\gamma_{\mathcal{L}}-u)(\xi\sin\gamma_{\mathcal{L}}-v)f_{\mathcal{L}}\rangle_{\mathcal{L}},\quad(u,v)^{T}=\bm{U},$

normalized by the free-molecular stress $\tau_{\infty}$ for a variety of
Knudsen numbers. The results (dots) are obtained with the two-node DVD-EQMOM
and DVD-DVM-II with 15 prescribed velocity directions. The shear stress
increases with Kn. The two-node DVD-EQMOM seems to overestimate the shear
stress for a wide range of Kn; By contrast, the DVD-DVM-II solutions agree
well with the DSMC results in the whole flow regimes.

Figure 5. Normalized stress of the Couette flow at different Knudsen numbers.
The two-node DVD-EQMOM and DVD-DVM-II are both solved with 15 discrete
directions $\bm{l}_{m}=(\cos\gamma_{m},\sin\gamma_{m})$ and
$\gamma_{m}=\frac{m-1}{15}\pi$ for $m=1,\dots,15$. For the DVD-DVM-II, the
velocity nodes in each direction are $\xi_{mk}=\pm(0.2k-0.1)$ for
$k=1,\dots,13$. The DSMC results are extracted from Ref. [15].

###  6.2. 1-D Riemann problems

Two 1-D Riemann problems from Refs. [7, 12] are solved in this subsection. The
initial macroscopic data are

  1. (1)

Problem (A): $\rho(0,x)=1$ and $\theta(0,x)=\frac{1}{6}$ for all
$x\in\mathbb{R}$; $U(0,x<0)=(1,0)^{T}$ and $U(0,x>0)=(-1,0)^{T}$.

  2. (2)

Problem (B): $\rho(0,x<0)=3.093$ and $\rho(0,x>0)=1$; $U(0,x)=(0,0)^{T}$ and
$\theta(0,x)=1$ for all $x\in\mathbb{R}$.

The corresponding initial distribution $f=f(t,x,\xi_{x},\xi_{y})$ are taken to
be in equilibrium determined by the macroscopic data above.

As for the collision term, two limiting cases are considered: (i)
$\tau=\infty$ and thus the collision term vanishes. In this case, the
distribution $f$ is a traveling wave and its analytical expression is given in
Ref. [7]; (ii) $\tau=0$ corresponds to the classical Euler equation for an
inviscid compressible fluid. Notice that the specific heat ratio in the
analytical expression [32] should be taken as 2 due to the 2-D velocity space
assumption.

In the DVDM simulation, 15 velocity directions are chosen to be
$\bm{l}_{m}=(\cos\gamma_{m},\sin\gamma_{m})$ with
$\gamma_{m}=\frac{m-1}{15}\pi$ for $m=1,\dots,15$. For the DVD-EQMOM, we set
$N_{m}^{\prime}=2$; for the DVD-DVM-II, the velocity nodes in each direction
are chosen as $\xi_{mk}=\pm(0.2k-0.1)$ for $k=1,\dots,13$. The computational
domain is taken to be $-2<x<2$ and is discretized into 400 uniform cells with
$\Delta x=0.01$. The time step ensures that the CFL number is less than 0.5.
For $\tau=0$, the collisions reset the variables (either discretized
distributions or moments) to be solved as the equilibrium states at the end of
each time step.

Figs. 6 & 7 exhibit the spatial distribution of macroscopic quantities at
$t=0.2$ for Problem (A) with $\tau=\infty$ (no collision) and $\tau=0$
(ultrafast collision), respectively. For $\tau=\infty$, the DVD-DVM-II results
agree reasonably well with the analytic solution, while the two-node EQMOM is
less accurate in the region of $-0.5<x<0.5$. In particular, the density at
$x=0$ is underestimated, and the simulated energy (or temperature) at $x=0$ is
greater than the true value. For $\tau=0$, both models well capture the shock
waves propagating towards both sides, but there are some inaccuracies in the
central region around $x=0$ (the contact surface). The inaccuracies may be
attributed to either the DVDM assumptions or numerical schemes. It is
interesting to remark that such a deviation around $x=0$ was also observed in
the two-node EQMOM solutions in Refs. [7, 18].

Figure 6. Macroscopic quantities at $t=0.2$ for the Riemann problem (A) with
$\tau=\infty$ (i.e., no collision).

Figure 7. Macroscopic quantities at $t=0.2$ for the Riemann problem (A) with
$\tau=0$ (i.e., ultrafast collision).

Figs. 8 & 9 illustrate the spatial distribution of macroscopic quantities at
$t=0.2$ for Problem (B) with $\tau=\infty$ (no collision) and $\tau=0$
(ultrafast collision), respectively. As in Problem (A), for $\tau=\infty$ the
DVD-DVM is sufficiently accurate in this free-moving situation, while the two-
node EQMOM introduces non-negligible error for the energy (or temperature)
profile. By contrast, in the Euler limit $\tau=0$, both models exhibit
correctly the rarefaction wave, contact discontinuity and shock wave. However,
the discontinuities are less sharper than the theoretical solution, especially
for $\rho$ and $u$, which may be mainly caused by the first-order schemes. In
comparision with the DVD-DVM, the two-node DVD-EQMOM is less accurate in
describing the right-moving shock wave located at $x\approx 0.35$.

Figure 8. Macroscopic quantities at $t=0.2$ for the Riemann problem (B) with
$\tau=\infty$ (i.e., no collision).

Figure 9. Macroscopic quantities at $t=0.2$ for the Riemann problem (B) with
$\tau=0$ (i.e., ultrafast collision).

######  Remark 6.1.

In solving Problem (B) with $\tau=0$ and 15 directions (results given in Fig.
9), the CPU times of the 2-node DVD-EQMOM and DVD-DVM-II (26 nodes in each
direction) are 120 s and 500 s, respectively. The CPU times were recorded by
running Matlab codes on Intel(R) Core(TM) i7-1065G7.

##  7\. Conclusions and Perspectives

In this paper, we present a novel discrete-velocity-direction model (DVDM)
with a minimum entropy principle. As a semi-continuous discretization model of
BGK-type, it assumes that the molecule velocity has a few prescribed
directions, but the velocity modulus is still continuous. The Maxwellian
equilibrium is defined as a minimizer of a discrete entropy subject to
conservation laws of density, momentum and energy. We show that the discrete
equilibrium is uniquely determined by macroscopic fluid quantities computed
with nonnegative distributions generating finite density and energy, implying
that the DVDM is well defined. Numerically, the discrete equilibrium with a
given macroscopic flow state can be computed efficiently by convex
optimization algorithms. Moreover, the model ensures positivity of the
solutions and has a proper version of $H$-theorem.

The proposed model has several advantages. Firstly, it provides a new way in
choosing discrete velocities for the computational practice of the
conventional discrete-velocity methodology. In this way, we replace the
continuous one-dimensional velocity moduli of DVDM with finite nodes in each
direction and derive new discrete-velocity models (called DVD-DVM).
Furthermore, the 1-D EQMOM can be implemented in each prescribed discrete
direction to derive multidimensional hyperbolic EQMOM conveniently. As
reported in the literature [7, 23], 1-D EQMOM is quite successful in a number
of applications, and it also has good mathematical properties [18]. However, a
satisfactory multidimensional version of EQMOM is not available for a long
time. In this sense, we have offered a solution to the multidimensional
extension of the 1-D EQMOM. To show the performance of the above models, we
simulate two benchmark flows, the Couette flow and 1-D Riemann problems, with
reasonable results at moderate computational costs.

Potential focuses of future work can be directed to: (i) studying the DVDM
mathematically, including the realizable condition and the structural
stability condition [18, 33]; (ii) simulating real-world flows with 3-D
velocity space and more complex boundary conditions with higher-order
numerical schemes.

## Acknowledgment

This work is supported by the National Natural Science Foundation of China
(Grant no. 51906122 and 12071246) and National Key Research and Development
Program of China (Grant no. 2021YFA0719200). The authors acknowledge Mr.
Jialiang Zhou, Mr. Qiqi Rao and Prof. Shuiqing Li at Tsinghua University for
valuable discussions.

##  Appendix A Realizabily of the DVDM

In this appendix, we show that for the two-dimensional DVDM, more macroscopic
states can be realized with more discrete-velocity directions. In this case,
the discrete-velocity directions can be expressed as
$\bm{l}_{m}=(\cos\gamma_{m},\sin\gamma_{m})^{T}$ with $\gamma_{m}\in[0,\pi)$
and therefore can be understood as the complex numbers $e^{i\gamma_{m}}$.
Without loss of generality, we assume $0\leq\gamma_{1}<\dots<\gamma_{N}<\pi$.
The main result here is the following lemma.

######  Lemma A.1.

For any $N$-tuple $(a_{m})_{m=1}^{N}\in\mathbb{R}$ such that

$\sum_{m=1}^{N}a_{m}e^{i\gamma_{m}}=1,$

we have

$\sum_{m=1}^{N}|a_{m}|\geq\frac{\sin\gamma_{1}+\sin\gamma_{N}}{\sin(\gamma_{N}-\gamma_{1})}.$
(A.1)

The minimum is attained when

$a_{1}=\frac{\sin\gamma_{N}}{\sin(\gamma_{N}-\gamma_{1})},\
a_{N}=\frac{-\sin\gamma_{1}}{\sin(\gamma_{N}-\gamma_{1})},\text{ and
}a_{2}=\dots=a_{N-1}=0.$

###### Proof.

Set $z_{m}=a_{m}e^{i\gamma_{m}}$. According to the assumption, we have
$\sum_{m}z_{m}=1$ and

$\sum_{m}|a_{m}|=\sum_{m}|z_{m}|\geq\left\lvert\sum_{m}z_{m}\right\rvert=1.$

Therefore the lemma obviously holds with $\gamma_{1}=0$.

For $\gamma_{1}>0$, since $\gamma_{1}\leq\gamma_{m}<\pi$ for all $m$, $z_{m}$
is not a real number unless $z_{m}=0$. Thus, from $\sum_{m}z_{m}=1$ it follows
that there must be some $z_{m}$ with positive imaginary part (denote their sum
by $z^{+}$) and some $z_{m}$ with negative imaginary part (denote their sum by
$z^{-}$). Obviously, we have $\gamma_{+}:=\arg z^{+}\in(0,\pi)$ and
$\gamma_{-}:=\arg z^{-}\in(-\pi,0)$. Note that $z^{+}+z^{-}=1$, meaning that
the three complex numbers $z^{+}$, $z^{-}$ and $1$ constitute a triangle
(after a proper shift of $z^{-}$) in the complex plane. Then the angle
$(\gamma_{+}-\gamma_{-})$ from $z^{-}$ to $z^{+}$ is in $(0,\pi)$. We thus see
from the law of sines:

$\frac{1}{\sin(\gamma_{+}-\gamma_{-})}=\frac{|z^{+}|}{\sin(-\gamma_{-})}=\frac{|z^{-}|}{\sin\gamma_{+}}$

that

$\sum_{m}|a_{m}|=\sum_{m}|z_{m}|\geq|z^{+}|+|z^{-}|=\frac{\sin\gamma_{+}-\sin\gamma_{-}}{\sin(\gamma_{+}-\gamma_{-})}=:S(\gamma_{+},\gamma_{-}),$

which is monotonically increasing with $\gamma_{+}$ and decreasing with
$\gamma_{-}$. This can be seen by computing the derivatives:

$\displaystyle\frac{\partial S}{\partial\gamma_{+}}$
$\displaystyle=\frac{\cos\gamma_{+}\sin(\gamma_{+}-\gamma_{-})-(\sin\gamma_{+}-\sin\gamma_{-})\cos(\gamma_{+}-\gamma_{-})}{\sin^{2}(\gamma_{+}-\gamma_{-})}$
$\displaystyle=\frac{\sin\gamma_{-}(\cos(\gamma_{+}-\gamma_{-})-1)}{\sin^{2}(\gamma_{+}-\gamma_{-})}\geq
0,$

and similarly $\frac{\partial S}{\partial\gamma_{-}}\leq 0$. Geometrically we
can easily see that $\gamma_{+}\geq\gamma_{1}$ and
$\gamma_{-}\leq\gamma_{N}-\pi$. Thus, the minimum of $|z^{+}|+|z^{-}|$ is
attained when $\gamma_{+}=\gamma_{1}$ and $\gamma_{-}=\gamma_{N}-\pi$. This
corresponds to $a_{2}=\dots=a_{N-1}=0$ and leads exactly to the RHS of Eq.
(A.1). ∎

With this lemma, we can explicitly write down the constraint on temperature
$\theta$. To do this, we denote by $\gamma_{\bm{U}}\in(-\pi,\pi]$ the argument
of the flow velocity $\bm{U}$ in Eq. (3.1). Then the second equality in Eq.
(3.1) can be rewritten as

$\sum_{m=1}^{N}\hat{\rho}_{m}\hat{u}_{m}e^{i\gamma_{m}}=\rho
Ue^{i\gamma_{\bm{U}}}$

or

$\sum_{m=1}^{N}s\hat{\rho}_{m}\hat{u}_{m}e^{i(\gamma_{m}-\eta)}=\rho U,$ (A.2)

where

$s=\left\\{\begin{aligned} 1,\quad&\text{if }0<\gamma_{\bm{U}}\leq\pi,\\\
-1,\quad&\text{if }-\pi<\gamma_{\bm{U}}\leq 0,\end{aligned}\right.$

and $\eta=\gamma_{\bm{U}}-\frac{s-1}{2}\pi\in(0,\pi]$.

Let $k$ ($1\leq k\leq N$) be such that

$\gamma_{k}<\eta\leq\gamma_{k+1}.$

For convenience, we define $\gamma_{N+1}=\pi+\gamma_{1}$. Then we set

$\tilde{\gamma}_{m}=\left\\{\begin{aligned} &\gamma_{m}+\pi-\eta,&1\leq m\leq
k,\\\ &\gamma_{m}-\eta,&k+1\leq m\leq N.\end{aligned}\right.$

It is clear that
$0\leq\tilde{\gamma}_{k+1}<\cdots<\tilde{\gamma}_{N}<\tilde{\gamma}_{1}<\cdots<\tilde{\gamma}_{k}<\pi$.
Now we apply the lemma above to Eq. (A.2) to get

$U^{2}\left(\frac{\sin\tilde{\gamma}_{k+1}+\sin\tilde{\gamma}_{k}}{\sin(\tilde{\gamma}_{k}-\tilde{\gamma}_{k+1})}\right)^{2}\leq\frac{1}{\rho^{2}}\left(\sum_{m=1}^{N}|\hat{\rho}_{m}\hat{u}_{m}|\right)^{2}<2E,$

where the last inequality is Eq. (3.2). The left-hand side is

$U^{2}\left(\frac{\sin(\gamma_{k+1}-\eta)+\sin(\eta-\gamma_{k})}{\sin(\gamma_{k+1}-\gamma_{k})}\right)^{2},$

which obviously goes to $U^{2}$ as $\gamma_{k+1}-\gamma_{k}\rightarrow 0$.
Moreover, the constraint for $\theta=(2E-U^{2})/2$ becomes

$\theta>U^{2}\frac{\sin(\gamma_{k+1}-\eta)\sin(\eta-\gamma_{k})}{1+\cos(\gamma_{k+1}-\gamma_{k})}.$
(A.3)

This lower bound approaches 0 when $\gamma_{k+1}-\gamma_{k}\to 0$.

## References

  * [1] S. Ansumali, I. V. Karlin and H. C. Öttinger, Minimal entropic kinetic models for hydrodynamics, Europhys. Lett. 63 (2003) 798–804.
  * [2] L. Armijo, Minimization of functions having Lipschitz continuous first partial derivatives, Pacific J. Math. 16 (1966) 1–3.
  * [3] N. Bellomo, D. Burini, G. Dosi, L. Gibelli, D. Knopoff, N. Outada, P. Terna and M.-E. Virgillito, What is life? A perspective of the mathematical kinetic theory of active particles, Math. Mod. Meth. Appl. Sci. 31 (2021) 1821–1866.
  * [4] N. Bellomo and R. Gatignol, From the Boltzmann equation to discretized kinetic models, in Lecture Notes on the Discretization of the Boltzmann Equation (World Scientific, 2003) pp. 1–16.
  * [5] P. Bhatnagar, E. Gross and 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.
  * [6] Z. Cai, Y. Fan and R. Li, A framework on moment model reduction for kinetic equation, SIAM J. Appl. Math., 75 (2015) 2001–2023.
  * [7] C. Chalons, R. Fox, F. Laurent, M. Massot and A. Vi$\acute{e}$, Multivariate Gaussian extended quadrature method of moments for turbulent disperse multiphase flow, Multiscale Model. Simul. 15 (2017) 1553–1583.
  * [8] Y. Chen, Q. Huang and W.-A. Yong, Discrete-velocity-direction models of BGK-type with minimum entropy: II. Flow simulations, in preparation.
  * [9] C.-K. Chu, Kinetic-theoretic description of the formation of a shock wave, Phys. Fluids 8 (1965) 12–22.
  * [10] H.B. Curry, The method of steepest descent for non-linear minimization problems, Quarterly of Applied Mathematics 2 (1944) 258–261.
  * [11] G. Dimarco, Q. Li, L. Pareschi and B. Yan, Numerical methods for plasma physics in collisional regimes, J. Plasma Phys., 81 (2015) 305810106.
  * [12] R. Fox, A quadrature-based third-order moment method for dilute gas-particle flows, J. Comput. Phys. 227 (2008) 6313–6350.
  * [13] S. Friedlander, Smoke, Dust, and Haze: Fundamentals of Aerosol Dynamics (2nd ed., Oxford University Press, 2000).
  * [14] R. Gatignol, Théorie Cinétique de Gaz à Répartition Discrète de Vitesses (Springer-Verlag, 1975).
  * [15] Z. Guo, K. Xu and R. Wang, Discrete unified gas kinetic scheme for all Knudsen number flows: Low-speed isothermal case, Phys. Rev. E 88 (2013) 033305.
  * [16] S. Harris, An Introduction to the Theory of the Boltzmann Equation (Dover Publications, 2004).
  * [17] L. Holway, New statistical models for kinetic theory: Methods of construction, Phys. Fluids, 9 (1966) 1658–1673.
  * [18] Q. Huang, S.Q. Li and W.-A. Yong, Stability analysis of quadrature-based moment methods for kinetic equations, SIAM J. Appl. Math. 80 (2020) 206–231.
  * [19] Q. Huang, P. Ma, Q. Gao and S.Q. Li, Ultrafine particle formation in pulverized coal, biomass, and waste combustion: Understanding the relationship with flame synthesis process, Energy Fuels 34 (2020) 1386–1395.
  * [20] W. Koller, Semi-continuous extended kinetic theory, In: Lecture Notes on the Discretization of the Boltzmann Equation (World Scientific, 2003) 97–132.
  * [21] C. Levermore, Moment closure hierarchies for kinetic theories, J. Statist. Phys., 83 (1996) 1021–1065.
  * [22] D. Li and D. Marchisio, Implementation of CHyQMOM in OpenFOAM for the simulation of non-equilibrium gas-particle flows under one-way and two-way coupling, Power Technol., 396 (2022) 765–784.
  * [23] D.Marchisio and R. Fox, Computational Models for Polydisperse Particulate and Multiphase Systems (Cambridge University Press, 2013).
  * [24] L. Mieussens, Discrete velocity model and implicit scheme for the BGK equation of rarefied gas dynamics, Math. Mod. Meth. Appl. Sci. 10 (2000) 1121–1149.
  * [25] L. Mieussens, Discrete-velocity models and numerical schemes for the Boltzmann-BGK equation in plane and axisymmetric geometries, J. Comput. Phys. 162 (2000) 429–466.
  * [26] L. Mieussens, Convergence of a discrete-velocity model for the Boltzmann-BGK equation, Computers and Mathematics with Applications 41 (2001) 83–96.
  * [27] S. Naris and D. Valougeorgis, The driven cavity flow over the whole range of the Knudsen number, Phys. Fluids 17 (2005) 097106.
  * [28] L. Preziosi and L. Rondoni, Discretization of the Boltzmann equation and the semicontinuous model, in Lecture Notes on the Discretization of the Boltzmann equation (World Scientific, 2003) pp. 59–95.
  * [29] E.M. Shakhov, Generalization of the Krook kinetic relaxation equation, Fluid Dyn., 3 (1968) 95–96.
  * [30] F. Sharipov, Rarefied Gas Dynamics: Fundamentals for Research and Practice (Wiley-VCH, 2015).
  * [31] F. Thüroff, C. Weber and E. Frey, Numerical treatment of the Boltzmann equation for self-propelled particle systems, Phys. Rev. X 4 (2014) 041030.
  * [32] E. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics (Springer-Verlag, 2009).
  * [33] W.-A. Yong, Singular perturbations of first-order hyperbolic systems with stiff source terms, J. Differential Equations 155 (1999) 89–132.
  * [34] C. Yuan, F. Laurent, R.O. Fox, An extended quadrature method of moments for population balance equations, J. Aerosol Sci. 51 (2012) 1–23.
  * [35] Z. Zhang, J. Xu, Z. Qi and G. Xi, A discrete velocity direction model for the Boltzmann equation and applications to micro gas flows, J. Comput. Phys. 227 (2008) 5256–5271.
  * [36] Z. Zhang, C. Peng and J. Xu, H theorem and sufficient conditions for the discrete velocity direction model, Mod. Phys. Lett. B 27 (2013) 1350007.
  * [37] Z. Zhang, W. Zhao, Q. Zhao, G. Lu and J. Xu, Inlet and outlet boundary conditions for the discrete velocity direction model, Mod. Phys. Lett. B 32 (2018) 1850048.