††thanks: Corresponding author

# Entropic equilibrium for the lattice Boltzmann method:
Hydrodynamics and numerical properties

S. A. Hosseini  Department of Mechanical and Process Engineering, ETH Zurich,
8092 Zurich, Switzerland  I. V. Karlin<EMAIL_ADDRESS>Department of
Mechanical and Process Engineering, ETH Zurich, 8092 Zurich, Switzerland

###### Abstract

The entropic lattice Boltzmann framework proposed the construction of the
discrete equilibrium by taking into consideration minimization of a discrete
entropy functional. The effect of this form of the discrete equilibrium on
properties of the resulting solver has been the topic of discussions in the
literature. Here we present a rigorous analysis of the hydrodynamics and
numerics of the entropic. In doing so we demonstrate that the entropic
equilibrium features unconditional linear stability, in contrast to the
conventional polynomial equilibrium. We reveal the mechanisms through which
unconditional linear stability is guaranteed, most notable of which the
adaptive normal modes propagation velocity and the positive-definite nature of
the dissipation rates of all eigen-modes. We further present a simple local
correction to considerably reduce the deviations in the effective bulk
viscosity.

###### pacs:

47.11.-j

##  I Introduction

The lattice Boltzmann method is a numerical method developed in the late
80’s/early 90’s, as an alternative to classical solvers for the Navier-Stokes
equations, initially, in the incompressible flow limit [1]. This approach
finds its roots in the kinetic theory of gases and is essentially a solver for
the phase-space discretized Boltzmann equation with the linear
Bhatnagar–Gross–Krook (BGK) approximation for the collision term [2]. The
relatively uncomplicated algorithm, low cost of discrete operations in the
time/space evolution equations, their locality stemming from the purely
hyperbolic nature of the corresponding equations -as compared to elliptic-
hyperbolic alternatives such as the Poisson–Navier–Stokes equations, and
properties such as low numerical dissipation and strictly conservative nature
have pushed the lattice Boltzmann method to the forefront of computational
fluid dynamics.
Early on after its first appearance in the literature, numerical stability
became a major point of discussion and research. In its original form, i.e.
the single relaxation time BGK model with a quadratic-in-velocity discrete
equilibrium, the lattice Boltzmann solver has been reported to be sensitive to
the viscosity non-dimensionalized by the time-step and grid sizes, also
referred to as the Fourier number, and the maximum non-dimensional velocity
[3]. Note that such limitations are not proper to the lattice Boltzmann method
and come with any other numerical approximation to a system of
hyperbolic/parabolic partial differential equations, as illustrated by the Lax
equivalence theorem for finite difference methods [4, 5] and the Lax-Richtmyer
stability condition [6]. Plethora of modifications to the original lattice BGK
(LBGK) have been proposed since, the majority of which can be categorized as
some form of multiple relaxation times putting in effect the burden of
instabilities onto relaxation rates of individual independent moments of the
distribution function [7, 8, 9, 10]. While such schemes have been successful,
to different and limited extent, in extending the operation range of the LBGK,
none of them has actually resulted in a scheme that is unconditionally stable
in the linear regime. Based on that observation one might be tempted to ask
whether focus should not be put on other components of the LBGK algorithm, for
instance the discrete equilibrium. While the continuous equilibrium
distribution function, i.e. Maxwell-Boltzmann distribution, is a minimizer of
entropy and the Boltzmann-BGK equation complies with the H-theorem, it is not
guaranteed that truncated expansions of the distribution function in a given
projection space will abide by these conditions, i.e. non-commutativity of
entropy minimization and discretization in phase-space. The entropic lattice
Boltzmann method [11] proposed a lattice Boltzmann realization guaranteeing
stability based on a construction of the equilibrium attractor taking into
account minimization of the discrete entropy. While the entropic lattice
Boltzmann method has been successfully used for a wide variety of large
Reynolds number simulations, see for instance [12, 13, 14], the question of
its stability as compared to other lattice Boltzmann realizations has been a
topic of debates in the literature, see for instance [15, 16, 17, 18], as best
illustrated by the following quote from [16]: ”… it is unclear theoretically
how the ELBE with a constant relaxation parameter $\tau$ can improve the
numerical stability of the LBGK scheme.”
This sentence summarizes, in a nutshell, the question that will be considered
in the present manuscript. Does the form of the discrete equilibrium
distribution function, for a fixed relaxation scheme, affect the linear
stability of the solver? Can the specific form of the discrete equilibrium
derived by minimization of discrete entropy guarantee unconditional linear
stability? And if so what are the mechanisms by which the entropic equilibrium
guarantees linear stability? The present work will address these questions
with a detailed analysis of different forms of the discrete equilibrium
distribution function. Note that the present work only addresses the question
of linear stability, as non-linear stability in the entropic lattice Boltzmann
method is brought about by another ingredient which is dynamic restriction of
the maximum path-length in the relaxation process.
The manuscript is organized as follows: after a brief review of basic concepts
of the lattice Boltzmann method and different discrete equilibria in section
II, the continuum-level dynamics of the equivalent macroscopic balance
equations will be analyzed. First dispersion properties will be dissected by
looking into the Euler-level equations in section III. Dissipation at the
Navier-Stokes level will then be discussed in section IV. In section V the
analyses will be extended to the entire wave-number space through the von
Neumann method and continuum limit results of sections IV and V will be
confirmed via numerical simulations. In addition, a simple and local
correction is proposed to considerably reduce the deviations in the effective
bulk viscosity for all equilibria. The manuscript will end with section VI
summarizing all observation and discussing larger impact of the present study.

##  II Basic concepts

In the remainder of the present article, unless stated otherwise, we will
consider single relaxation time lattice Boltzmann models with the following
discrete time-evolution equation:

$f_{i}(\bm{r}+\bm{c}_{i}\delta t,t+\delta
t)-f_{i}(\bm{r},t)=2\beta\left(f_{i}^{\rm
eq}(\rho,\bm{u})-f_{i}(\bm{r},t)\right).$ (1)

Here $f_{i}$ are the discrete distribution functions, $\rm{r}$ the position in
space, $t$ time, $\delta t$ the time-step size, $\rho$ the fluid density,
$\bm{u}$ the velocity and $\beta$ the relaxation frequency defined as [19,
20]:

$\beta=\frac{\delta t}{2\nu/\varsigma^{2}+\delta t},$ (2)

where $\nu$ is the fluid kinematic viscosity and $\varsigma$ the lattice sound
speed:

$\varsigma=\frac{\delta r}{\sqrt{3}\delta t},$ (3)

and $\delta r$ the grid size. $f_{i}^{\rm eq}$ are the equilibrium
distribution functions which will be the focus of this work. Here all analyses
and discussions will consider first-neighbor lattices, in 1-D,

$c_{i}\in\frac{\delta r}{\delta t}\\{-1,0,1\\},$ (4)

and corresponding tensorial products in 2- and 3-D. The weights associated
with this lattice are, in 1-D,

$w_{i}\in\left\\{\frac{1}{6},\frac{2}{3},\frac{1}{6}\right\\}.$ (5)

Weights for 2- and 3-D lattices can be obtained as products of 1-D weights.
Details about the equilibrium distribution functions, starting with the
entropic construction are given in the next subsections.

###  II.1 Discrete entropic equilibrium construction

In the context of the entropic lattice Boltzmann method the discrete
equilibrium state is found as the minimizer of a Lyapunov functional $H$ [21,
11]:

$H=\sum_{i=1}^{Q}h_{i}(f_{i}),$ (6)

where $h_{i}$ are convex functions, under mass and momentum conservation
constraints:

$\displaystyle\sum_{i=1}^{Q}f_{i}^{\rm eq}=\rho,$ (7)
$\displaystyle\sum_{i=1}^{Q}\bm{c}_{i}f_{i}^{\rm eq}=\rho\bm{u},$ (8)

which is different from the polynomial construction which also imposes a
constraint on the second-order moment of the equilibrium distribution
function. Introducing the corresponding Lagrange multipliers,
$\lambda_{0},\lambda_{\alpha}$:

$\delta\sum_{i=1}^{Q}\left(h_{i}(f_{i})-\lambda_{0}f_{i}-\lambda_{\alpha}c_{i\alpha}f_{i}\right)=0,$
(9)

which yields

$h_{i}^{\prime}(f_{i}^{\rm eq})=\lambda_{0}+\lambda_{\alpha}c_{i\alpha}.$ (10)

Defining the inverse of $h_{i}^{\prime}(f_{i}^{\rm eq})$ as
$\mu_{i}=\left[{h_{i}^{\prime}(f_{i}^{\rm eq})}\right]^{-1}$, which must exist
due to the convexity of $h_{i}$ the formal solution of the minimization
problem reads:

$f_{i}^{\rm eq}=\mu_{i}\left(\lambda_{0}+\lambda_{\alpha}c_{i\alpha}\right).$
(11)

Taking the $H$-function to be [22]:

$H=\sum_{i=1}^{Q}f_{i}\ln(f_{i}/w_{i}),$ (12)

and considering stencils with $3^{D}$ discrete velocities the corresponding
equilibria can be expressed as:

$f_{i}^{\rm
eq}=w_{i}\exp\left(\lambda_{0}\right)\prod_{\alpha=1}^{D}\exp\left(c_{i\alpha}\lambda_{\alpha}\right).$
(13)

Introducing the following changes of variables,
$X=\exp\left(-\lambda_{0}\right)$ and
$Z_{\alpha}=\exp\left(\lambda_{\alpha}\right)$ the equilibrium can be re-
written as:

$f_{i}^{\rm eq}=w_{i}X^{-1}\prod_{\alpha=1}^{D}Z_{\alpha}^{c_{i\alpha}}.$ (14)

Here $X$ and $Z_{\alpha}$, i.e. the corresponding Lagrange multipliers, can be
obtained by writing the constraints on the moments:

$\displaystyle\rho X$
$\displaystyle=\sum_{i=1}^{Q}w_{i}\prod_{\alpha=1}^{D}Z_{\alpha}^{c_{i\alpha}},$
(15a) $\displaystyle\rho u_{\beta}X$
$\displaystyle=\sum_{i=1}^{Q}w_{i}c_{i\beta}\prod_{\alpha=1}^{D}Z_{\alpha}^{c_{i\alpha}}.$
(15b)

Solving this system of equations results in:

$Z_{\alpha}=\frac{2u_{\alpha}+\sqrt{{\left(\frac{u_{\alpha}}{\varsigma}\right)}^{2}+1}}{1-u_{\alpha}},$
(16)
$X^{-1}=\rho\prod_{\alpha=1}^{D}\left(2-\sqrt{{\left(\frac{u_{\alpha}}{\varsigma}\right)}^{2}+1}\right),$
(17)

leading to the following final expression for the discrete entropic
equilibrium [22]:

$f_{i}^{\rm
eq}=w_{i}\rho\prod_{\alpha=1}^{D}\left(2-\sqrt{{\left(\frac{u_{\alpha}}{\varsigma}\right)}^{2}+1}\right){\left(\frac{2u_{\alpha}+\sqrt{{\left(\frac{u_{\alpha}}{\varsigma}\right)}^{2}+1}}{1-u_{\alpha}}\right)}^{c_{i\alpha}}.$
(18)

###  II.2 Polynomial discrete equilibria

Another class of discrete equilibria, widely used in discrete velocity
approaches, are the polynomial equilibria. The most well-known realization of
polynomial equilibria for isothermal flows is the second-order polynomial
equilibrium defined as [3]:

$f_{i}^{\rm
eq}=w_{i}\rho\left(1+\frac{\bm{c}_{i}\cdot\bm{u}}{\varsigma^{2}}+\frac{{\left(\bm{c}_{i}\cdot\bm{u}\right)}^{2}}{2\varsigma^{4}}-\frac{\bm{u}^{2}}{2\varsigma^{2}}\right).$
(19)

While initially derived via a second-order Taylor expansion of the Maxwell-
Boltzmann distribution around Ma=0, it can also be obtained as a second-order
Hermite expansion of the Maxwell-Boltzmann distribution function.
This form of the equilibrium is known to lead to Galilean-variant errors in
the off-diagonal components of the viscous stress tensor, which scale with
order three in Mach number when the local temperature coincides with the
lattice reference temperature and order one when the local temperature is
different from the lattice reference. Attempts at improving upon the second-
order polynomial equilibrium, especially for compressible flows, have led to
the product form equilibrium:

$f_{i}^{\rm
eq}=\rho\prod_{\alpha=1}^{D}{\left(1-\varsigma^{2}-u_{\alpha}^{2}\right)}^{1-\lvert
c_{i\alpha}\lvert}{\left(\frac{c_{i\alpha}u_{\alpha}+\varsigma^{2}+u_{\alpha}^{2}}{2}\right)}^{\lvert
c_{i\alpha}\lvert}.$ (20)

This form of the equilibrium was obtained in [23] by _guiding_ the
minimization problem of the previous section by explicitly adding a constraint
on diagonal components of the second-order equilibrium moment tensor. This
form of the equilibrium can also be obtained, in 2-D, via a fourth order
Hermite expansion [24], or the moment matching method [25].

##  III Euler-level dynamics: dispersion

In this section we will discuss Euler level dynamics of the hydrodynamic
equations corresponding to each discrete equilibrium distribution function.

###  III.1 Conservation equations

Applying the Chapmann-Enkog analysis to the entropic lattice Boltzmann
equations, at order $\varepsilon$ one finds the following conservation laws:

$\displaystyle\partial_{t}^{(1)}\rho+\bm{\nabla}\cdot\rho\bm{u}$
$\displaystyle=0,$ (21a)
$\displaystyle\partial_{t}^{(1)}\rho\bm{u}+\bm{\nabla}\cdot\rho\bm{u}\otimes\bm{u}+\bm{\nabla}\cdot\rho\varsigma^{2}\bm{I}+\bm{\nabla}\cdot\bm{P}^{*}$
$\displaystyle=0.$ (21b)

Here $\bm{P}^{*}$ is a diagonal matrix defined as:

$P^{*}_{\alpha\alpha}=\widetilde{\Pi}_{\alpha\alpha}^{\rm
eq}-\rho\varsigma^{2}=\sum_{i=1}^{Q}{(\bm{c}_{i}-\bm{u})}^{2}f^{\rm
eq}_{i}-\rho\varsigma^{2},$ (22)

Using Eq. (18), $P^{*}_{\alpha\alpha}$ is readily computed:

$P^{*}_{\alpha\alpha}=-\rho\varsigma^{2}\left[{(u_{\alpha}/\varsigma)}^{2}-2\sqrt{{(u_{\alpha}/\varsigma)}^{2}+1}+2\right],$
(23)

while for Eqs. (19) and (20) $P^{*}_{\alpha\alpha}=0$. As defined here,
$P^{*}_{\alpha\alpha}$ is the deviation of the pressure from the
Maxwell–Boltzmann pressure, which is both flow-dependent and anisotropic for
the entropic equilibrium. Note that in the limit of
$u_{\alpha}\rightarrow\pm\delta r/\delta t$ the deviation,
$P^{*}_{\alpha\alpha}$ goes to $-\rho\varsigma^{2}$ resulting in a vanishing
total thermodynamic pressure. The trace of the deviation is illustrated in
Fig. 1 via the following normalized variable:

$\delta=\frac{P^{*}_{xx}+P^{*}_{yy}}{2\rho\varsigma^{2}}.$ (24)

Figure 1: Illustration of deviations in trace of pressure tensor for the
entropic equilibrium. Left: Normalized deviation in 2-D. Right: Normalized
error in 1-D.

Based on Fig. 1, as is well known, the deviation of pressure implied by the
entropic equilibrium are immaterial for nearly-incompressible flow
simulations. However, as we shall see in the analysis below, precisely these
deviations have crucial implications on the well-posedness and stability of
LBGK with the entropic equilibrium.

###  III.2 Characteristics analysis: Eigen-modes

In the case of polynomial equilibria of order equal to or larger than two the
second-order central moment, representing pressure, is only a function of
local density and temperature:

$\widetilde{\Pi}_{\alpha\alpha}=\rho\varsigma^{2},$ (25)

and the corresponding sound speed, $c_{s}$, can be readily computed via:

$c_{s}=u_{\alpha}\pm\sqrt{\frac{\partial\widetilde{\Pi}_{\alpha\alpha}}{\partial\rho}\bigg{|}_{T={\rm
cst}}}=u_{\alpha}\pm\varsigma.$ (26)

In the case of the entropic equilibrium the equilibrium pressure is also
function of local velocity making computation of the sound speed less evident.
To compute an analytical expression for the sound speed with the entropic
equilibrium, we use the method of characteristics, briefly introduced in the
next section.

####  III.2.1 Characteristics analysis formalism

Consider a system of conservation equations consisting of first-order partial
differential equations of the following form [26]:

$\partial_{t}\Phi+\partial_{x}\mathcal{F}+\mathcal{C}=0,$ (27)

where $\Phi$ is the vector of conserved variables, $\mathcal{F}$ the vector of
fluxes and $\mathcal{C}$ represents non-homogeneous terms in the equations.
This system, which is represented above with the conservative form, can also
be represented using the primitive form as:

$\partial_{t}\phi+\mathcal{A}\partial_{x}\phi+\mathcal{C}^{\prime}=0,$ (28)

where $\phi$ is the vector of primitive variable and $\mathcal{A}$ is a matrix
of size $n\times n$ -for a system with $n$ independent variables. Note that
unlike conserved variables, the choice of primitive variable is not unique.
The two descriptions are related to eachother via:

$\displaystyle\partial_{t}\Phi$ $\displaystyle=\mathcal{P}\partial_{t}\phi,$
(29a) $\displaystyle\partial_{x}\mathcal{F}$
$\displaystyle=\mathcal{Q}\partial_{x}\phi,$ (29b)

with:

$\displaystyle\mathcal{P}_{mn}$
$\displaystyle=\frac{\partial\Phi_{m}}{\partial\phi_{n}}$ (30a)
$\displaystyle\mathcal{Q}_{mn}$
$\displaystyle=\frac{\partial\mathcal{F}_{m}}{\partial\phi_{n}},$ (30b)

and:

$\displaystyle\mathcal{A}$ $\displaystyle=\mathcal{P}^{-1}\mathcal{Q},$ (31a)
$\displaystyle\mathcal{C}$
$\displaystyle=\mathcal{P}^{-1}\mathcal{C}^{\prime}.$ (31b)

Let us now introduce the left, $\bm{l}$, and right, $\bm{r}$, eigen-vectors of
$\mathcal{A}$ such that:

$\displaystyle l_{m}\mathcal{A}$ $\displaystyle=\lambda_{m}l_{m},$ (32a)
$\displaystyle\mathcal{A}r_{m}$ $\displaystyle=\lambda_{m}r_{m}.$ (32b)

where $\lambda_{i}$ are the eigen-values of $\mathcal{A}$. Then a diagonal
matrix $\Lambda$ can be obtained via:

$\Lambda=\mathcal{S}\mathcal{A}\mathcal{S}^{-1},$ (33)

where the rows of $\mathcal{S}$ are the left eigen-vectors and the columns of
$\mathcal{S}^{-1}$ are the right eigen-vectors. Multiplying Eq. (28) by
$\mathcal{S}$ one obtains the characteristics form of the system:

$l_{m}\partial_{t}\phi+\lambda_{m}l_{m}\partial_{x}\phi+l_{m}\mathcal{C}^{\prime}=0,$
(34)

which upon introduction of a new variable
$d\mathcal{V}_{m}=l_{m}d\phi+l_{m}\mathcal{C}^{\prime}dt$ reduces to a system
of wave equations with velocities $\lambda_{m}$:

$\partial_{t}\mathcal{V}_{m}+\lambda_{m}\partial_{x}\mathcal{V}_{m}=0.$ (35)

This analysis can be readily extended to multiple physical dimensions,
considering propagation along the $\alpha$-axis, by re-writing Eq. (28) as:

$\partial_{t}\phi+\mathcal{A}_{\alpha}\partial_{\alpha}\phi+\mathcal{C}^{\prime}=0,$
(36)

where $\alpha\in\\{x,y,z\\}$.

####  III.2.2 Characteristics: 1-D system

Let us first apply this analysis to the system of hyperbolic conservation
equations recovered by the D1Q3 lattice, shown in Eqs. (21a) and (21b). The
primitive form of this system will results in:

$\phi=\begin{bmatrix}\rho\\\ u_{x}\end{bmatrix}$ (37)

and:

$\mathcal{A}=\begin{bmatrix}u_{x}&\rho\\\
\varsigma^{2}+\frac{\partial_{\rho}P^{*}_{xx}}{\rho}&u_{x}+\frac{\partial_{u_{x}}P^{*}_{xx}}{\rho}\end{bmatrix}.$
(38)

Applying the characteristics analysis to this system the following eigen-
values are recovered:

$\displaystyle c_{s}^{+}$
$\displaystyle=u_{x}+\frac{\partial_{u_{x}}P^{*}_{xx}}{2\rho}+\sqrt{{\left(\frac{\partial_{u_{x}}P^{*}_{xx}}{2\rho}\right)}^{2}+{\varsigma}^{2}+\partial_{\rho}P^{*}_{xx}},$
(39a) $\displaystyle c_{s}^{-}$
$\displaystyle=u_{x}+\frac{\partial_{u_{x}}P^{*}_{xx}}{2\rho}-\sqrt{{\left(\frac{\partial_{u_{x}}P^{*}_{xx}}{2\rho}\right)}^{2}+{\varsigma}^{2}+\partial_{\rho}P^{*}_{xx}}.$
(39b)

Using the left eigen-vectors, readily obtained as right eigen-vectors of
$\mathcal{A}^{\dagger}$, the following wave system is recovered:

$\displaystyle\left[\frac{c_{s}^{+}-u_{x}}{\rho}\partial_{t}\rho+\partial_{t}u_{x}\right]+c_{s}^{+}\left[\frac{c_{s}^{+}-u_{x}}{\rho}\partial_{x}\rho+\partial_{x}u_{x}\right]$
$\displaystyle=0,$ (40a)
$\displaystyle\left[\frac{c_{s}^{-}-u_{x}}{\rho}\partial_{t}\rho+\partial_{t}u_{x}\right]+c_{s}^{-}\left[\frac{c_{s}^{-}-u_{x}}{\rho}\partial_{x}\rho+\partial_{x}u_{x}\right]$
$\displaystyle=0.$ (40b)

Note that for polynomial equilibria of Eqs. (19) and (20), i.e.
$P^{*}_{xx}=0$, one recovers, in accord with (26):

$\displaystyle c_{s}^{+}$ $\displaystyle=u_{x}+\varsigma,$ (41a)
$\displaystyle c_{s}^{-}$ $\displaystyle=u_{x}-\varsigma,$ (41b)

and

$\displaystyle\left[\frac{\varsigma}{\rho}\partial_{t}\rho+\partial_{t}u_{x}\right]+\left(u_{x}+\varsigma\right)\left[\frac{\varsigma}{\rho}\partial_{x}\rho+\partial_{x}u_{x}\right]$
$\displaystyle=0,$ (42a)
$\displaystyle\left[\frac{\varsigma}{\rho}\partial_{t}\rho-\partial_{t}u_{x}\right]+\left(u_{x}-\varsigma\right)\left[\frac{\varsigma}{\rho}\partial_{x}\rho-\partial_{x}u_{x}\right]$
$\displaystyle=0.$ (42b)

We will discuss the specific case of the entropic equilibrium in the next
paragraphs.

####  III.2.3 Characteristics: multi-dimensional systems

Now for systems of dimension larger than one, here to illustrate a 2-D system:

$\phi=\begin{bmatrix}\rho\\\ u_{x}\\\ u_{y}\end{bmatrix}$ (43)

and:

$\mathcal{A}_{x}=\begin{bmatrix}u_{x}&\rho&0\\\
\varsigma^{2}+\frac{\partial_{\rho}P^{*}_{xx}}{\rho}&u_{x}+\frac{\partial_{u_{x}}P^{*}_{xx}}{\rho}&\frac{\partial_{u_{y}}P^{*}_{xx}}{\rho}\\\
0&0&u_{x}\end{bmatrix}.$ (44)

and:

$\mathcal{A}_{y}=\begin{bmatrix}u_{y}&0&\rho\\\ 0&u_{y}&0\\\
\varsigma^{2}+\frac{\partial_{\rho}P^{*}_{y}}{\rho}&\frac{\partial_{u_{x}}P^{*}_{yy}}{\rho}&u_{y}+\frac{\partial_{u_{y}}P^{*}_{yy}}{\rho}\end{bmatrix}.$
(45)

and considering propagation along the $x-$axis one recovers the two acoustic
modes of Eqs. (40a) and (40b) and an additional _shear_ mode propagating at
speed $u_{x}$. Applying these eigen-values and the left side eigen-vectors the
following wave system is recovered:

$\displaystyle\partial_{t}u_{y}+u_{x}\partial_{x}u_{y}$ $\displaystyle=0,$
(46a)
$\displaystyle\left[\frac{c_{s}^{+}-u_{x}}{\rho}\partial_{t}\rho+\partial_{t}u_{x}\right]+c_{s}^{+}\left[\frac{c_{s}^{+}-u_{x}}{\rho}\partial_{x}\rho+\partial_{x}u_{x}\right]$
$\displaystyle=0,$ (46b)
$\displaystyle\left[\frac{c_{s}^{-}-u_{x}}{\rho}\partial_{t}\rho+\partial_{t}u_{x}\right]+c_{s}^{-}\left[\frac{c_{s}^{-}-u_{x}}{\rho}\partial_{x}\rho+\partial_{x}u_{x}\right]$
$\displaystyle=0.$ (46c)

Along the $y$-axis the wave system changes into:

$\displaystyle\partial_{t}u_{x}+u_{y}\partial_{y}u_{x}$ $\displaystyle=0,$
(47a)
$\displaystyle\left[\frac{c_{s}^{+}-u_{y}}{\rho}\partial_{t}\rho+\partial_{t}u_{y}\right]+c_{s}^{+}\left[\frac{c_{s}^{+}-u_{y}}{\rho}\partial_{y}\rho+\partial_{y}u_{y}\right]$
$\displaystyle=0,$ (47b)
$\displaystyle\left[\frac{c_{s}^{-}-u_{y}}{\rho}\partial_{t}\rho+\partial_{t}u_{y}\right]+c_{s}^{-}\left[\frac{c_{s}^{-}-u_{y}}{\rho}\partial_{y}\rho+\partial_{y}u_{y}\right]$
$\displaystyle=0.$ (47c)

where the acoustic propagation modes in $y-$direction are now:

$\displaystyle c_{s}^{+}$
$\displaystyle=u_{y}+\frac{\partial_{u_{y}}P^{*}_{yy}}{2\rho}+\sqrt{{\left(\frac{\partial_{u_{y}}P^{*}_{yy}}{2\rho}\right)}^{2}+{\varsigma}^{2}+\partial_{\rho}P^{*}_{yy}},$
(48a) $\displaystyle c_{s}^{-}$
$\displaystyle=u_{y}+\frac{\partial_{u_{y}}P^{*}_{yy}}{2\rho}-\sqrt{{\left(\frac{\partial_{u_{y}}P^{*}_{yy}}{2\rho}\right)}^{2}+{\varsigma}^{2}+\partial_{\rho}P^{*}_{yy}}.$
(48b)

Analysis of the propagation speed of different eigen-modes shows that acoustic
modes will propagate isotropically only in the limit of:

$\partial_{\rho}P_{xx}^{*}=\partial_{\rho}P_{yy}^{*},$ (49)

and

$\partial_{u_{x}}P_{xx}^{*}=\partial_{u_{y}}P_{yy}^{*}.$ (50)

A special solution of these conditions is $P_{xx}^{*}=P_{yy}^{*}=0$ which
corresponds to the polynomial equilibria explicitly enforcing the second-order
moment of the equilibrium distribution function, i.e. Eqs. (19) and (20).
However these equilibria come with limitations that the entropic equilibrium
does not have. This point will be discussed in the next paragraph.

###  III.3 Stabilization of entropic equilibrium via sound speed

Plugging the entropic equilibrium into the eigen-values obtained for the
general formulation two non-symmetrical sound propagation speeds are
recovered:

$\displaystyle c_{s}^{\rm e+}$
$\displaystyle=\frac{u_{x}+\varsigma\,\sqrt{2\,\sqrt{{\left(\frac{u_{x}}{\varsigma}\right)}^{2}+1}-1}}{\sqrt{{\left(\frac{u_{x}}{\varsigma}\right)}^{2}+1}},$
(51a) $\displaystyle c_{s}^{\rm e-}$
$\displaystyle=\frac{u_{x}-\varsigma\,\sqrt{2\,\sqrt{{\left(\frac{u_{x}}{\varsigma}\right)}^{2}+1}-1}}{\sqrt{{\left(\frac{u_{x}}{\varsigma}\right)}^{2}+1}}.$
(51b)

The propagation speed of these two eigen-modes are shown in Fig. 2 as a
function of speed.

Figure 2: (Left) Sound speed for entropic, Eqs. (51a) and (51b), and
polynomial equilibria, Eqs. (41a) and (41b), as a function of velocity
$u_{x}$. (Right) Comparison of the speed of fastest propagating eigen-modes:
(blue lines) polynomial and (red lines) entropic equilibria.

The behavior of entropic sound speed shows an interesting property of the
entropic model pointing already to a (potentially) unconditional linear
stability. Having re-written the equivalent Euler system in terms of coupled
wave equations, note that for the solver to access information required to
form the solution in time the numerical domain of dependence of any point in
space and time must include the analytical domain of dependence, i.e. the
initial conditions have an effect on the exact value of the solution at that
point [27]. Simply put, the fastest eigen-modes in the system _cannot_
propagate faster than the lattice links. Only considering three _physical_
eigen-modes, i.e. $u_{x}$, $c_{s}^{+}$ and $c_{s}^{-}$, one arrives at the
following condition on linear stability:

$\max(\lvert u_{x}\lvert,\lvert c_{s}^{+}\lvert,\lvert
c_{s}^{-}\lvert)\leq\frac{\delta r}{\delta t}.$ (52)

For the polynomial equilibria, given that sound speed is constant one recovers
the following maximum tolerable velocity:

$\lvert u_{x}^{\rm max}\lvert=\frac{\delta r}{\delta
t}-\varsigma=0.4226\frac{\delta r}{\delta t},$ (53)

which as will be shown in the next section through stability analyses, is
indeed the maximum reachable velocity. For the entropic equilibrium on the
other hand, it is observed that the sound speed self-adjusts as a function of
local velocity to guarantee Eq. (52) is always satisfied, leading to:

$\lvert u_{x}^{\rm max}\lvert=\frac{\delta r}{\delta t}.$ (54)

At the higher/lower end of the velocity spectrum,

$\displaystyle\lim_{u_{x}\rightarrow-\delta r/\delta t}c_{s}^{e+}$
$\displaystyle=\frac{\delta r}{\delta t}.$ (55a)
$\displaystyle\lim_{u_{x}\rightarrow-\delta r/\delta t}c_{s}^{e+}$
$\displaystyle=0,$ (55b) $\displaystyle\lim_{u_{x}\rightarrow\delta r/\delta
t}c_{s}^{e-}$ $\displaystyle=0.$ (55c)
$\displaystyle\lim_{u_{x}\rightarrow-\delta r/\delta t}c_{s}^{e-}$
$\displaystyle=-\frac{\delta r}{\delta t}.$ (55d)

Another point worth noting is that if one was to define the Mach number as the
ratio of flow speed to the normal mode propagation speed, i.e.

$\displaystyle{\rm Ma}^{+}$ $\displaystyle=\frac{u_{x}}{c_{s}^{+}-u_{x}},$
(56a) $\displaystyle{\rm Ma}^{-}$
$\displaystyle=\frac{u_{x}}{c_{s}^{-}-u_{x}},$ (56b)

the Mach number would evolve as shown in Fig. 3, which in effect means that
the entropic equilibrium leads to a scheme that is stable in the limit of
${\rm Ma}\rightarrow\infty$.

Figure 3: Mach number as a function of speed, $u_{x}$ as defined in Eqs. (56a)
and (56b) respectively in red and blue.

Note that the deviation in the propagation speed of acoustic modes scales with
the non-dimensional velocity with order three. This behavior is illustrated in
Fig. 4.

Figure 4: Normalized deviation of entropic equilibrium sound speed as a
function of non-dimensional velocity.

Overall the following points could be made based on the Euler-level analyses:

  * •

The entropic equilibrium pressure tensor admits deviations from the Maxwell-
Boltzmann pressure tensor that are not isotropic.

  * •

Deviations are negligible for non-dimensional velocities as large as $0.5$
which is well above the domain of validity of the weakly compressible flow
assumption.

  * •

The deviations scale as $\propto{\left(\frac{u\delta t}{\delta
r}\right)}^{4}$, see Fig. 1.

  * •

The entropic equilibrium leads to a non-isotropic sound speed.

  * •

In the limit of low non-dimensional flow velocities, this non-isotropic
deviation scales out with $\propto{\left(\frac{u\delta t}{\delta
r}\right)}^{3}$, see Fig. 4.

  * •

The _flow velocity-aware_ nature of the sound speed allows for a system that
adjusts its fastest propagating eigen-mode speed so as to always guarantee the
condition of Eq. (52) is satisfied, see Fig. 2.

  * •

Defining the Mach number as the ratio of propagation speed of shear to normal
modes, it is observed that entropic equilibrium can go to ${\rm
Ma}^{\pm}\rightarrow\infty$ in the limit of $u_{x}\rightarrow\pm\delta
r/\delta t$, see Fig. 3.

The slightly modified equilibrium pressure also affects dissipation at the
Navier-Stokes level. The altered dissipation behavior is discussed in the next
section.

##  IV Navier-Stokes-level dynamics: dissipation

###  IV.1 1-D lattice: Bulk viscosity

To illustrate the impact of the entropic equilibrium on dissipation let us
first start with a 1-D system with three discrete velocities. In 1-D there is
only acoustic modes, contrary to 2- and 3-D cases where shear modes are also
present. Performing the classical Chapman–Enskog analysis, at order
$\varepsilon^{2}$ the following momentum balance equation is recovered:

$\partial_{t}^{(2)}\rho
u_{x}+\partial_{x}\left(1-\beta\right)\Pi_{xx}^{(1)}=0,$ (57)

where after some algebra the non-equilibrium stress tensor $\Pi_{xx}^{(1)}$
for any equilibrium pressure of the form:

$P_{xx}=\rho\varsigma^{2}+P^{*}_{xx},$ (58)

results in:

$\Pi_{xx}^{(1)}=-\frac{1}{2\beta}\left[2A\rho\varsigma^{2}\partial_{x}u_{x}+B\partial_{x}\rho\right],$
(59)

with

$A=\left(1-\frac{3}{2}\frac{u_{x}^{2}}{\varsigma^{2}}-\frac{3}{2\rho\varsigma^{2}}u_{x}\partial_{u_{x}}P^{*}_{xx}-\frac{{(\partial_{u_{x}}P^{*}_{xx})}^{2}}{2\rho^{2}\varsigma^{2}}-\frac{\partial_{\rho}P^{*}_{xx}}{2\varsigma^{2}}\right),$
(60)

and

$B=-3u_{x}\partial_{\rho}P^{*}_{xx}-\frac{\varsigma^{2}}{\rho}\partial_{u_{x}}P^{*}_{xx}-\frac{\partial_{u_{x}}P^{*}_{xx}\partial_{\rho}P^{*}_{xx}}{\rho}-u_{x}^{3}.$
(61)

For the sake of clarity, in the remainder of the article, we will refer to
(60) as the effective viscosity and (61) as the compressibility error at the
Navier-Stokes level. For the case of the polynomial equilibria, i.e. both
second-order and product form, one gets:

$A=\left(1-\frac{3}{2}\frac{u_{x}^{2}}{\varsigma^{2}}\right),$ (62)

and:

$B=-u_{x}^{3},$ (63)

which, neglecting the second term, indicates that the corresponding partial
differential equation is only dissipative for:

$\lvert u_{x}\lvert\leq\varsigma\sqrt{\frac{2}{3}}.$ (64)

Note that this maximum velocity is larger than the one imposed by CFL
condition, i.e. Eq. (53). In the case of the entropic equilibrium:

$A=1-\frac{3}{2}\frac{u_{x}^{2}}{\varsigma^{2}}+\frac{{(u_{x}/\varsigma)}^{2}+3\,{(u_{x}/\varsigma)}^{4}-2\,\sqrt{{(u_{x}/\varsigma)}^{2}+1}+2}{2{(u_{x}/\varsigma)}^{2}+2}$
(65)

and:

$B=0.$ (66)

Clearly for the entropic equilibrium in the limit of $u_{x}\rightarrow\delta
r/\delta t$, the effective viscosity $A$ goes to zero. The effective
viscosities of Eqs. (62) and (65) are compared and shown in Fig. 5.

Figure 5: Comparison of the effective bulk viscosity: (blue) polynomial and
(black) entropic equilibria.

Furthermore, it is interesting to note that while the entropic equilibrium
guarantees a positive effective viscosity it maintains the second-order
convergence to the nominal viscosity, i.e. $A=1$, in the limit of $u_{x}\delta
t/\delta r\rightarrow 0$, as illustrated in Fig. 6.

Figure 6: Scaling of the error in effective 1-D viscosity, i.e. $1-A$: (blue)
polynomial and (black) entropic equilibria.

Going even further and comparing the compressibility errors, i.e. $B$, it is
observed that the entropic equilibrium, contrary to polynomial equilibrium,
does not have such errors. For the polynomial equilibria this error scales out
at order three.
While two types of error were shown to exist in the Navier-Stokes level
dissipation term, it is still unclear how each term would affect the
dissipation of the eigen-modes derived in section III. This issue will be
discussed in the next section.

###  IV.2 Analysis of linearized Navier-Stokes equations in the limit of
vanishing wave-number

In the previous section we observed that two types of errors affecting the
Navier-Stokes-level dissipation could be observed: One we called effective
viscosity and another one we introduced as compressiblity error. While the
role of the former is quite clear, even from its name, the _real_ effect of
the latter on viscosity is rather abstract. To clarify the role of each error
we will analyze the corresponding linearized equations in the limit of
vanishing wave-number, i.e. $k_{x}\rightarrow 0$. For the sake of clarity we
will first start with the simple case of polynomial equilibria. To linearize
the system of partial differential equations at the Euler+Navier-Stokes level
we separate velocity and density into a constant contribution and a
perturbation,

$\displaystyle\rho$ $\displaystyle=\bar{\rho}+\rho^{\prime},$ (67a)
$\displaystyle u_{x}$ $\displaystyle=\bar{u}_{x}+u_{x}^{\prime}.$ (67b)

Introducing these into the system of partial differential equations and only
keeping first-order perturbations we obtain the linearized system:

$\displaystyle\partial_{t}\rho^{\prime}$
$\displaystyle=-\bar{\rho}\partial_{x}u_{x}^{\prime}-\bar{u}_{x}\partial_{x}\rho^{\prime},$
(68a) $\displaystyle\partial_{t}u_{x}^{\prime}$
$\displaystyle=-\bar{u}_{x}\partial_{x}u_{x}^{\prime}-\frac{\varsigma^{2}}{\bar{\rho}}\partial_{x}\rho^{\prime}+2\nu
A\partial_{x}^{2}u_{x}^{\prime}+\frac{\nu
B}{\bar{\rho}^{2}\varsigma^{2}}\partial_{x}^{2}\rho^{\prime},$ (68b)

where $B=B(\bar{\rho},\bar{u}_{x})$ and $A=A(\bar{\rho},\bar{u}_{x})$. Note
the last term in Eq. (68b) stemming from the compressibility error. Next we
consider the perturbations to be monochromatic plane waves, i.e.

$\displaystyle\rho^{\prime}$
$\displaystyle=\hat{\rho}\exp{\left(\mathrm{i}(\omega t-k_{x}x)\right)},$
(69a) $\displaystyle u_{x}^{\prime}$
$\displaystyle=\hat{u}_{x}\exp{\left(\mathrm{i}(\omega t-k_{x}x)\right)},$
(69b)

where $\mathrm{i}=\sqrt{-1}$. Introducing the perturbations into Eqs. (68a)
and (68b) we end up with the following system:

$\displaystyle\omega\hat{\rho}$
$\displaystyle=\mathrm{i}\bar{\rho}k_{x}\hat{u}_{x}+\sqrt{-1}k_{x}\bar{u}_{x}\hat{\rho},$
(70a) $\displaystyle\omega\hat{u}_{x}$
$\displaystyle=\mathrm{i}k_{x}\bar{u}_{x}\hat{u}_{x}+\frac{\varsigma^{2}}{\bar{\rho}}k_{x}\hat{\rho}-2\nu
Ak_{x}^{2}\hat{u}_{x}-\frac{\nu
B}{\bar{\rho}^{2}\varsigma^{2}}k_{x}^{2}\hat{\rho}.$ (70b)

Finally, we solve the corresponding eigen-value problem and Taylor-expand
corresponding solutions around $k_{x}=0$; Keeping terms of order $k_{x}$ and
$k_{x}^{2}$ we get the following eigen-values:

$\displaystyle\omega_{c_{s}^{+}}$
$\displaystyle=(\bar{u}_{x}+\varsigma)k_{x}+\mathrm{i}\nu\left(A+\frac{B}{2\bar{\rho}\varsigma^{3}}\right)k_{x}^{2},$
(71a) $\displaystyle\omega_{c_{s}^{-}}$
$\displaystyle=(\bar{u}_{x}-\varsigma)k_{x}+\mathrm{i}\nu\left(A-\frac{B}{2\bar{\rho}\varsigma^{3}}\right)k_{x}^{2}.$
(71b)

Looking at the dissipation terms, i.e. $\propto k_{x}^{2}$, interesting
observations on the role of each type of error could be made: (a) The
existence of the compressibility error results in anisotropy in the
dissipation of acoustic modes causing $c_{s}^{+}$ and $c_{s}^{-}$ to dissipate
at different rates even though the dispersion is isotropic for polynomial
equilibria. (b) The average of the dissipation rates of the two acoustic modes
is modulated by the effective viscosity coefficient, $A$.
Extending this analysis to the more general form of the equilibrium as in
previous sections the following eigen-values are recovered:

$\displaystyle\omega_{c_{s}^{+}}$ $\displaystyle=c_{s}^{+}k_{x}$
$\displaystyle+\mathrm{i}\nu
A\left(1+\frac{\varsigma^{2}\partial_{\bar{u}_{x}}P_{xx}^{*}+B/A}{2\bar{\rho}\varsigma^{2}\sqrt{\varsigma^{2}+\partial_{\bar{\rho}}P_{xx}^{*}+(\partial_{\bar{u}_{x}}P_{xx}^{*}/2\bar{\rho})^{2}}}\right)k_{x}^{2},$
(72a) $\displaystyle\omega_{c_{s}^{-}}$ $\displaystyle=c_{s}^{-}k_{x}$
$\displaystyle+\mathrm{i}\nu
A\left(1-\frac{\varsigma^{2}\partial_{\bar{u}_{x}}P_{xx}^{*}+B/A}{2\bar{\rho}\sqrt{\varsigma^{2}+\partial_{\bar{\rho}}P_{xx}^{*}+(\partial_{\bar{u}_{x}}P_{xx}^{*}/2\bar{\rho})^{2}}}\right)k_{x}^{2},$
(72b)

where $c_{s}^{+}$ and $c_{s}^{-}$ are those derived in Eqs. (48a) and (48b)
and $P_{xx}^{*}=P_{xx}^{*}(\bar{\rho},\bar{u}_{x})$. Plugging in the entropic
equilibrium:

$\displaystyle\omega_{c_{s}^{+}}$ $\displaystyle=c_{s}^{e+}k_{x}$
$\displaystyle+\mathrm{i}\nu
A\left(1+\frac{1-\varsigma\bar{u}_{x}\left(\sqrt{{\left(\frac{\bar{u}_{x}}{\varsigma}\right)}^{2}+1}-1\right)}{\sqrt{2\sqrt{{\left(\frac{\bar{u}_{x}}{\varsigma}\right)}^{2}+1}-1}}\right)k_{x}^{2},$
(73a) $\displaystyle\omega_{c_{s}^{-}}$ $\displaystyle=c_{s}^{e-}k_{x}$
$\displaystyle+\mathrm{i}\nu
A\left(1-\frac{1-\varsigma\bar{u}_{x}\left(\sqrt{{\left(\frac{\bar{u}_{x}}{\varsigma}\right)}^{2}+1}-1\right)}{\sqrt{2\sqrt{{\left(\frac{\bar{u}_{x}}{\varsigma}\right)}^{2}+1}-1}}\right)k_{x}^{2},$
(73b)

where $c_{s}^{e+}$ and $c_{s}^{e-}$ were derived in Eqs. (51a) and (51b). The
interesting point worth noting here is that even though for the entropic
equilibrium it was shown that $B=0$, writing the equations in terms of the
eigen-modes of the system one recovers a compressibility-like error term.
Now that we have a closed-form solution for the total dissipation rate of the
linearized hydrodynamic equations, noted as $A^{\prime}$, defined in Eqs.
(73a) and (73b) as:

$\omega_{c_{s}}=c_{s}^{e}k_{x}+\mathrm{i}\nu A{{}^{\prime}}k_{x}^{2},$ (74)

well-posedeness of the corresponding wave system can be assessed by looking at
the sign of the dissipation rates. The results are shown in Fig. 7.

Figure 7: Overall dissipation rate of linearized hydrodynamic equations for
(left) polynomial and (right) entropic equilibria. Red lines show dissipation
of $c_{s}^{+}$ mode and blue lines that of $c_{s}^{-}$ mode.

It is observed that the entropic equilibrium always remains dissipating while
the polynomial equilibrium keeps that property within a certain range of
velocity:

$\frac{1}{\sqrt{3}}-1\leq\frac{u_{x}\delta t}{\delta r}\leq
1-\frac{1}{\sqrt{3}},$ (75)

which interestingly enough matches exactly the CFL condition on eigen-modes as
shown in Fig. 2.
The next section will discuss the dissipation rates of shear modes.

###  IV.3 Multi-dimensional lattices: Shear viscosity

The multiscale analysis of the previous section can be extended to multiple
dimensions, detailed in B. Decomposing the non-equilibrium stress tensor into
two contributions as:

$\Pi_{2}^{(1)}=-\frac{1}{2\beta}\left(\bm{S}^{*}+\bm{D}^{*}\right),$ (76)

where:

$\bm{S}^{*}={\left[(\rho\varsigma^{2}\bm{I}+\bm{P}^{*})\cdot\bm{\nabla}\bm{u}\right]}+{\left[(\rho\varsigma^{2}\bm{I}+\bm{P}^{*})\cdot\bm{\nabla}\bm{u}\right]}^{\dagger}+\bm{E}^{s}.$
(77)

Here $\bm{E}^{s}$ regroups errors specific to the second-order polynomial,
detailed in Eq. (126). Focusing on off-diagonal components and following the
presentation of the previous section, let us introduce an effective viscosity,
$\bm{A}^{s}$ such that:

$S^{*}_{\alpha\beta}+E_{\alpha\beta}^{s}=\rho\varsigma^{2}\left[A_{\alpha\beta}^{s}\partial_{\alpha}u_{\beta}+A_{\beta\alpha}^{s}\partial_{\beta}u_{\alpha}\right]+B_{\alpha}^{s}\partial_{\alpha}\rho+B_{\beta}^{s}\partial_{\beta}\rho.$
(78)

In the case of the product form equilibrium:

$A_{\alpha\beta}^{s}=A_{\beta\alpha}^{s}=1,$ (79)

while for the second order polynomial:

$A_{\alpha\beta}^{s}=1-\frac{u_{\alpha}}{\varsigma^{2}}(u_{\alpha}+2u_{\beta}),$
(80)

and the entropic equilibrium, Eq. (18):

$A_{\alpha\beta}^{s}=1+\frac{P_{\alpha\alpha}^{*}}{\rho\varsigma^{2}}=2\sqrt{{\left(\frac{u_{\alpha}}{\varsigma}\right)}^{2}+1}-{\left(\frac{u_{\alpha}}{\varsigma}\right)}^{2}-1.$
(81)

The dependence of $A_{xy}^{s}$ on $u_{x}\delta t/\delta r$ with $u_{y}=0$ for
all equilibria are illustrated in Fig. 8.

Figure 8: Comparison of scaling of the effective viscosity in shear components
of the viscous stress tensor for the (black) entropic and (blue) second-order
polynomial equilibria.

The deviations in the entropic equilibrium are shown to scale out at order
four while for the second order polynomial equilibrium the order is only two.
Furthermore, the effective viscosity $A_{xy}^{s}$ for the entropic equilibrium
is positive definite for $\lvert u_{x}\lvert\in[0,\delta r/\delta t]$, while
for the second order polynomial equilibrium it is only positive for:

$\lvert u_{x}\lvert\leq\varsigma.$ (82)

This point is illustrated in Fig. 9.

Figure 9: Comparison of positivity domain of the effective viscosity in shear
components of the viscous stress tensor for the (black) entropic and (blue)
second-order polynomial equilibria.

Based on the analysis of the Navier-Stokes level equations the following
observations could be made:

  * •

Contrary to polynomial equilibria, the Navier-Stokes-level viscosity of normal
modes recovered by the entropic discrete equilibrium is positive definite, for
$-\delta r/\delta t\leq u_{x}\leq\delta r/\delta t$, see Fig. 5. Polynomial
equilibria lose positivity following Eq. (64).

  * •

The deviation in the viscosity tied to compressibility, i.e.
$\propto\partial_{x}\rho$, for the polynomial equilibria scales with
$\propto{\left(\frac{u_{x}\delta t}{\delta r}\right)}^{3}$, see Eq. (63). The
entropic equilibrium on the other hand does not admit such errors, see Eq.
(66).

  * •

An analysis of the dissipation of eigen modes of different models reveals that
entropic normal eigen modes are indeed subject to a compressibility-like
error, see Eqs. (73a) and (73b).

  * •

Compressibility errors cause the two acoustic modes $c_{s}^{+}$ and
$c_{s}^{-}$ to dissipate at different rates, even in the limit of
$k_{x}\rightarrow 0$.

  * •

Following the trend set in previous sections, the polynomial equilibria result
in an overall dissipation rate that is only positive for both modes in $\lvert
u_{x}\lvert\in[0\text{ }\frac{\delta r}{\delta t}-\varsigma]$. For the
entropic equilibrium on the other hand the overall dissipation rate of the
linearized hydrodynamic equations is positive definite for all modes in
$\lvert u_{x}\lvert\in[0\text{ }\delta r/\delta t]$.

  * •

Looking at Fig. 7, it is seen that the anisotropy in dissipation rate of the
two acoustic modes is much less pronounced for the entropic equilibrium.

  * •

As for normal modes, the Navier-Stokes level dissipation rate of shear modes
in the entropic equilibrium is positive definite. The second-order polynomial
equilibrium on the other hand is only positive under the condition of Eq.
(82). Note that for a velocity vector with non-zero components in the $y-$ and
$z-$directions this condition becomes more restrictive.

  * •

The deviations in the shear dissipation rate scale out with
$\propto{\left(\frac{u_{\alpha}\delta t}{\delta r}\right)}^{4}$ while for the
second order polynomial equilibrium they scale as
$\propto{\left(\frac{u_{\alpha}\delta t}{\delta r}\right)}^{2}$.

##  V Numerical applications

###  V.1 Linear stability analysis

####  V.1.1 von Neumann spectral analysis: introduction

In the context of the von Neumann linear analysis, the discrete time-evolution
system of equations is expanded around a reference state
$\bar{f}_{i}\left(\bar{\rho},\bar{u}\right)$ via a first-order Taylor
expansion:

$f_{i}\approx\bar{f}_{i}+f^{{}^{\prime}}_{i},$ (83)

where $f^{{}^{\prime}}_{i}$ is the linear perturbation. Defining the discrete
collision operator as:

$\Omega_{i}=2\beta\left(f_{i}^{\rm eq}-f_{i}\right),$ (84)

the linearized operator is:

$\Omega_{i}(f_{i})\approx\Omega_{i}\rvert_{\bar{f}_{i}}+\mathcal{J}_{ij}f^{{}^{\prime}}_{j},$
(85)

where $\mathcal{J}_{ij}$ is the Jacobian of the collision operator evaluated
about $\bar{f}_{j}$, i.e

$\mathcal{J}_{ij}=\partial_{f_{j}}\Omega_{i}\rvert_{\bar{f}_{j}}.$ (86)

Placing back these expressions into the discrete time-evolution equations one
obtains [28]:

$f_{i}^{{}^{\prime}}\left(\bm{r}+\bm{c}_{i}\delta t,t+\delta
t\right)=\left(\delta_{ij}+\mathcal{J}_{ij}\right)f_{j}^{{}^{\prime}}\left(\bm{r},t\right).$
(87)

Detailed expression for the Jacobian of the entropic equilibrium is given in
Appendix C. The last step of the von Neumann analysis is to assume that
perturbations $f^{\prime}_{i}$ are monochromatic plane waves of the form:

$f^{\prime}_{i}=F_{i}\exp{[\mathrm{i}(\bm{k}\cdot\bm{r}-\omega_{i}t)]},$ (88)

where $F_{i}$ is the wave amplitude, $\mathrm{i}$ is the imaginary unit,
$||\bm{k}||=k$ is the wave-number, and $\omega$ is the complex time frequency
of the wave. $k$ is related to the wave-length of $f^{\prime}_{i}$, whereas
$\Im(\omega)$ and $\Re(\omega)$ are related to its attenuation and propagation
speed. Introducing these perturbations into Eq. (87) one obtains the following
eigenvalue problem of size $Q$ :

$\bm{MF}=\exp{(-\mathrm{i}\omega_{i})}\bm{F},$ (89)

where $\bm{F}$ is the eigenvector composed of all amplitudes. $\bm{M}$ is the
matrix associated to Eq. (87). This matrix can be expressed as :

$\bm{M}=\bm{E}\left[\bm{\delta}+\mathcal{J}\right],$ (90)

with

$E_{ij}=\exp[-\mathrm{i}(\delta t\bm{c}_{i}\cdot\bm{k})]\delta_{ij}.$ (91)

Note that the matrix $\bm{M}$ and the eigenvalue problem (89) depend on the
mean flow ($\bar{\rho},\lvert\lvert\bar{\bm{u}}\lvert\lvert$), the wave-number
($k_{x}$ and $k_{y}$) and the relaxation frequency $\beta$. This means that
for each set of these parameters the eigenvalue problem needs to be solved to
obtain the corresponding values of $\Re{(\omega)}$ and $\Im{(\omega)}$. Doing
so, the spectral properties can be obtained.

####  V.1.2 Linear stability domain

To assess the linear stability domain via the von Neumann approach, for a
given set of properties
$(\bar{\rho},\lvert\lvert\bar{\bm{u}}\lvert\lvert,\nu\delta t/\delta r^{2})$
the eigen-value problem of Eq. (89) is solved for $k_{x}\in[0\text{ }\pi]$ and
$k_{y}\in[0\text{ }\pi]$ and
$\theta=\cos^{-1}\left(\frac{u_{x}}{\lvert\lvert\bm{u}\lvert\lvert}\right)\in[0\text{
}\pi/2]$. The solver is considered to be linear-stable for the set of
parameters if it guarantees negative attenuation rate for all modes and values
of $k_{x},k_{y}$ and $\theta$. We first quantify the stability domain of the
discrete solvers by finding the maximum non-dimensional velocity for which
linear stability is guaranteed for a range of non-dimensional viscosities.
Here we consider $\nu\delta t/\delta r^{2}\in\\{10^{-6},5\times
10^{-6},10^{-5},5\times 10^{-5},10^{-4},5\times 10^{-4},10^{-3},5\times
10^{-3},10^{-2},5\times 10^{-2},0.1,0.5,1,1.1,1.2\\}$. The stability domains
of three different equilibria, namely entropic, second order polynomial and
product form are shown in Fig. 10.

Figure 10: Comparison of linear stability domain of different equilibrium
distribution functions: (red with diamond markers) second order polynomial,
(blue with square markers) product form and (black with circular markers)
entropic.

It can be noted that the entropic equilibrium had, by far, the widest domain
of stability guaranteeing unconditional stability which is unaffected by the
value of the non-dimensional viscosity [29]. For the second-order polynomial
and product form equilibria on the other hand, in the limit of vanishing non-
dimensional viscosities the solvers become unconditionally unstable. On the
other end of the spectrum, i.e. large non-dimensional viscosities the maximum
speed encounters a threshold which is exactly the one shown in Eq. (53).
To better illustrate the stability behavior of the different equilibria, and
especially to better show the unconditional linear stability of the entropic
equilibrium, directional (orientation of velocity vector) stability domains
are shown in Figs. 11, 12 and 13.

Figure 11: Directional stability domain of the entropic equilibrium for two
different non-dimensional viscosities: (dashed black) $10^{-5}$and (dotted
red) $0.1$.

For the entropic equilibrium, the solver linear stability covers the entire
lattice and confirms observations from characteristics analysis, i.e.
unconditional stability.

Figure 12: Directional stability domain of the polynomial equilibrium for two
different non-dimensional viscosities: (dashed black) $10^{-5}$and (dotted
red) $0.1$.

For the second-order polynomial for large non-dimensional viscosities the
linear stability domains is close to isotropic and limited by the condition of
Eq. (53).

Figure 13: Directional stability domain of the product form equilibrium for
two different non-dimensional viscosities: (dashed black) $10^{-5}$and (dotted
red) $0.1$.

The product form equilibrium exhibits a directional behavior similar to the
entropic one for large non-dimensional viscosities, however with a much
smaller tolerated maximum velocity. The largest reachable non-dimensional
velocity is reached when the velocity vector is oriented in the direction of
one of the diagonal links of the lattice,

$\lvert\lvert\bm{u}\lvert\lvert=\sqrt{2}\left(\frac{\delta r}{\delta
t}-\varsigma\right)=0.598\frac{\delta r}{\delta t}.$ (92)

####  V.1.3 Spectral dispersion/dissipation

In this section we mainly aim at confirming two observations made via
theoretical analyses of the hydrodynamic equations: (a) local-velocity
awareness of the entropic sound speed and (b) the effect of different types of
error on the effective normal dissipation rate. To that end we will compare
results from von Neumann analysis of the full discrete system to previously
obtained analytical expressions.
First we look into the dispersion of normal modes. The spectral propagation
speed is computed as $\Im(\omega)/k_{x}$. The calculations have been carried
out for velocities $u_{x}\delta r/\delta t\in[-1\,\,1]$ and are compared to
hydrodynamic limit predictions from the characteristics analysis. The results
are shown in Fig. 14.

Figure 14: Speed of normal modes for the entropic equilibrium as a function of
velocity, (black markers) as obtained from the spectral analysis of the full
discrete system and (red line) from the characteristics analysis, i.e. Eqs.
(51a) and (51b). The insert on the right hand side shows (black) the typical
discrete spectrum. Red lines show dissipation of $c_{s}^{+}$ mode and blue
lines that of $c_{s}^{-}$ mode.

The comparison shows that the analytical expressions derived in Eqs. (51a) and
(51b) exactly match the numerical results from spectral analysis at
$k_{x}\rightarrow 0$.
The second parameter to be validate in this section is the effective
dissipation rates of normal modes. The dissipation rate is extracted as
$-\Re(\omega)\delta t/k_{x}^{2}\delta r^{2}$ from the full spectral analysis
and compared to analytical expressions derived in Eqs. (73a) and (73b). The
calculations have been carried out for velocities $u_{x}\delta r/\delta
t\in[-1\,\,1]$ and are compared to hydrodynamic limit predictions in Fig. 15.

Figure 15: Effective overall viscosity, $A^{\prime}$, for the entropic
equilibrium as a function of velocity, (black markers) as obtained from the
spectral analysis of the full discrete system and (red line) from the
corresponding Navier-Stokes linearized equations, i.e. Eqs. (73a) and (73b).
The insert on the right hand side shows (black) the typical discrete spectrum.
Red lines show dissipation of $c_{s}^{+}$ mode and blue lines that of
$c_{s}^{-}$ mode.

The comparisons show exact agreement and confirm the analytical results
derived in previous sections.

###  V.2 Correction of leading-order bulk viscosity error

####  V.2.1 Rescaling bulk viscosity

In the previous section, it was shown that at the Navier-Stokes level the bulk
viscosity had two Galilean-variant errors, see Eq. (60) and (61). The first
term, i.e. Eq. (60) can be corrected by redefining the relaxation frequency
as:

$\beta=\frac{\delta t}{\frac{2\nu}{A\varsigma^{2}}+\delta t},$ (93)

which, for instance, for polynomial equilibria reduces to:

$\beta=\frac{\delta
t}{\frac{2\nu}{\left(\varsigma^{2}-3u_{x}^{2}/2\right)}+\delta t}.$ (94)

In doing so, and before checking the effect of this correction two points have
to be considered: (a) The effect of this redefinition of the relaxation rate
on the second error term, i.e. Eq. (61), and (b) changes in the now-velocity-
dependent and variable relaxation rate. The changes induced in the relaxation
rate by this correction are illustrated in Fig. 16.

Figure 16: Variation in the redefined relaxation rate, Eq. (93) as a function
of local velocity: (blue) polynomial equilibrium and (black) entropic
equilibrium. Here non-dimensional viscosity is taken to be $\nu\delta t/\delta
r^{2}=0.1$.

The key observation here is that for the polynomial equilibria the redefined
relaxation rate grows much faster than for the entropic equilibrium, and that
$\beta\in[0\text{ }1]$ is not always guaranteed. For this condition to be
satisfied one must have:

$\frac{2\nu}{\varsigma^{2}A}\geq 0,$ (95)

which boils down to a condition of positivity of $A$ since both $\nu$ and
$\varsigma$ are positive. As shown previously in Fig. 5 this condition is not
always satisfied for the polynomial and product form equilibria. For the
entropic equilibrium on the other hand, because $A$ is positive definite, the
redefined relaxation rate satisfies $\beta\in[0\text{ }1]$ for $\lvert
u_{x}\lvert<\delta r/\delta t$.
The effect of this correction will be checked via numerical simulations in the
next section.

####  V.2.2 Validation via dissipation of isothermal pressure waves

In the limit of small velocities and density variations, the contributions
from non-linear terms in the Navier-Stokes equation can be ignored and
acoustic waves can be modeled through the linear theory including losses,
which reads [30]:

$\partial_{t}^{2}u_{x}=c_{s}^{2}\partial_{x}^{2}u_{x}+\left(\frac{D+1}{D}\nu+\frac{\eta}{\rho}\right)\partial_{t}\partial_{x}^{2}u_{x},$
(96)

where the bulk viscosity $\eta$ in the context of the classical LB model is
fixed at $\frac{D-1}{D}\rho\nu$. The exact solution of this equation can be
shown to be of the form $u_{x}\propto\exp\left({\sqrt{-1}k_{x}x+\sigma
t}\right)$ where:

$\sigma=-\left(\frac{2}{3}\nu+\frac{1}{2}\eta\right)k_{x}^{2}+\sqrt{-1}k_{x}c_{s}\sqrt{1-{\left(\frac{D+1}{D}\nu+\frac{\eta}{\rho}\right)}^{2}\frac{k_{x}^{2}}{c_{s}^{2}}}.$
(97)

The total energy, here defined as
$E=\int_{x}\frac{1}{2}u^{2}+c_{s}^{2}\left(\rho-\rho_{0}\right)dx$ will
therefore decay in time as:

$E(t)=E(0)\exp{\left[-\left(\frac{D+1}{D}\nu+\frac{\eta}{\rho}\right)\frac{t}{k_{x}^{2}}\right]},$
(98)

with $k_{x}=\frac{2\pi}{L}$, where $L$ is the wavelength. To evaluate the
effective bulk viscosity of the considered numerical schemes, we initialize a
wave-function in a domain of size $N_{x}$ with periodic boundary conditions
as:

$\rho=\rho_{0}+\delta\rho\sin\left(\frac{2\pi x}{N_{x}}\right),$ (99)

for different values of initial velocity $U_{0}\in[-0.3\text{ }0.3]$. The
waves are then left to evolve, and energy is stored for each time-step. The
resulting time-evolution of the energy is then fitted with a function of the
form of (98) to extract the effective viscosity. This process is illustrated
in Fig. 17, with results obtained from a sample simulation.

Figure 17: Evolution of total energy as a function of time as obtained from
two simulations using the entropic equilibrium, one with the regular
definition of $\beta$ from Eq. (2) and one with the redefined $\beta$ from Eq.
(93). $\nu\delta t/\delta r^{2}$ is set to $0.05$.

Simulations were ran for different convection speeds for both entropic and
product form equilibria both with and without correction of the relaxation
frequencies. For the entropic equilibrium the corrected relaxation frequencies
were computed according to Eq. (93) while for the product form Eq. (94) was
used.

Figure 18: Effective viscosity as a function of non-dimensional velocity as
obtained from the 1-D acoustic wave dissipation case for (left) product form
and (right) entropic equilibria. Red square symbols show results with $\beta$
defined as in Eq. (2) and black circular symbols results with $\beta$ as
defined in Eq. (93).

The obtained results are illustrated in Fig. 18. Based on these results a
number of observations can be made; The first one is that in the non-corrected
forms, the product form equilibrium induces more deviations into the effective
viscosity as compared to the entropic equilibrium. Second, the rescaling of
the relaxation frequency does reduce the error in the effective viscosity,
however as expected it does not fully eliminate it.

Finally the outcomes of this final section can be summarized as follows:

  * •

Linear spectral analysis confirmed that the entropic equilibrium guarantees
unconditional stability for any wave-number $\bm{k}$ for $\lvert
u_{\alpha}\lvert\leq\delta r/\delta t$, see Fig. 11. The family of polynomial
equilibria, both second order and product form have a more limited linear
stability domain and are affected by the choice of the non-dimensional
viscosity, see Fig. 12, 13 and 10.

  * •

The errors in effective viscosity, in 1-D, from Eq. (60) can be, to a great
extend, reduced via a simple rescaling of the relaxation frequency as given by
Eq. (93). This is confirmed by Fig. 18.

  * •

The rescaling strategy applied to the entropic equilibrium has a number of
advantages over its product form counter part: For the product form it can
result in $\beta\notin[0\text{ }1]$, see Eq. (95). The corrected entropic
model exhibits much less pronounced deviations as compared to the product
form.

##  VI Conclusions and discussion

Construction of discrete distribution functions in reduced kinetic models,
whether it be moments-based methods such as the Grad system or discrete
velocity methods such as the lattice Boltzmann method is crucial to both
accuracy and stability of the system. While the classical approach, used for
instance for classical polynomial equilibria in the lattice Boltzmann,
consists in matching moments of interest in the equilibrium to their
continuous counter-parts. The entropic construction provides an approach that
satisfies target moments constraints within the range of operation and
asymptotically vanishing sound speed which, as demonstrated in details here,
leads to unconditional linear stability. Our detailed analysis of the entropic
construction further showed that this asymptotic regularization, in the case
of the entropic equilibrium, also leads to positive-definite Navier-Stokes-
level dissipation rates. A simple local approach was also introduced to
further reduce deviation in the effective viscosity of the model. The
correction was also shown to maintain positive-definiteness only for the
entropic equilibrium.

###### Acknowledgements.

This work was supported by European Research Council (ERC) Advanced Grant
834763-PonD. Computational resources at the Swiss National Super Computing
Center CSCS were provided under the grant s1066.

## References

  * Succi [2002] S. Succi,  _The Lattice Boltzmann Equation for Fluid Dynamics and Beyond_ (Clarendon, 2002).
  * Bhatnagar  _et al._ [1954] P. L. Bhatnagar, E. P. Gross, and M. Krook, A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems, Physical Review 94, 511 (1954).
  * Kruger  _et al._ [2017] T. Kruger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M. Viggen,  _The Lattice Boltzmann Method: Principles and Practice_, Graduate Texts in Physics (Springer International Publishing, Cham, 2017).
  * Smith [1985] G. D. Smith,  _Numerical solution of partial differential equations: finite difference methods_ , 3rd ed., Oxford applied mathematics and computing science series (Clarendon Press ; Oxford University Press, Oxford [Oxfordshire] : New York, 1985).
  * Strikwerda [2004] J. C. Strikwerda,  _Finite difference schemes and partial differential equations_ , 2nd ed. (Society for Industrial and Applied Mathematics, Philadelphia, 2004).
  * Lax and Richtmyer [1956] P. D. Lax and R. D. Richtmyer, Survey of the stability of linear finite difference equations, Communications on Pure and Applied Mathematics 9, 267 (1956).
  * Latt and Chopard [2006] J. Latt and B. Chopard, Lattice Boltzmann method with regularized pre-collision distribution functions, Mathematics and Computers in Simulation 72, 165 (2006), iSBN: 0378-4754 Publisher: Elsevier.
  * Lallemand and Luo [2003] P. Lallemand and L.-S. Luo, Theory of the lattice Boltzmann method: Acoustic and thermal properties in two and three dimensions, Physical Review E 68, 036706 (2003).
  * Geier  _et al._ [2006] M. Geier, A. Greiner, and J. G. Korvink, Cascaded digital lattice Boltzmann automata for high Reynolds number flow, Physical Review E 73, 066705 (2006), publisher: APS.
  * Geier  _et al._ [2015] M. Geier, M. Schönherr, A. Pasquali, and M. Krafczyk, The cumulant lattice Boltzmann equation in three dimensions: Theory and validation, Computers & Mathematics with Applications 70, 507 (2015).
  * Karlin  _et al._ [1998] I. V. Karlin, A. N. Gorban, S. Succi, and V. Boffi, Maximum Entropy Principle for Lattice Kinetic Equations, Physical Review Letters 81, 6 (1998).
  * Mazloomi  _et al._ [2015] A. Mazloomi, S. Chikatamarla, and I. Karlin, Entropic Lattice Boltzmann Method for Multiphase Flows, Physical review Letters 114, 174502 (2015).
  * Chikatamarla  _et al._ [2006] S. S. Chikatamarla, S. Ansumali, and I. V. Karlin, Entropic lattice Boltzmann models for hydrodynamics in three dimensions, Physical review letters 97, 010201 (2006), publisher: APS.
  * Prasianakis and Karlin [2008] N. I. Prasianakis and I. V. Karlin, Lattice Boltzmann method for simulation of compressible flows on standard lattices, Physical Review E 78, 016704 (2008).
  * Yong and Luo [2003] W.-A. Yong and L.-S. Luo, Nonexistence of H theorems for the athermal lattice Boltzmann models with polynomial equilibria, Physical Review E 67, 051105 (2003).
  * Luo  _et al._ [2011] L.-S. Luo, W. Liao, X. Chen, Y. Peng, and W. Zhang, Numerics of the lattice Boltzmann method: Effects of collision models on the lattice Boltzmann simulations, Physical Review E 83, 056710 (2011).
  * Karlin  _et al._ [2011] I. V. Karlin, S. Succi, and S. S. Chikatamarla, Comment on “Numerics of the lattice Boltzmann method: Effects of collision models on the lattice Boltzmann simulations”, Physical Review E 84, 068701 (2011).
  * Luo [2012] L.-S. Luo, Reply to “Comment on ‘Numerics of the lattice Boltzmann method: Effects of collision models on the lattice Boltzmann simulations’”, Physical Review E 86, 048701 (2012).
  * Ansumali and Karlin [2000] S. Ansumali and I. V. Karlin, Stabilization of the lattice Boltzmann method by the H theorem: A numerical test, Physical Review E 62, 7999 (2000), publisher: APS.
  * Ansumali and Karlin [2002] S. Ansumali and I. V. Karlin, Entropy Function Approach to the Lattice Boltzmann Method, Journal of Statistical Physics 107, 291 (2002).
  * Karlin  _et al._ [1999] I. V. Karlin, A. Ferrante, and H. C. Öttinger, Perfect entropy functions of the Lattice Boltzmann method, Europhysics Letters (EPL) 47, 182 (1999).
  * Ansumali  _et al._ [2003] S. Ansumali, I. V. Karlin, and H. C. Öttinger, Minimal entropic kinetic models for hydrodynamics, EPL (Europhysics Letters) 63, 798 (2003), iSBN: 0295-5075 Publisher: IOP Publishing.
  * Prasianakis and Karlin [2007] N. I. Prasianakis and I. V. Karlin, Lattice Boltzmann method for thermal flow simulation on standard lattices, Physical Review E 76, 016702 (2007).
  * Shan and He [1998] X. Shan and X. He, Discretization of the Velocity Space in the Solution of the Boltzmann Equation, Physical Review Letters 80, 65 (1998).
  * Hosseini  _et al._ [2020] S. A. Hosseini, N. Darabiha, and D. Thévenin, Compressibility in lattice Boltzmann on standard stencils: effects of deviation from reference temperature, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 378, 20190399 (2020).
  * Thompson [1987] K. W. Thompson, Time dependent boundary conditions for hyperbolic systems, Journal of Computational Physics 68, 1 (1987).
  * Courant  _et al._ [1928] R. Courant, K. Friedrichs, and H. Lewy, Uber die partiellen Differenzengleichungen der mathematischen Physik, Mathematische Annalen 100, 32 (1928).
  * Worthing  _et al._ [1997] R. A. Worthing, J. Mozer, and G. Seeley, Stability of lattice Boltzmann methods in hydrodynamic regimes, Physical Review E 56, 2243 (1997).
  * Hosseini [2020] S. A. Hosseini,  _Development of a lattice Boltzmann-based numerical method for the simulation of reacting flows_ , Ph.D. thesis, Universite Paris-Saclay/Otto-von-Guericke-Universitat Magdeburg (2020).
  * Kinsler [2000] L. E. Kinsler, ed.,  _Fundamentals of acoustics_ , 4th ed. (Wiley, New York, 2000).

##  Appendix A One-dimensional multi-scale analysis: Bulk viscosity

Using a Taylor expansion around $(r,t)$:

$\delta t\mathcal{D}_{t}f_{i}+\frac{{\delta
t}^{2}}{2}{\mathcal{D}_{t}}^{2}f_{i}+{O}({\delta
t}^{3})=2\beta\left(f_{i}^{\rm eq}-f_{i}\right),$ (100)

where we have only retained terms up to order two and
$\mathcal{D}_{t}=\partial_{t}+c_{i}\partial_{x}$. Then introducing
characteristic flow size $\mathcal{L}$ and velocity $\mathcal{U}$ the equation
is made non-dimensional as:

$\left(\frac{\delta
r}{\mathcal{L}}\right)\mathcal{D}_{t}^{\prime}f_{i}+\frac{1}{2}{\left(\frac{\delta
r}{\mathcal{L}}\right)}^{2}{\mathcal{D}_{t}^{\prime}}^{2}f_{i}=2\beta\left(f_{i}^{\rm
eq}-f_{i}\right),$ (101)

where primed variables denote non-dimensional form and

$\mathcal{D}_{t}^{\prime}=\frac{\mathcal{U}}{c}\left(\partial_{t}^{\prime}+c_{i}^{\prime}\partial_{x}^{\prime}\right),$
(102)

where $c=\delta r/\delta t$. Assuming acoustic, i.e.
$\frac{\mathcal{U}}{c}\sim 1$ and hydrodynamic, i.e. $\frac{\delta
r}{\mathcal{L}}\sim\varepsilon$, scaling and dropping the primes for the sake
of readability:

$\varepsilon\mathcal{D}_{t}f_{i}+\frac{1}{2}\varepsilon^{2}{\mathcal{D}_{t}}^{2}f_{i}+{O}(\varepsilon^{3})=2\beta\left(f_{i}^{\rm
eq}-f_{i}\right).$ (103)

Then introducing multi-scale expansions:

$\displaystyle f_{i}$ $\displaystyle=f_{i}^{(0)}+\varepsilon
f_{i}^{(1)}+\varepsilon^{2}f_{i}^{(2)}+O(\varepsilon^{3}),$ (104a)
$\displaystyle\partial_{t}$
$\displaystyle=\varepsilon\partial_{t}^{(1)}+\varepsilon^{2}\partial_{t}^{(2)}+O(\varepsilon^{3}),$
(104b) $\displaystyle\partial_{x}$ $\displaystyle=\varepsilon\partial_{x},$
(104c)

the following equations are recovered at scales $\varepsilon$ and
$\varepsilon^{2}$:

$\displaystyle\varepsilon$
$\displaystyle:\mathcal{D}_{t}^{(1)}f_{i}^{(0)}=-2\beta f_{i}^{(1)},$ (105a)
$\displaystyle\varepsilon^{2}$
$\displaystyle:\partial_{t}^{(2)}f_{i}^{(0)}+\mathcal{D}_{t}^{(1)}\left(1-\beta\right)f_{i}^{(1)}=-2\beta
f_{i}^{(2)},$ (105b)

with $f_{i}^{(0)}=f_{i}^{\rm eq}$. Taking the moments of the Chapman-Enskog-
expanded equation at order $\varepsilon$:

$\displaystyle\partial_{t}^{(1)}\rho+\partial_{x}\rho u_{x}$
$\displaystyle=0,$ (106a) $\displaystyle\partial_{t}^{(1)}\rho
u_{x}+\partial_{x}\rho
u_{x}^{2}+\partial_{x}\rho\varsigma^{2}+\partial_{x}P^{*}_{xx}$
$\displaystyle=0.$ (106b)

At order $\varepsilon^{2}$ the continuity equation is:

$\partial_{t}^{(2)}\rho=0.$ (107)

Summing up Eqs. 106a and 107 we recover the continuity equation as:

$\partial_{t}\rho+\partial_{x}\rho u_{x}+\mathcal{O}(\varepsilon^{3})=0,$
(108)

For the momentum equations we have:

$\partial_{t}^{(2)}\rho
u_{x}+\partial_{x}\left(1-\beta\right)\Pi_{xx}^{(1)}=0.$ (109)

where after some algebra the non-equilibrium stress tensor $\Pi_{xx}^{(1)}$
for any diagonal equilibrium pressure of the form:

$P_{\alpha\alpha}=\rho\varsigma^{2}+P^{*}_{\alpha\alpha},$ (110)

results in:

$\Pi_{xx}^{(1)}=-\frac{1}{2\beta}\left[2A\rho\varsigma^{2}\partial_{x}u_{x}+B\partial_{x}\rho\right],$
(111)

with

$A=\left(1-\frac{3}{2}\frac{u_{x}^{2}}{\varsigma^{2}}-\frac{3}{2\rho\varsigma^{2}}u_{x}\partial_{u_{x}}P^{*}_{xx}-\frac{{(\partial_{u_{x}}P^{*}_{xx})}^{2}}{2\rho^{2}\varsigma^{2}}-\frac{\partial_{\rho}P^{*}_{xx}}{2\varsigma^{2}}\right),$
(112)

and

$B=-3u_{x}\partial_{\rho}P^{*}_{xx}-\frac{\varsigma^{2}}{\rho}\partial_{u_{x}}P^{*}_{xx}-\frac{\partial_{u_{x}}P^{*}_{xx}\partial_{\rho}P^{*}_{xx}}{\rho}-u_{x}^{3}.$
(113)

##  Appendix B Multi-dimensional case

For $D>1$, using the same procedure as the 1-D case we have the following
system of equations at different orders in $\varepsilon$:

$\displaystyle\varepsilon$
$\displaystyle:\mathcal{D}_{t}^{(1)}f_{i}^{(0)}=-2\beta f_{i}^{(1)},$ (114a)
$\displaystyle\varepsilon^{2}$
$\displaystyle:\partial_{t}^{(2)}f_{i}^{(0)}+\mathcal{D}_{t}^{(1)}\left(1-\beta\right)f_{i}^{(1)}=-2\beta
f_{i}^{(2)},$ (114b)

where now $\mathcal{D}^{(1)}=\partial_{t}^{(1)}+\bm{c}_{i}\cdot\bm{\nabla}$,
which at order $\varepsilon$ lead to:

$\displaystyle\partial_{t}^{(1)}\rho+\bm{\nabla}\cdot\rho\bm{u}$
$\displaystyle=0,$ (115a)
$\displaystyle\partial_{t}^{(1)}\rho\bm{u}+\bm{\nabla}\cdot\rho\bm{u}\otimes\bm{u}+\bm{\nabla}\cdot\rho\varsigma^{2}\bm{I}+\bm{\nabla}\cdot\bm{P}^{*}$
$\displaystyle=0.$ (115b)

where the tensor $\bm{P}^{*}$ is diagonal with elements
$P^{*}_{\alpha\alpha}$. As for the 1-D case at order $\varepsilon^{2}$ the
continuity equation is:

$\partial_{t}^{(2)}\rho=0.$ (116)

For the momentum equations we have:

$\partial_{t}^{(2)}\rho\bm{u}+\bm{\nabla}\cdot\left(1-\beta\right)\Pi_{2}^{(1)}=0,$
(117)

with:

$\Pi_{2}^{(1)}=-\frac{1}{2\beta}\left[\partial_{t}^{(1)}\Pi_{2}^{(0)}+\bm{\nabla}\cdot\Pi_{3}^{(0)}\right].$
(118)

with:

$\Pi_{2}^{(1)}=-\frac{1}{2\beta}\left[\partial_{t}^{(1)}\Pi_{2}^{(0)}+\bm{\nabla}\cdot\Pi_{3}^{(0)}\right].$
(119)

where summation over $\gamma$ is implied. After some algebra one ends up with:

$\partial_{t}^{(1)}\Pi_{\alpha\beta}^{(0)}=-\partial_{\gamma}\rho
u_{\alpha}u_{\beta}u_{\gamma}\\\
-u_{\alpha}\partial_{\beta}(P_{\beta\beta}^{*}+\rho\varsigma^{2})-u_{\beta}\partial_{\alpha}(P_{\alpha\alpha}^{*}+\rho\varsigma^{2})\\\
+\delta_{\alpha\beta}\partial_{t}^{(1)}(P_{\alpha\alpha}^{*}+\rho\varsigma^{2}),$
(120)

while for the third-order equilibrium moments:

$\partial_{\gamma}\Pi^{(0)}_{\alpha\beta\gamma}=\partial_{\gamma}\rho
u_{\alpha}u_{\beta}u_{\gamma}+\left[\partial_{\gamma}u_{\gamma}(\rho\varsigma^{2}\delta_{\alpha\beta}+P^{*}_{\alpha\beta})\right]_{\rm
perm}\\\
-\delta_{\alpha\beta\gamma}(\left[\partial_{\gamma}u_{\gamma}P_{\alpha\beta}^{*}\right]_{\rm
perm}+\partial_{\gamma}\rho u_{\alpha}u_{\beta}u_{\gamma})\\\
-E_{\alpha\beta\gamma},$ (121)

where $E_{\alpha\beta\gamma}$ is a term representing errors specific to the
second order polynomial equilibrium. For both product form and entropic
equilibria $E_{\alpha\beta\gamma}=0$ while for the second order polynomial
equilibrium:

$E_{\alpha\beta\gamma}=(1-\delta_{\alpha\beta\gamma})\delta_{\alpha\beta}\left(\partial_{\gamma}\rho
u_{\alpha}u_{\beta}u_{\gamma}+\partial_{\gamma}u_{\gamma}P_{\alpha\beta}^{*}\right).$
(122)

Combining Eqs. (120) and (121):

$\Pi_{2}^{(1)}=-\frac{1}{2\beta}\left(\bm{S}^{*}+\bm{D}^{*}\right),$ (123)

with:

$S_{\alpha\beta}^{*}=(\rho\varsigma^{2}+P^{*}_{\alpha\alpha})\partial_{\alpha}u_{\beta}+(\rho\varsigma^{2}+P^{*}_{\beta\beta})\partial_{\beta}u_{\alpha}+E^{s}_{\alpha\beta},$
(124)

and:

$D_{\alpha\alpha}^{*}=-3\left(\rho
u_{\alpha}^{2}+P_{\alpha\alpha}^{*}+u_{\alpha}\partial_{u_{\alpha}}P_{\alpha\alpha}^{*}+\frac{{\left(\partial_{u_{\alpha}}P_{\alpha\alpha}^{*}\right)}^{2}}{3\rho}\right)\partial_{\alpha}u_{\alpha}\\\
+\left(P_{\alpha\alpha}^{*}-\rho\partial_{\rho}P_{\alpha\alpha}^{*}\right)\partial_{\gamma}u_{\gamma}\\\
-\left(u_{\alpha}^{3}+3u_{\alpha}\partial_{\rho}P_{\alpha\alpha}^{*}+\frac{\varsigma^{2}}{\rho}\partial_{u_{\alpha}}P_{\alpha\alpha}^{*}+\frac{1}{\rho}\partial_{\rho}P_{\alpha\alpha}^{*}\partial_{u_{\alpha}}P_{\alpha\alpha}^{*}\right)\partial_{\alpha}\rho\\\
+E^{d}_{\alpha\beta}.$ (125)

Here:

$E^{s}_{\alpha\beta}=-\left(P^{*}_{\alpha\alpha}+\rho
u_{\alpha}^{2}\right)\partial_{\alpha}u_{\beta}-\left(P^{*}_{\beta\beta}+\rho
u_{\beta}^{2}\right)\partial_{\beta}u_{\alpha}\\\
-u_{\alpha}^{2}u_{\beta}\partial_{\alpha}\rho-
u_{\beta}^{2}u_{\alpha}\partial_{\beta}\rho-2\rho
u_{\alpha}u_{\beta}(\partial_{\alpha}u_{\beta}+\partial_{\beta}u_{\alpha}),$
(126)

and

$E^{d}_{\alpha\alpha}=3\left(\rho
u_{\alpha}^{2}+P_{\alpha\alpha}^{*}+u_{\alpha}\partial_{u_{\alpha}}P_{\alpha\alpha}^{*}+\frac{{\left(\partial_{u_{\alpha}}P_{\alpha\alpha}^{*}\right)}^{2}}{3\rho}\right)\partial_{\alpha}u_{\alpha}\\\
-\left(P_{\alpha\alpha}^{*}-\rho\partial_{\rho}P_{\alpha\alpha}^{*}\right)\partial_{\gamma}u_{\gamma}\\\
+\left(u_{\alpha}^{3}+3u_{\alpha}\partial_{\rho}P_{\alpha\alpha}^{*}+\frac{\varsigma^{2}}{\rho}\partial_{u_{\alpha}}P_{\alpha\alpha}^{*}+\frac{1}{\rho}\partial_{\rho}P_{\alpha\alpha}^{*}\partial_{u_{\alpha}}P_{\alpha\alpha}^{*}\right)\partial_{\alpha}\rho\\\
+8\rho
u_{\alpha}^{2}\partial_{\alpha}u_{\alpha}+2u_{\alpha}^{3}\partial_{\alpha}\rho-\rho
u_{\alpha}^{2}\partial_{\gamma}u_{\gamma}-u_{\alpha}^{2}u_{\gamma}\partial_{\gamma}\rho-2u_{\alpha}u_{\gamma}\partial_{\gamma}u_{\alpha},$
(127)

where summation over $\gamma$ is implied in Eqs. (125), (126) and (127). Note
that in 1-D, the diagonal terms in Eqs. (124) and (125) reduce to Eqs. (112)
and (113).

##  Appendix C Jacobian of entropic equilibrium

Re-writing the entropic equilibrium as:

$f_{i}^{\rm eq}=w_{i}\rho\prod_{\alpha=x,y}\Lambda_{\alpha}\Gamma_{i\alpha},$
(128)

the Jacobian can be expressed as:

$\frac{\partial f_{i}^{\rm eq}}{\partial f_{j}}=f_{i}^{\rm
eq}\left(\frac{\partial\ln\rho}{\partial
f_{j}}+\sum_{\alpha=x,y}\frac{\partial\ln\Lambda_{\alpha}}{\partial
f_{j}}+\frac{\partial\ln\Gamma_{i\alpha}}{\partial f_{j}}\right),$ (129)

with:

$\displaystyle\Lambda_{\alpha}$ $\displaystyle=2-S_{\alpha},$ (130a)
$\displaystyle\Gamma_{i\alpha}$
$\displaystyle={\left(\frac{2u_{\alpha}+S_{\alpha}}{1-u_{\alpha}}\right)}^{c_{i\alpha}},$
(130b) $\displaystyle S_{\alpha}$
$\displaystyle=\sqrt{{(u_{\alpha}/\varsigma)}^{2}+1},$ (130c)

and :

$\displaystyle\frac{\partial S_{\alpha}}{\partial f_{j}}$
$\displaystyle=\frac{1}{\rho}\left[\frac{c_{j\alpha}u_{\alpha}/\varsigma^{2}+1}{S_{\alpha}}-S_{\alpha}\right],$
(131a) $\displaystyle\frac{\partial\Lambda_{\alpha}}{\partial f_{j}}$
$\displaystyle=-\frac{1}{\rho}\left[\frac{c_{j\alpha}u_{\alpha}/\varsigma^{2}+1}{S_{\alpha}}-S_{\alpha}\right],$
(131b) $\displaystyle\frac{\partial\Gamma_{i\alpha}}{\partial f_{j}}$
$\displaystyle=c_{i\alpha}\left[\frac{2c_{i\beta}+S_{\alpha}+\frac{\partial
S_{\alpha}}{\partial
f_{j}}}{\rho\left(1-u_{\alpha}\right)}-\frac{\left(2u_{\alpha}+S_{\alpha}\right)\left(1-c_{i\beta}\right)}{\rho{\left(1-u_{\alpha}\right)}^{2}}\right]$
$\displaystyle{\left(\frac{u_{\alpha}+S_{\alpha}}{1-u_{\alpha}}\right)}^{c_{i\alpha}-1}.$
(131c)