**2. Equilibrium for DV Boltzmann Kinetic Model and the Euler Equations**

The local equilibrium of the model (1) is obtained as a minimum of the *H* functional with the constraints corresponding to the conservation laws; it has the following form (Formula (5) in [43])

$$f\_i^{eq} = \exp\left(a + b \cdot c\_i + dc\_i^2\right), \quad i = 1 \ldots N\_\prime \tag{2}$$

where the coefficients *ai*, *bi*, *di* depend on the density, flow velocityand temperature *ρ*, *u*, *θ* and "·" defines scalar product. In this paragraph, we consider the particle's dynamics in *D* spatial dimensions. We assume that the local equilibrium is close to the absolute equilibrium with the density *ρ*<sup>0</sup> = 1 flow velocity *u*<sup>0</sup> = 0 and the temperature *θ*0; then, one can write down the absolute equilibrium denoted as *wi* in the form

$$w\_i = \exp(a^0 + d^0 c\_i^2), \quad i = 1 \dots N\_{\prime}$$

where *a*<sup>0</sup> = *a*0(*ρ*0, *θ*0), *d*<sup>0</sup> = *d*0(*ρ*0, *θ*0), we also term *wi* as lattice weights. The conservation laws for mass, momentum, energy yield the following equations for the lattice velocities and the lattice weights

$$
\sum\_{i} w\_{i} = 1, \quad \sum\_{i} w\_{i} \mathbf{c}\_{i} = 0, \quad \sum\_{i} w\_{i} \mathbf{c}\_{i} \mathbf{c}\_{i} = \theta\_{0} \delta\_{i} \tag{3}
$$

note that *cici* is a tensor with elements *ci*,*<sup>α</sup>ci*,*β*, *α*, *β* = 1 ... *D*, and that *D* is the number of spatial dimensions. For the coefficients, we have the expression

$$a\_i = a\_i^0 + \Delta a\_\prime \quad b\_i = \Delta b\_\prime \quad d\_i = d\_i^0 + \Delta d\_\prime$$

where Δ*a*, **Δ***b*, Δ*d* are small quantities. Similarly to the previous studies [43] we expand the expression (2) on Δ*a*, **Δ***b*, Δ*d*, and one has

$$f\_i^{eq} = w\_i \left(1 + \Delta a + \frac{1}{2} \Delta a^2 + o(\Delta a^2)\right) \times$$

$$\times \left(1 + \Delta b \cdot c\_i + \frac{1}{2} \Delta b \Delta b : c\_i c\_i + o(\Delta b^2)\right) \left(1 + \Delta d c\_i^2 + \frac{1}{2} \Delta d^2 c\_i^4 + o(\Delta d^2)\right) =$$

$$= w\_i \left(1 + \Delta a + \Delta b \cdot c\_i + \Delta d c\_i^2 +$$

$$+ \frac{1}{2} \Delta a^2 + \Delta a \Delta b \cdot c\_i + \Delta a \Delta d c\_i^2 + \frac{1}{2} \Delta b \Delta b : c\_i c\_i + \Delta d \Delta b \cdot c\_i c\_i^2 + \frac{1}{2} \Delta d^2 c\_i^4\right) + o(\Delta^2) \tag{4}$$

where *o*(Δ2) stands for *o*(Δ*a*2), *o*(**Δ***b*2), *o*(Δ*d*2), also the operator ":" is tensor convolution. Next, we assume that

$$
\rho = 1 + \Delta \rho\_{\prime} \quad \mu = \Delta \mu\_{\prime} \quad \theta = \theta\_0 + \Delta \theta\_{\prime},
$$

where Δ*ρ*, Δ*u*, Δ*θ* are small. Taking (4) into account, for the first local equilibrium moments, we derive the following equations

$$\sum\_{i} f\_i^{eq} = 1 + \Delta \rho = 1 + \Delta a + D\theta\_0 \Delta d + \frac{1}{2} \Delta a^2 + D\theta\_0 \Delta a \Delta d + \frac{\theta\_0}{2} \Delta b^2 + \frac{m\_4}{2} \Delta d^2,\tag{5}$$

$$\sum\_{i} f\_i^{eq} c\_i = \Delta \mu + \Delta \rho \Delta \mu = \theta\_0 \Delta b + \theta\_0 \Delta a \Delta b + D^{-1} m\_4 \Delta d \Delta b,\tag{6}$$

$$\sum\_{i} f\_i^{eq} c\_i^2 = D\theta\_0 + D\theta\_0 \Delta \rho + D\Delta \theta + \Delta u^2 + D\Delta \rho \Delta \theta = $$

$$= D\theta\_0 \left( 1 + \Delta a + \frac{1}{2} \Delta a^2 \right) + m\_4 \left( \Delta d + \Delta a \Delta d + \frac{1}{2D} \Delta b^2 \right) + \frac{m\_6}{2} \Delta d^2,\tag{7}$$

and we omitted third-order terms and used the definitions

$$m\_4 = \sum\_i w\_i c\_{i'}^4, \quad m\_6 = \sum\_i w\_i c\_{i'}^6. \tag{8}$$

We assume that Δ*ρ*, Δ*u*, Δ*θ* are of the same order of smallness, which we define as *O*(Δ). Then, we seek solutions to the Equations (5)–(7) in the form

$$
\Delta a = \Delta a\_{\rm lin} + \Delta a\_{\rm noml}, \quad \Delta b = \Delta b\_{\rm lin} + \Delta b\_{\rm nonl}, \quad \Delta d = \Delta d\_{\rm lin} + \Delta d\_{\rm nonl}.
$$

where the terms Δ*alin*, **Δ***blin*, Δ*dlin* are solutions to the linearized Equations (5)–(7) of order *O*(Δ) and Δ*anonl*, **Δ***bnonl*, Δ*dnonl* are nonlinear corrections of order *O*(Δ2). First, from the linearized Equations (5)–(7), one has

$$
\Delta \rho = \Delta a\_{\rm lin} + D \theta\_0 \Delta d\_{\rm lin}, \quad \Delta u = \theta\_0 \Delta b\_{\rm lin}, \quad D \theta\_0 \Delta \rho + D \Delta \theta = D \theta\_0 \Delta a\_{\rm lin} + m\_4 \Delta d\_{\rm lim},
$$

these equations have the solutions

$$
\Delta a\_{lin} = \Delta \rho - \frac{D^2 \theta\_0 \Delta \theta}{m\_4 - D^2 \theta\_0^2}, \quad \Delta b\_{lin} = \frac{\Delta \mathbf{u}}{\theta\_0}, \quad \Delta d\_{lin} = \frac{D \Delta \theta}{m\_4 - D^2 \theta\_0^2}. \tag{9}
$$

Next, we find the nonlinear corrections Δ*anonl*, **Δ***bnonl*, Δ*dnonl* from the Equations (5)–(7). It would be convenient to start with the Equation (6), which can be rewritten as

$$
\Delta\rho\Delta u = \theta\_0 \Delta a\_{lin} \Delta b\_{lin} + D^{-1} m\_4 \Delta d\_{lin} \Delta b\_{lin} + \theta\_0 \Delta b\_{nonl,\*}
$$

from the last equation, we immediately obtain

$$
\Delta b\_{nom} = -\frac{\Delta \theta \Delta \mathfrak{u}}{\theta\_0^2}.\tag{10}
$$

Consideration of the Equations (5) and (7) yields

$$
\Delta a\_{nonl} + D\theta\_0 \Delta d\_{nonl} + \frac{1}{2} \Delta a\_{lin}^2 + D\theta\_0 \Delta a\_{lin} \Delta d\_{lin} + \frac{1}{2} \theta\_0 \Delta b\_{lin}^2 + \frac{m\_4}{2} \Delta d\_{lin}^2 = 0,
$$

*Dθ*0Δ*anonl* + *m*4Δ*dnonl* + *Dθ*<sup>0</sup> <sup>2</sup> <sup>Δ</sup>*a*<sup>2</sup> *lin* <sup>+</sup> *<sup>m</sup>*4Δ*alin*Δ*dlin* <sup>+</sup> *<sup>m</sup>*<sup>4</sup> 2*D***Δ***b*<sup>2</sup> *lin* <sup>+</sup> *<sup>m</sup>*<sup>6</sup> <sup>2</sup> <sup>Δ</sup>*d*<sup>2</sup> *lin* <sup>=</sup> <sup>Δ</sup>*u*<sup>2</sup> <sup>+</sup> *<sup>D</sup>*Δ*ρ*Δ*θ*,

by applying (9) we get the solutions

$$\Delta a\_{nvol} = -\frac{\Delta \rho^2}{2} - \frac{D\theta\_0 \Delta \mathbf{u}^2}{(m\_4 - D^2 \theta\_0^2)} - \frac{D^4 \theta\_0^2 \Delta \theta^2}{2(m\_4 - D^2 \theta\_0^2)^2} + \frac{(D\theta\_0 m\_6 - m\_4^2)D^2 \Delta \theta^2}{2(m\_4 - D^2 \theta\_0^2)^3},\tag{11}$$

$$
\Delta d\_{nonl} = -\frac{\Delta u^2}{2D\theta\_0^2} + \frac{\Delta u^2}{(m\_4 - D^2\theta\_0^2)} + \frac{D^3 \theta\_0 \Delta \theta^2}{(m\_4 - D^2\theta\_0^2)^2} - \frac{(m\_6 - D\theta\_0 m\_4)D^2 \Delta \theta^2}{2(m\_4 - D^2\theta\_0^2)^3}.\tag{12}
$$

The combination of (9) and (10)–(12) leads to the following expression for *f eq i*

**Proposition 1.** *The DV local equilibrium f eq <sup>i</sup> in the form (2) can be expressed as*

$$f\_i^{eq} = w\_i(k\_0 + k\_1 \cdot \mathbf{c}\_i + k\_2 \mathbf{c}\_i^2 + k\_3 : \mathbf{c}\_i \mathbf{c}\_i + k\_4 \cdot \mathbf{c}\_i \mathbf{c}\_i^2 + k\_5 \mathbf{c}\_i^4) + O(\Delta^3), \quad i = 1 \ldots N,\tag{13}$$

*where*

$$k\_0 = 1 + \Delta \rho - \frac{D\theta\_0}{(m\_4 - D^2 \theta\_0^2)} (D\Delta \theta + D\Delta \rho \Delta \theta + \Delta \mu^2) + \frac{D\theta\_0 m\_6 - m\_4^2}{2(m\_4 - D^2 \theta\_0^2)^3} D^2 \Delta \theta^2,\tag{14}$$

$$k\_1 = \frac{\Delta \mu}{\theta\_0} + \frac{\Delta \rho \Delta \mu}{\theta\_0} - \frac{m\_4}{\theta\_0^2 (m\_4 - D^2 \theta\_0^2)} \Delta \theta \Delta \mu,\tag{15}$$

$$k\_{2} = -\frac{\Delta\mathfrak{u}^{2}}{2D\theta\_{0}^{2}} + \frac{1}{(m\_{4} - D^{2}\theta\_{0}^{2})} (D\Delta\theta + D\Delta\rho\Delta\theta + \Delta\mathfrak{u}^{2}) - \frac{m\_{6} - D\theta\_{0}m\_{4}}{2(m\_{4} - D^{2}\theta\_{0}^{2})^{3}} D^{2}\Delta\theta^{2},\tag{16}$$

$$k\_3 = \frac{\Delta\mu\Delta u}{2\theta\_0^2},\tag{17}$$

$$\mathbf{k} \mathbf{k} = \frac{D \Delta \theta \Delta \mathbf{u}}{\theta\_0 (m\_4 - D^2 \theta\_0^2)},\tag{18}$$

$$k\_5 = \frac{D^2 \Delta \theta^2}{2(m\_4 - D^2 \theta\_0^2)^2} \tag{19}$$

*and* Δ*ρ*, Δ*u* Δ*θ are small density, flow velocity and temperature variations of order O*(Δ)*; moreover, the moments m*<sup>4</sup> *and m*<sup>6</sup> *are defined by (8), and the absolute equilibrium wi satisfies the conditions*

$$
\sum\_{i} w\_i = 1, \quad \sum\_{i} w\_i \mathbf{c}\_i = 0, \quad \sum\_{i} w\_i \mathbf{c}\_i \mathbf{c}\_i = \theta\_0 \delta\_i
$$

*where θ*<sup>0</sup> *is the reference temperature.*

By applying (13)–(19), the pressure tensor *P* with the components *Pαβ*, *α*, *β* = 1 ... *D* at the local equilibrium can be derived

$$P = \sum\_{i} f\_i^{eq} (\mathbf{c}\_i - \Delta \mathbf{u})(\mathbf{c}\_i - \Delta \mathbf{u}) = \sum\_{i} f\_i^{eq} (\mathbf{c}\_i \mathbf{c}\_i - 2\mathbf{c}\_i \Delta \mathbf{u} + \Delta \mathbf{u} \Delta \mathbf{u}) = 0$$

$$= (\theta\_0 k\_0 + D^{-1} m\_4 k\_2 + D^{-1} m\_6 k\_5) \delta + \mathbf{k}\_3 : \mathbf{R} - \Delta \mathbf{u} \Delta \mathbf{u} + O(\Delta^3) = 0$$

$$= \rho \theta \delta + \frac{2D \theta\_0^2 - m\_4}{2D^2 \theta\_0^2} \Delta \mathbf{u}^2 \delta + \frac{\Delta \mathbf{u} \Delta \mathbf{u}}{2 \theta\_0^2} : \mathbf{R} - \Delta \mathbf{u} \Delta \mathbf{u} + O(\Delta^3),$$

where *ρθ* = (1 + Δ*ρ*)(*θ*<sup>0</sup> + Δ*θ*) and

$$\mathcal{R} = \sum\_{i} w\_i c\_i c\_i c\_i c\_i c\_i.$$

Now, let us assume that *R* is isotropic tensor; in such a case, its components can be written in the form (Formula (69) in [37])

$$R\_{a\beta\lambda\gamma} = \frac{m\_4}{D(D+2)} (\delta\_{a\beta}\delta\_{\lambda\gamma} + \delta\_{a\gamma}\delta\_{\beta\lambda} + \delta\_{a\lambda}\delta\_{\beta\gamma}),\tag{20}$$

one can see that the tensor *P* equalling *ρθδ* + *O*(Δ3) is obtained, if

$$m\_4 = D(D+2)\theta\_0^2. \tag{21}$$

Compared to *P* for the local Maxwell distribution, in the DV approach, the error *O*(Δ3) is observed; therefore, we can conclude that the hydrodynamics (mass and momentum equations) at the Euler level of accuracy is recovered with the errors of order *O*(Δ3) if the conditions (3), (20) and (21) are satisfied.

Finally, we consider the heat flow *q* at the level of the Euler equations

$$2\boldsymbol{\eta} = \sum\_{i} f\_i^{eq} (\mathbf{c}\_i - \Delta \boldsymbol{u})^2 (\mathbf{c}\_i - \Delta \boldsymbol{u}) = \sum\_{i} f\_i^{eq} c\_i^2 \mathbf{c}\_i - 2\Delta \boldsymbol{u} \cdot \mathbf{P} - 2\rho E \Delta \boldsymbol{u} + O(\Delta^3),$$

where *E* = (*ρ*/2)(*Dθ* + Δ*u*2), *θ* = *θ*<sup>0</sup> + Δ*θ*, *ρ* = 1 + Δ*ρ*, applying (13)–(19) we obtain

$$2\mathfrak{q} = \frac{1}{D\theta\_0} \left( m\_4 + m\_4 \Delta \rho + \frac{D\theta\_0 m\_6 - m\_4^2}{\theta\_0 (m\_4 - D^2 \theta\_0^2)} \Delta \theta \right) \Delta \mathfrak{u} - 1$$

$$-2(\theta\_0 + \Delta\theta + \theta\_0 \Delta\rho)\Delta\mathfrak{u}\cdot\mathcal{S} - (D\theta\_0 + D\theta\_0 \Delta\rho + D\Delta\theta)\Delta\mathfrak{u} + O(\Delta^3) = $$

$$= \frac{1}{D\theta\_0} (m\_4 - D(D+2)\theta\_0^2)(\Delta\mathfrak{u} + \Delta\rho\Delta\mathfrak{u}) + \left(\frac{D\theta\_0 m\_6 - m\_4^2}{D\theta\_0^2 (m\_4 - D^2\theta\_0^2)} - (D+2)\right)\Delta\theta\Delta\mathfrak{u} + O(\Delta^3),$$

one can see that the terms proportional Δ*u*, Δ*ρ*Δ*u* are eliminated if *m*<sup>4</sup> satisfies (21), in addition, the second term can be removed if *m*<sup>6</sup> = *D*(*D* + 2)(*D* + 4)*θ*<sup>3</sup> <sup>0</sup> or we are restricted by the isothermal flows Δ*θ* = 0; in such a case, the heat flow (which equals zero for the Euler equations) is only of order *O*(Δ3). In the present study, we assume that the temperature variations are negligible, Δ*θ* = 0.

#### **3. Navier–Stokes Equations**

In order to obtain the Navier–Stokes equations, one needs to find the corrections to the pressure tensor corresponding to the viscous terms. This can be performed by applying the Chapman–Enskog expansion for DV Boltzmann model [29]. Then, following the previous results [29], we assume that the solution to (1) can be expressed in the form

*fi* = *f eq <sup>i</sup>* + *f* (1) *<sup>i</sup>* + *<sup>O</sup>*(*Kn*2), where *<sup>f</sup>* (1) *<sup>i</sup>* are of order *O*(*Kn*) and *Kn* is the Knudsen number. At the limit of small Mach numbers, the equations for *f* (1) *<sup>i</sup>* read as (Equation (19) in [43])

$$\frac{df\_i^{eq}}{dt} = \sum\_{jkl}^{N} A\_{kl}^{ij} (w\_k f\_l^{(1)} + w\_l f\_k^{(1)} - w\_i f\_j^{(1)} - w\_j f\_i^{(1)}), \quad i = 1 \ldots N,\tag{22}$$

one can see from (22) that the solutions *f* (1) *<sup>i</sup>* are determined by the concrete DV Boltzmann model, i.e., *f* (1) *<sup>i</sup>* depend on *<sup>A</sup>ij kl*. The solution to the linear Equations (22) can be obtained as (Formula (22) in [43])

$$f\_i^{(1)} = a\_i \mathbf{Q}\_i : \frac{\partial}{\partial r} \Delta u + b\_i \text{div}(\Delta u), \tag{23}$$

where *Q<sup>i</sup>* is a second-order tensor whose exact form we will discuss further, and *ai*, *bi* are numerical coefficients.

#### **4. Spurious Invariants**

For a collision, in which the particles with the velocities *ci*, *cj* turn into the particles with the velocities *ck*, *cl*, we introduce the following **reaction vector** [40,41,44]

$$\mathfrak{e} = (\ldots, \widehat{1^{\lambda}}, \ldots, \widehat{-1^{\lambda}}, \ldots, \widehat{-1^{\lambda}}, \ldots, \widehat{-1^{\lambda}}, \ldots) \in \mathbb{R}^{N}, \quad A\_{kl}^{ij} > 0,$$

where the entries denoted by dots equal zero. Assume that we have *p* linearly independent reaction vectors *es*,*s* = 1 ... *p* . We denote a matrix consisting of all reaction row vectors *e<sup>s</sup>* as the **collision matrix**

$$\mathbb{C} = (e\_1; e\_2; \dots; e\_p) \in \mathbb{R}^p \times \mathbb{R}^N.$$

Note that the collision invariants *<sup>ϕ</sup>*(*c***1**,... *cN*) <sup>∈</sup> *<sup>R</sup><sup>N</sup>* are defined by the relation [28]

$$
\varphi\_i + \varphi\_j = \varphi\_k + \varphi\_{l\prime} \quad A\_{kl}^{ij} > 0,
$$

this condition can also be rewritten in the following form [44,45]

$$
\mathfrak{p} \cdot \mathfrak{e}\_{\mathfrak{s}} = 0, \quad \mathfrak{s} = 1 \dots p,\tag{24}
$$

i.e., the linear subspace spanned by the invariants is orthogonal to the subspace spanned by the reaction vectors. The condition (24) can be applied for the detection of spurious invariants:

**Proposition 2.** *Assume that, for some DV Boltzmann models, the number of linearly independent physical collision invariants equals q, then additional invariants do not exist if [40,41,44,45]*

$$\operatorname{rank}(\mathbf{C}) = N - q\_{\prime\prime}$$

*where N is the number of the discrete velocities.*
