# Essentially entropic lattice Boltzmann model: Theory and simulations

Mohammad Atif  Jawaharlal Nehru Centre for Advanced Scientific Research,
Jakkur, Bangalore 560064, India  Praveen Kumar Kolluru  Jawaharlal Nehru
Centre for Advanced Scientific Research, Jakkur, Bangalore 560064, India
Santosh Ansumali<EMAIL_ADDRESS>Jawaharlal Nehru Centre for Advanced
Scientific Research, Jakkur, Bangalore 560064, India  SankhyaSutra Labs
Limited, Bangalore, India

###### Abstract

We present a detailed description of the essentially entropic lattice
Boltzmann model. The entropic lattice Boltzmann model guarantees unconditional
numerical stability by iteratively solving the nonlinear entropy evolution
equation. In this paper we explain the construction of closed-form analytic
solutions to this equation. We demonstrate that near equilibrium this exact
solution reduces to the standard lattice Boltzmann model. We consider a few
test cases to show that the exact solution does not exhibit any significant
deviation from the iterative solution. We also extend the analytical solution
for the ES-BGK model to remove the limitation on the Prandtl number for heat
transfer problems. The simplicity of the exact solution removes the
computational overhead and algorithmic complexity associated with the entropic
lattice Boltzmann models.

The lattice Boltzmann model (LBM) is an efficient kinetic formulation of the
nonlinear hydrodynamic phenomena on a lattice designed to capture the physics
of macroscopic flow (Frisch  _et al._ , 1986; Chen  _et al._ , 1992; Ansumali
_et al._ , 2003; Yudistiawan  _et al._ , 2010; Adhikari  _et al._ , 2005;
Mazloomi  _et al._ , 2015; Kolluru  _et al._ , 2020a). The Navier-Stokes
dynamics emerges as the hydrodynamic limit of this kinetic model which
performs simple microscale operations on the populations of fictitious
particles (Higuera  _et al._ , 1989; Qian  _et al._ , 1992; Benzi  _et al._ ,
1992). The discrete equilibrium in LBM is chosen such that the macroscopic
constraints are satisfied (McNamara and Zanetti, 1988; Qian  _et al._ , 1992;
Benzi  _et al._ , 1992). Historically, the top-down approach of choosing the
discrete equilibrium distribution from the macroscopic dynamics emerged as a
computationally attractive alternative to the Boolean particle dynamics of the
lattice gas model (Frisch  _et al._ , 1986; McNamara and Zanetti, 1988;
Higuera  _et al._ , 1989). However, this top-down approach lost a few
desirable features of the lattice gas such as the unconditional numerical
stability, the $H$ theorem and consequently the faithful representation of
microscopic Boltzmann dynamics (Karlin  _et al._ , 1999; Succi  _et al._ ,
2002). It was soon realized that the lack of a discrete time $H$ theorem
results in the growth of numerical instabilities (Boghosian  _et al._ , 2001;
Karlin  _et al._ , 1999; Succi  _et al._ , 2002).

The entropic lattice Boltzmann model (ELBM) emerged as an alternate
methodology to restore the $H$ theorem for discrete space-time evolution
(Karlin  _et al._ , 1998; Wagner, 1998; Karlin  _et al._ , 1999; Chen and
Teixeira, 2000; Boghosian  _et al._ , 2001; Succi  _et al._ , 2002; Ansumali
_et al._ , 2003; Boghosian  _et al._ , 2003). It was considered a paradigm
shift for computational fluid dynamics because the numerical stability of a
hydrodynamic solver was ensured by compliance with the thermodynamics at the
discrete time level (Succi  _et al._ , 2002). Currently, the ELBM is accepted
as a viable tool for simulation of turbulence, multiphase flows, as well as
microflows due to its unconditional numerical stability, and has shown
remarkable improvement over the traditional LBM (Ansumali  _et al._ , 2006;
Aidun and Clausen, 2010; Chikatamarla and Karlin, 2013; Mazloomi  _et al._ ,
2015; Atif  _et al._ , 2017). The additional step in ELBM, known as the
entropic involution step, involves a numerical search for the discrete path
length corresponding to jump to a mirror state on the isentropic surface.
Considerable efforts have been made to ensure the correctness and efficient
implementation of this step (Ansumali and Karlin, 2000, 2002a; Tosi  _et al._
, 2006; Chikatamarla  _et al._ , 2006; Brownlee  _et al._ , 2007; Gorban and
Packwood, 2012). However, there is scope for a better theoretical
understanding of the ELBM if one is able to obtain a closed form expression
for the discrete path length. For example:

  * •

The variable discrete path length could be understood as an adaptive implicit
modeling of the unresolved scales of the flow via the thermodynamic route, and
may provide a new insight into the subgrid modeling of turbulence.

  * •

It should enhance the efficiency of the ELBM by avoiding a numerical search
for the path length.

  * •

It will resolve the ambiguities in the implementation of ELBM. It should be
noted that for some rare events, the details of which are discussed in Sec.
II, the entropic involution step has no solution, and hence there is no unique
definition of the path length (Gorban and Packwood, 2012).

In Ref. (Atif  _et al._ , 2017), the authors reformulated the ELBM and
obtained a closed form analytical solution for the discrete path length
$\alpha$. This was achieved by relaxing the entropy equality condition used in
ELBM and replace it with the constraint that entropy must increase within a
discrete time step. The analytical form of $\alpha$ was found as the root of a
quadratic equation $-a\alpha^{2}+b\alpha-c$, where the coefficients $a,b,c$
are given in Eq. (46). The near equilibrium limit of this exact solution is
the standard LBGK value of $\alpha=2$. Its simplicity removes the
computational overhead and algorithmic complexity associated with ELBM. In
this paper, we discuss the theory of the entropic lattice Boltzmann model and
explain the construction of the closed form analytic solution for the discrete
path length in detail. We also demonstrate that the exact solution exhibits no
significant deviation from the iterative ELBM solution by considering a few
canonical setups. This paper is organized as follows: In Sec. I, we briefly
review the entropic lattice Boltzmann model. In Sec. II, we describe the
entropic involution step in its traditional form and derive its near-
equilibrium limit. In Sec. III, we explain the methodology to construct exact
solutions for the path length. In Sec. IV, we perform a detailed comparison of
the our solution with ELBM and BGK values of path length. In Sec. V we derive
the analytical solution to the path length for the ES-BGK model. Finally, in
Sec. VI we derive the expression for turbulent viscosity corresponding to the
exact solution of the path length.

##  I Entropic lattice Boltzmann model

In this section, we introduce the LBM and its entropic formulation in $D$
dimensions. In LBM one defines a set of discrete velocities ${\bf c}_{i}$,
$i=1,\cdots,N$ such that they form links of a space-filling lattice (Succi,
2001), and at every lattice node ${\bm{x}}$ and time $t$ a set of discrete
populations $f({\bm{c}}_{i},{\bm{x}},t)\equiv f_{i}$. Here, the set of
populations $f_{i}$ is understood as a vector
$\bm{f}=\\{f_{1},f_{2},\cdots,f_{N}\\}$ in the $N$ dimensional vector space,
where $N$ is the number of discrete populations. We define the bilinear action
between two functions of discrete velocities $\phi$ and $\psi$ as

$\left<\phi,\psi\right>=\sum_{i=1}^{N}\phi_{i}\psi_{i}.$ (1)

Analogous to continuous kinetic theory, the hydrodynamic variables such as the
mass density $\rho$, velocity $\mathbf{u}$, and the scaled temperature
$\theta$ are defined as

$\rho=\left<f,1\right>,\quad\rho{\bm{u}}=\left<f,{\bm{c}}\right>,\quad\rho
u^{2}+D\rho\theta=\left<f,{\bm{c}}^{2}\right>.$ (2)

Similarly, the $H$ function for hydrodynamics is taken in Boltzmann form as
(Karlin  _et al._ , 1999; Ansumali  _et al._ , 2003; Ansumali and Karlin,
2005)

$H[f]=\left<f,\log\frac{f}{w}-1\right>,$ (3)

with weights $w_{i}>0$. The population $\bm{f}({\bm{x}}+{\bm{c}}_{i}\Delta
t,t+\Delta t)$ after a time step $\Delta t$ starting from $\bm{f}({\bm{x}},t)$
is written as two step process:

  1. 1.

The discrete free-flight as

$\bm{f}({\bm{x}}+{\bm{c}}_{i}\Delta t,t+\Delta t)=\bm{f}^{*}({\bm{x}},t),$ (4)

which shifts the populations from one lattice node to another. Similar to the
free flight of molecules, this step preserves the entropy globally, i.e.,
$\sum_{\bm{x}}H[f({\bm{x}}+{\bm{c}}_{i}\Delta t,t+\Delta
t)]=\sum_{\bm{x}}H[f]$ (see Ref. (Wagner, 1998) for a detailed proof).

  2. 2.

The collisional relaxation towards the discrete equilibrium as

$\bm{f}^{*}({\bm{x}},t)=\bm{f}({\bm{x}},t)+\alpha\beta\left[\bm{f}^{\rm
eq}(\mathcal{M}^{\rm slow}({\bm{x}},t))-\bm{f}({\bm{x}},t)\right],$ (5)

typically modeled by a single relaxation model of Bhatnagar-Gross-Krook (BGK)
(Bhatnagar  _et al._ , 1954) with mean free time $\tau$. Here,
$\mathcal{M}^{\rm
slow}({\bm{x}},t)=\\{\rho({\bm{x}},t),\bm{u}({\bm{x}},t),\theta({\bm{x}},t)\\}$
are the collisional invariants ($\theta({\bm{x}},t)\notin\mathcal{M}^{\rm
slow}({\bm{x}},t)$ for isothermal LBM). For the standard LBGK, $\alpha=2$, and
the dimensionless discrete relaxation parameter $\beta={\Delta
t}/{(2\tau+\Delta t)}$ is bounded in the interval $0<\beta<1$. Notice that
$\beta=1$ implies $\tau=0$, and as the kinematic viscosity $\nu=\tau\theta$,
$\beta=1$ implies that there is no dissipation in the system. For a typical
LBM simulation the operating range is an over-relaxation regime of $\Delta
t/\tau\gg 1$ where $\beta\rightarrow 1$. In the standard LBM, this regime of
$\beta\rightarrow 1$ encounters numerical instability, which is resolved in
the ELBM by treating $\alpha$ as a variable which is evaluated at each point
and time step such that the $H$ theorem is satisfied. This is discussed in
detail in Sections II-III.

To recapitulate, the discrete free-flight that represents the convection
process leads to no dissipation, hence no entropy production (Wagner, 1998).
The collisional relaxation, however, has non-zero entropy production due to
relaxation of the populations towards the equilibrium but is entirely local in
position space.

Historically, the discrete isothermal equilibrium at a reference temperature
$\theta_{0}$ was chosen as (Qian  _et al._ , 1992)

$f_{i}^{\rm
eq}=w_{i}\rho\left[1+\frac{u_{\alpha}c_{\alpha}}{\theta_{0}}+\frac{u_{\alpha}u_{\beta}}{2\theta_{0}^{2}}\left(c_{\alpha}c_{\beta}-\theta_{0}\delta_{\alpha\beta}\right)\right],$
(6)

which was sufficient to recover the Navier-Stokes dynamics upto
$\mathcal{O}(u^{2})$, provided that the moments of the weights $w_{i}$ satisfy

$\left<w,1\right>=1,\,\left<w,c_{\alpha}c_{\beta}\right>=\theta_{0}\delta_{\alpha\beta},\left<w,c_{\alpha}c_{\beta}c_{\gamma}c_{\kappa}\right>=\theta_{0}^{2}\Delta_{\alpha\beta\gamma\kappa},$
(7)

where
$\Delta_{\alpha\beta\gamma\kappa}=\delta_{\alpha\beta}\delta_{\gamma\kappa}+\delta_{\alpha\gamma}\delta_{\beta\kappa}+\delta_{\alpha\kappa}\delta_{\beta\gamma}$.
However, this polynomial form of discrete equilibrium permits the populations
to attain negative values thus making the simulations numerically unstable
(Karlin  _et al._ , 1999; Succi  _et al._ , 2002). A method that resolves the
issue of nonpositive form of equilibrium distribution is to construct the
discrete equilibrium $\bm{f}^{\rm eq}$ as the minimizer of the convex $H$
function under the constraint that the mass density, the momentum density, and
the energy density (ignored for isothermal scenarios) are conserved (Karlin
_et al._ , 1999; Boghosian  _et al._ , 2001; Atif  _et al._ , 2018; Kolluru
_et al._ , 2020b). The discrete entropic equilibrium thus obtained is of the
form

$f_{i}^{\rm eq}=w_{i}\rho\exp\left(-\mu-\zeta_{\alpha}c_{i\alpha}-\gamma
c_{i}^{2}\right),$ (8)

where $\mu,\zeta_{\alpha},\gamma$ are the Lagrange multipliers. For the $D1Q3$
model, the discrete entropic isothermal equilibrium in the explicit form is

$f_{\pm 1}^{\rm
eq}=\frac{\rho}{6}\,\varUpsilon\left[\frac{2{u}_{\alpha}+\sqrt{1+3{u}_{\alpha}^{2}}}{1-{u}_{\alpha}}\right]^{\pm
1},\quad f_{0}^{\rm eq}=\frac{4\rho}{6}\,\varUpsilon,$ (9)

where $\varUpsilon=2-\sqrt{1+3{u}^{2}}$. For the higher-dimensional extensions
of $D1Q3$, i.e., $D2Q9,D3Q27,$ the generalized expression of the discrete
entropic isothermal equilibrium is (Ansumali  _et al._ , 2003)

$f_{i}^{\rm
eq}=w_{i}\rho\prod_{\alpha=1}^{D}\varUpsilon\left[\frac{2{u}_{\alpha}+\sqrt{1+3{u}_{\alpha}^{2}}}{1-{u}_{\alpha}}\right]^{c_{i\alpha}/\sqrt{3\theta_{0}}}.$
(10)

The above entropic equilibrium can be compared with Eq. (6) by performing a
series expansion around $u=0$. The expansion up to ${\cal O}(u^{3})$ is

$\displaystyle\begin{split}f_{i}^{\rm
eq}=w_{i}\rho\bigg{[}1+\frac{u_{\alpha}c_{\alpha}}{\theta_{0}}+\frac{u_{\alpha}u_{\beta}}{2\theta_{0}^{2}}\left(c_{\alpha}c_{\beta}-\theta_{0}\delta_{\alpha\beta}\right)\\\
+\frac{1}{6\theta_{0}^{3}}\left(u_{\alpha}u_{\beta}u_{\gamma}c_{\alpha}c_{\beta}c_{\gamma}-\theta_{0}u^{2}u_{\alpha}c_{\alpha}\right)\bigg{]},\end{split}$
(11)

which matches the historically employed equilibrium from Eq. (6) till ${\cal
O}(u^{2})$. The errors in the higher moments such as viscous stress and heat
flux is of ${\cal O}(u^{4})$ and ${\cal O}(u^{3})$ respectively (Ansumali,
2004). As for most higher-order models, the Lagrange multipliers cannot be
evaluated in explicit form and need to be found numerically. The series form
can be used as an alternative for simulations at low Mach numbers (${\rm Ma}$)
defined as ${\rm Ma}=u/c_{s}$, where $c_{s}$ is the sound speed.

##  II The entropic involution

The existence of the entropy function $H$ accompanied with the entropic
equilibrium derived in a variational fashion provides an opportunity for
creating a nonlinearly stable numerical method (Karlin  _et al._ , 1999; Succi
_et al._ , 2002; Boghosian  _et al._ , 2001). As the advection process [Eq.
(4)] does not lead to entropy production (Chen and Teixeira, 2000), a
nonlinearly stable LBM can be achieved by making the collisional relaxation to
equilibrium [Eq. (5)] adhere to the $H$ theorem, or in other words, by
ensuring that there is nonpositive entropy production during the collision
(Karlin  _et al._ , 1999).

The physical domain is discretized into grid points, at each of which we
define a set of $N$ populations $\bm{f}=\\{f_{0},f_{1},\cdots f_{N-1}\\}$.
Each point has an entropy level $H$ associated with it. For example, at a grid
point with set of populations $\bm{f}^{+}=\\{f^{+}_{0},f^{+}_{1},\cdots
f^{+}_{N-1}\\}$, from Eq. (3), $H[f^{+}]$ is a scalar. The equilibrium
$\bm{f}^{\rm eq}$ is the point with the least value of $H$, as, by
construction, it is the minimizer of the convex entropy function $H$ under the
relevant constraints.

The collision step given by Eq. (5) is understood in geometric terms as
follows: in an $N$ dimensional phase space, starting from the pre-collisional
state $\bm{f}$, one covers a distance (path length) $\alpha\beta$ in the
direction of $\bm{f}^{\rm eq}-\bm{f}$ to reach the post-collisional state
$\bm{f}^{*}$, i.e.,

$\bm{f}^{*}=\bm{f}+\alpha\beta[\bm{f}^{\rm eq}-\bm{f}].$ (12)

Here, for convenience we have dropped the position and time coordinates
$\bm{x},t$ as the collision step is local in position space and instantaneous.
We first consider the $D1Q2$ lattice as an example to visualize the phase
space and discuss the entropic collisional dynamics. This one dimensional
lattice has only two populations $f_{1},f_{-1}$ with discrete velocities
$+1,-1$ respectively (see Fig. 1). Due to the lack of enough degrees of
freedom, the $D1Q2$ lattice does not conserve momentum and hence cannot model
hydrodynamics. The mass density ($\rho=f_{1}+f_{-1}$) is a conserved moment,
and the momentum density ($\rho u=f_{1}-f_{-1}$) becomes a nonconserved
moment. These two constraints can be inverted to obtain the relations

$f_{1}=\frac{\rho+\rho u}{2},\quad f_{-1}=\frac{\rho-\rho u}{2}.$ (13)

Figure 2 represents the isoentropic contours in the vector space for the
$D1Q2$ lattice. The criterion of mass conservation $f_{1}+f_{-1}=\rho$
dictates that the collisional dynamics for $\rho=1$ is restricted on the
straight line in the figure. The equilibrium is given by

${\bm{f}}^{\rm eq}\equiv\\{f_{1}^{\rm eq},f_{-1}^{\rm
eq}\\}=\\{\rho/2,\rho/2\\}.$ (14)

It can be seen from the Fig. 2 (bottom) that near the equilibrium the
isoentropy contours are almost circular. This property of the $H$ function
$(H=f_{1}\log f_{1}+f_{2}\log f_{2}-\rho-\rho\log 2)$ is valid for the higher
dimensional lattices as well.

Another model we consider is D1Q3, which will be used later for illustrating
the concepts of entropic involution. For the $D1Q3$ lattice, the populations
are $\\{f_{-1},f_{0},f_{1}\\}$ with discrete velocities $\\{-1,0,+1\\}$
respectively. The mass conservation constraint requires that
$f_{-1}+f_{0}+f_{1}=\rho$, a plane on which the entire discrete dynamics is
constrained (see Fig. 4). The equilibrium for the $D1Q3$ lattice is given by
Eq. (9). The conserved moments are the mass density $\rho=f_{-1}+f_{0}+f_{1}$
and momentum density $\rho u=f_{1}-f_{-1}$, whereas the nonconserved moment is
the stress $\sigma_{xx}=f_{1}+f_{-1}-f^{\rm eq}_{1}-f^{\rm eq}_{-1}$. These
three constraints can be inverted to obtain the relations

$\displaystyle\begin{split}&\tilde{f}_{-1}\equiv\frac{f_{-1}}{\rho}=\frac{\tilde{f}^{\rm
eq}_{1}+\tilde{f}^{\rm eq}_{-1}+\tilde{\sigma}_{xx}-u}{2},\\\
&\tilde{f}_{0}\equiv\frac{f_{0}}{\rho}=1-\tilde{\sigma}_{xx}-\tilde{f}^{\rm
eq}_{1}-\tilde{f}^{\rm eq}_{-1},\\\
&\tilde{f}_{1}\equiv\frac{f_{1}}{\rho}=\frac{\tilde{f}^{\rm
eq}_{1}+\tilde{f}^{\rm eq}_{-1}+\tilde{\sigma}_{xx}+u}{2},\end{split}$ (15)

where $\tilde{\sigma}_{xx}=\sigma_{xx}/\rho,\tilde{f}^{\rm eq}_{i}=f^{\rm
eq}_{i}/\rho$.

Figure 1: Discrete velocities in a D1Q2 model. This one dimensional lattice
has only two populations $f_{1},f_{-1}$ with discrete velocities $+1,-1$
respectively, and cannot model hydrodynamics due to the lack of enough degrees
of freedom.

Figure 2: Isoentropy contours for a $D1Q2$ lattice. It can be seen from zoomed
figure (bottom) that near the equilibrium the isoentropy contours become
almost circular.

We now define a mirror state

$\bm{f}^{\rm mirror}=\bm{f}+\alpha(\bm{f}^{\rm eq}-\bm{f}),$ (16)

which is essentially $\bm{f}^{*}$ from Eq. (12) with $\beta=1$. Here, we
remind that $\beta=1$ is a zero dissipation state, therefore, the mirror state
$\bm{f}^{\rm mirror}$ lies at the same entropy as the initial state $\bm{f}$,
i.e.,

$H[\bm{f}^{\rm mirror}]=H[\bm{f}].$ (17)

The aim of the entropic involution step is to find the $\alpha$ corresponding
to the mirror state. Note that all the states $\bm{f},\bm{f}^{*},\bm{f}^{\rm
mirror}$ are at a higher entropy level than $\bm{f}^{\rm eq}$. Hence, starting
from $\bm{f}$ and moving in the direction of $\bm{f}^{\rm eq}-\bm{f}$, the
value of $H$ decreases till the equilibrium state, after which it begins to
rise. The maximum allowable path length that could be covered is $\alpha$,
after which $H$ increases beyond its pre-collisional state, and the $H$
theorem is violated. This is depicted in Fig. 3 for the D1Q2 lattice.

Figure 3: Entropic collisional dynamics for $D1Q2$ lattice. Note that the pre-
collisional state $\bm{f}$ and the mirror state ${\bm{f}}^{\rm mirror}$ are at
the same entropy level.

Figure 4: Top: The polytope of positivity for the $D1Q3$ lattice is a
triangular section of the plane inside which all the populations are positive,
and outside of which one or more populations become negative. Bottom:
Representation of a pre-collisional state $\bm{f}$ for which the mirror state
is not defined.

There exists an important structure in the distribution functions space – the
polytope of positivity (Gorban and Packwood, 2012). It is the region inside
which all the populations are positive but outside of which one or more
populations become negative. The shaded triangular region in Fig. 4 (top) is
the polytope of positivity for the $D1Q3$ lattice. The entropic involution
does not yield a solution when the isoentropic surfaces are partially outside
the polytope of positivity. This is due to the presence of the logarithm in
the entropy function which is undefined when one of the populations is
negative. Figure 4 (bottom) shows a pre-collisional state $\bm{f}$ for which
the mirror state lies outside the triangle, hence cannot be defined.

In LBGK, the path length is fixed to a constant value of $\alpha_{\rm
LBGK}=2$. The ELBM introduces the concept of the state dependent $\alpha$
(Karlin  _et al._ , 1999), evaluated numerically by solving the nonlinear
equation [Eq. (17)] (Ansumali and Karlin, 2002a; Tosi  _et al._ , 2006;
Chikatamarla  _et al._ , 2006). Once the path length $\alpha$ and therefore
the mirror state are known, the post-collisional state is found by the linear
contraction

$\bm{f}^{*}=\bm{f}^{\rm mirror}-\alpha(1-\beta)[\bm{f}^{\rm
eq}-\bm{f}]=\bm{f}+\alpha\beta[\bm{f}^{\rm eq}-\bm{f}].$ (18)

Since $0<\beta<1$, it is guaranteed that $H[\bm{f}^{*}]<H[\bm{f}^{\rm
mirror}]$. To summarize, the ELBM ensures adherence to the $H$ theorem in the
collision by first “over-relaxing” the populations to an equal entropy (zero
dissipation) mirror state followed by adding dissipation, thus, ensuring a
nonpositive entropy production (Karlin  _et al._ , 1999).

Next, we discuss the near equilibrium limit of the entropic involution. In a
well resolved simulation, the departure of populations from the equilibrium is
small and the entropic involution step yields the solution $\alpha=\alpha_{\rm
LBGK}=2$. To demonstrate this, we define the dimensionless departure from the
equilibrium as

$\displaystyle x_{i}=\frac{f_{i}^{\rm eq}}{f_{i}}-1.$ (19)

As the populations $f_{i},f_{i}^{\rm eq}$ are positive, $x_{i}\in(-1,\infty)$.
Here, the lower limit is due to the extreme case of $f_{i}^{\rm eq}\rightarrow
0$, whereas the upper limit is due to $f_{i}\rightarrow 0$. Further, we
introduce a decomposition of distributions $f_{i}$ in terms of the departure
from equilibrium as (Gorban  _et al._ , 1996)

$\Omega^{+}=\\{f_{i}:x_{i}\geq 0\\},\quad\Omega^{-}=\\{f_{i}:-1<x_{i}<0\\}.$
(20)

This asymmetry in the range of $x$ is crucial in the subsequent derivation of
the exact solution. With this decomposition, we also partition the bilinear
action into two partial contributions

$\left<f,\psi\right>_{\Omega^{\pm}}=\sum_{f_{i}\in\Omega^{\pm}}f_{i}\psi_{i}.$
(21)

The path length $\alpha$ is the root of the equation

$\displaystyle\Delta H\equiv H[\bm{f}^{\rm mirror}]-H[\bm{f}]=0,$ (22)

which is simplified to obtain (see Appendix A for a detailed derivation)

$\displaystyle\begin{split}H[\bm{f}^{\rm
mirror}]-H[\bm{f}]=\left<f,\left(1+\alpha x\right)\log{\left(1+\alpha
x\right)}\right>\\\ -\alpha\left<f,x\log(1+x)\right>.\end{split}$ (23)

In a well resolved simulation, the dimensionless departure of populations from
the equilibrium is small, i.e., $|x_{i}|\ll 1$. Therefore, expanding the above
equation about $x_{i}=0$ via Taylor series one obtains

$H[\bm{f}^{\rm
mirror}]-H[\bm{f}]=\alpha\left(\frac{\alpha}{2}-1\right)\left<f,x^{2}\right>+O(x^{3}).$
(24)

Thus, for small departure from the equilibrium, the non-trivial root of
$H[\bm{f}^{\rm mirror}]-H[\bm{f}]=0$ is $\alpha=2$. Hence, in the limit
$x_{i}\rightarrow 0$, the ELBM reduces to the LBGK.

We now derive the expanded form of Eq. (23) for the $D1Q2$ lattice. As stated
earlier, the $D1Q2$ lattice lacks the degrees of freedom to model
hydrodynamics, however, it is simple enough to show the analytical form of
$H[\bm{f}^{\rm mirror}]-H[\bm{f}]$. The Eq. (23) for the $D1Q2$ lattice can be
expanded to obtain

$\displaystyle\begin{split}&\Delta H\equiv H[\bm{f}^{\rm mirror}]-H[\bm{f}]\\\
&=f_{1}(1+\alpha x_{1})\log(1+\alpha x_{1})-\alpha f_{1}x_{1}\log(1+x_{1})\\\
&+f_{-1}(1+\alpha x_{-1})\log(1+\alpha x_{-1})-\alpha
f_{-1}x_{-1}\log(1+x_{-1}).\end{split}$ (25)

For this lattice, $f^{\rm eq}_{1}=f^{\rm eq}_{-1}=\rho/2$, therefore,
$x_{1}=\rho/(2f_{1})-1,x_{-1}=\rho/(2f_{-1})-1$, substituting which in the
above equation along with Eq. (13) yields

$\displaystyle\begin{split}&\frac{\Delta H}{\rho}=\left[\frac{1+u-\alpha
u}{2}\right]\log\left[\frac{1+u-\alpha u}{1+u}\right]\\\
&+\left[\frac{1-u+\alpha u}{2}\right]\log\left[\frac{1-u+\alpha
u}{1-u}\right]+\frac{\alpha u}{2}\log\left[\frac{1-u}{1+u}\right].\end{split}$
(26)

It is seen from the above equation that the solution of $\Delta H=0$ is
independent of $\rho$. It can also be verified that $\alpha=2$ is a nontrivial
solution (this is due to the symmetric nature of $D1Q2$ and is not the case
for $D1Q3$ and other higher dimensional lattices). Figure 5 shows that the
solution for $\Delta H=0$ remains $\alpha=2$ at all values of $u$.

Figure 5: The solution for $\Delta H=0$ remains $\alpha=2$ at all values of
$u$ for the $D1Q2$ lattice (this is not the case for $D1Q3$ and other higher
lattices).

Next, we derive the expanded form of Eq. (23) for the $D1Q3$ lattice. We
define $\tilde{x}_{i}$ as the $x_{i}$ for the $D1Q3$ model which are
calculated by substituting the equilibrium from Eq. (15) into Eq. (19) as

$\displaystyle\begin{split}&\tilde{x}_{-1}=\frac{2\tilde{f}^{\rm
eq}_{-1}}{\tilde{f}^{\rm eq}_{1}+\tilde{f}^{\rm
eq}_{-1}+\tilde{\sigma}_{xx}-u}-1,\\\ &\tilde{x}_{0}=\frac{\tilde{f}^{\rm
eq}_{0}}{1-\tilde{\sigma}_{xx}-\tilde{f}^{\rm eq}_{1}-\tilde{f}^{\rm
eq}_{-1}}-1,\\\ &\tilde{x}_{1}=\frac{2\tilde{f}^{\rm eq}_{1}}{\tilde{f}^{\rm
eq}_{1}+\tilde{f}^{\rm eq}_{-1}+\tilde{\sigma}_{xx}+u}-1.\end{split}$ (27)

The above $\tilde{x}_{i}$ are substituted in Eq. (23) to obtain the entropy
evolution for $D1Q3$ as

$\displaystyle\begin{split}&\frac{\Delta
H}{\rho}=\tilde{f}_{1}\left[(1+\alpha\tilde{x}_{1})\log(1+\alpha\tilde{x}_{1})-\alpha\tilde{x}_{1}\log(1+\tilde{x}_{1})\right]\\\
&+\tilde{f}_{-1}\left[(1+\alpha\tilde{x}_{-1})\log(1+\alpha\tilde{x}_{-1})-\alpha\tilde{x}_{-1}\log(1+\tilde{x}_{-1})\right]\\\
&+\tilde{f}_{0}\left[(1+\alpha\tilde{x}_{0})\log(1+\alpha\tilde{x}_{0})-\alpha\tilde{x}_{0}\log(1+\tilde{x}_{0})\right],\end{split}$
(28)

which is then solved using Newton-Raphson scheme for the path length $\alpha$.
This path length is dependent on $\tilde{\sigma}_{xx}$ and $u$ of the initial
state $\bm{f}$. Figure 6 plots the values of $\alpha$ for various
$u,\tilde{\sigma}_{xx}$. It can be seen that the region corresponding to the
LBGK value of 2, becomes thinner as $|\bm{u}|$ increases, and that the
deviation of $\alpha$ from the LBGK value becomes larger as
$|\bm{\tilde{}}{\sigma}_{xx}|$ increases. Figure 6 (bottom) plots the path
length as a function of $\bm{\tilde{}}{\sigma}_{xx}$ for various values of the
velocity $|\bm{u}|$. The shaded portion of the Fig. 6 (top) represents the
regions (typically with large moments) where the initial state is well defined
(lies within the polytope of positivity), whereas the mirror state lies
outside the polytope of positivity, thus, for such cases, the entropic
involution shows indeterminacy. It should be noted that these events are rare
and even if one encounters such cases it is known how to construct the path
length (Ansumali and Karlin, 2002a; Mazloomi M.  _et al._ , 2015).

Figure 6: The heat map of $\alpha$ corresponding to $\Delta H=0$ at various
values of $u,\tilde{\sigma}_{xx}$ for the $D1Q3$ lattice. The shaded region
represents the part of moment space where the mirror state lies outside the
polytope of positivity.

We now discuss the significance of over-relaxation in the entropic involution
step over the under-relaxation. A numerical scheme via the first order Euler
discretization of the Boltzmann BGK equation is possible. It reads as

$\displaystyle f(\bm{x}+\bm{c}\Delta t,t+\Delta t)$
$\displaystyle=f(\bm{x},t)+\frac{\Delta t}{\tau}\left[f^{\rm
eq}-f(\bm{x},t)\right]$ $\displaystyle=\left(1-\frac{\Delta
t}{\tau}\right)f(\bm{x},t)+\frac{\Delta t}{\tau}f^{\rm eq},$ (29)

and exhibits unconditional numerical stability if $\Delta t\ll\tau$. The $H$
theorem for this scheme is trivially satisfied as the post-collisional state
is a convex combination of the pre-collisional state and the equilibrium
state. This is called an under-relaxing scheme as the discrete dynamics never
crosses over the equilibrium state and corresponds to $\alpha<1$. However, for
many practical applications the relevant time scales are multiple orders of
magnitude greater than $\Delta t$. Therefore, for faster convergence it is
required to have numerical scheme which permits large time steps, i.e.,
$\Delta t\gg\tau$ is desirable (which correspond to $\alpha>1$). The over-
relaxation of the populations to a mirror state is thus an important feature
of the discrete dynamics as it allows one to achieve large time steps.

##  III Exact solution to the path length: Essentially entropic lattice
Boltzmann model

As discussed in the previous section, the discrete path length $\alpha$ is
available as the nontrivial root of Eq. (23). This equation is highly
nonlinear and is typically solved by a combination of bisection and Newton-
Raphson method (Ansumali and Karlin, 2000, 2002b). Considerable efforts have
been put in to ensure that the correct solution is obtained in an efficient
manner (Ansumali and Karlin, 2002a; Tosi  _et al._ , 2006; Chikatamarla  _et
al._ , 2006; Brownlee  _et al._ , 2007). In this section, we present an
alternate construction of ELBM where the discrete path length $\alpha$ is
known in explicit form without any indeterminacy. The key idea is to obtain
$\alpha$ by directly considering the natural criterion of monotonic decrease
of $H$ with time (Atif  _et al._ , 2017). This implies solving an inequality

$\Delta H\equiv H[\bm{f}^{*}]-H[\bm{f}]<0.$ (30)

The above inequality, by construction, accepts multiple solutions. For
example, when $\alpha\leq 1$ the inequality is trivially satisfied as the new
state is a convex combination of the old state and the equilibrium (Wagner,
1998). However, one is interested in an over-relaxed collision, where the new
state is no longer a convex combination of the old state and equilibrium. This
corresponds to the real solutions of Eq. (30) in the range
$1<\alpha<\alpha^{\rm max}$, where $\alpha^{\rm max}=-1/\left(\beta x_{i}^{\rm
min}\right)$ is maximum possible pass-length corresponding to an edge of the
polytope of positivity beyond which the populations become negative (Karlin
_et al._ , 1999). Among the multiple solutions of the inequality, we are
looking for the maximal path length $\alpha$ such that ${\Delta H\rightarrow
0}$. As is the case with ELBM, the solution should reduce to standard LBM
close to equilibrium ($\alpha=2$). Indeed, the present methodology is valid
for both discrete velocity models of LBM as well as the continuous in velocity
Boltzmann-BGK equations, where the summation in the inner products needs to be
replaced by appropriate integrals.

The general idea behind obtaining an analytical expression for the path length
$\alpha$ is as follows: we intend to split $\Delta H$ into two parts,

$\Delta H=H(\alpha)+H^{(B)},$ (31)

where $H^{(B)}$ is chosen such that it is nonpositive, and $H(\alpha)=0$ is an
easily solvable polynomial whose root is the path length $\alpha$. The
discrete-time $H$ theorem is satisfied as $H^{(B)}$ is nonpositive and
contributes to the entropy production, i.e.,

$\Delta H=H^{(B)}\leq 0.$ (32)

A word of caution is in order here. As stated earlier, the inequality $\Delta
H\leq 0$ by construction accepts multiple solutions. These solutions are not
identical but differ in two ways:

  1. 1.

Not all the solutions reduce to LBGK $(\alpha_{\rm LBGK}=2)$ in the limit of
$x_{i}\rightarrow 0$. Our interest is only in the solutions that reduce to the
standard LBM for $x_{i}\rightarrow 0$.

  2. 2.

The entropy production corresponding to each solution dictates its dissipative
nature, i.e., as the magnitude of $H^{(B)}$ increases the dynamics becomes
more and more dissipative. This is the reason why we are interested in the
solution such that ${\Delta H\rightarrow 0}.$ This point will be elucidated in
the forthcoming section, where we derive two expressions for $\alpha$, one of
which is more dissipative than the other.

Following the procedure detailed in Appendix A the Eq. (30) is rewritten as

$\displaystyle\Delta H=$
$\displaystyle\left<f,\left(1+\hat{x}\right)\log{\left(1+\hat{x}\right)}\right>-\alpha\beta\left<f,x\log(1+x)\right>,$
(33)

where $\hat{x}=\alpha\beta x$. Under the decomposition given by Eq. (20), the
above equation becomes

$\displaystyle\begin{split}\Delta
H&=\left<f,\left(1+\hat{x}\right)\log{\left(1+\hat{x}\right)}\right>_{\Omega^{-}}+\left<f,\left(1+\hat{x}\right)\log{\left(1+\hat{x}\right)}\right>_{\Omega^{+}}\\\
&-\alpha\beta\left<f,x\log(1+x)\right>_{\Omega^{-}}-\alpha\beta\left<f,x\log(1+x)\right>_{\Omega^{+}}.\end{split}$
(34)

We now derive two solutions to $\Delta H\leq 0$ by splitting Eq. (34) into a
polynomial and an entropy production term as in Eq. (31). These solutions
require bounds on the logarithm. The lower order solution is constructed by
exploiting the loose bounds, whereas the higher order solution is derived by
exploiting the sharper bounds (see Appendix B for details on the bounds of
logarithm). Both the solutions are shown to reduce to the LBGK value of 2 for
$x_{i}\rightarrow 0$.

###  III.1 Lower order solution

In this section, we find the path length by exploiting the loose bounds on the
logarithms [Eqs. (71),(74),(76)]. Upon adding and subtracting the term
$\left<f,{\cal A}_{1}+{\cal A}_{2}-{\cal A}_{3}\right>$ from Eq. (34), it is
written as

$\displaystyle\begin{split}\Delta
H=\left<f,\left(1+\hat{x}\right)\log{\left(1+\hat{x}\right)}-{\cal
A}_{1}\right>_{\Omega^{-}}+\left<f,{\cal A}_{1}\right>_{\Omega^{-}}\\\
+\left<f,{\left(1+\hat{x}\right)\log{\left(1+\hat{x}\right)}-{\cal
A}_{2}}\right>_{\Omega^{+}}+\left<f,{\cal A}_{2}\right>_{\Omega^{+}}\\\
-\alpha\beta\left<f,{x\log(1+x)-{\cal A}_{3}}\right>-\alpha\beta\left<f,{\cal
A}_{3}\right>,\end{split}$ (35)

where

$\displaystyle\begin{split}{\cal
A}_{1}=\hat{x}+\frac{\hat{x}^{2}}{2}-\frac{\hat{x}^{3}}{2},\,{\cal
A}_{2}=\hat{x}+\frac{\hat{x}^{2}}{2},\,{\cal
A}_{3}=\frac{2x^{2}}{2+x}.\end{split}$ (36)

Now, identifying that
$\left<f,x\right>_{\Omega^{+}}+\left<f,x\right>_{\Omega^{-}}=\left<f,x\right>=0$
due to conservation laws, Eq. (35) is written in a compact form as

$\displaystyle\Delta H$ $\displaystyle=\alpha\beta H_{1}(\alpha)+H_{1}^{(B)},$
(37)

where

$\displaystyle\begin{split}H_{1}^{(B)}=-{\left<f,G_{1}(\hat{x})\right>_{\Omega_{-}}}-{\left<f,G_{2}(\hat{x})\right>_{\Omega_{+}}}-{\alpha\beta\left<f,G_{3}(x)\right>}\\\
+{\alpha^{2}\beta(\beta-1)\left<f,\frac{x^{2}}{2}\right>}-{\alpha^{3}\beta(\beta^{2}-1)\left<f,\frac{x^{3}}{2}\right>_{\Omega^{-}}},\end{split}$
(38)

and

$H_{1}(\alpha)=-\alpha^{2}a_{1}+\alpha b_{1}-c_{1},$ (39)

with

$\quad
a_{1}=\left<f,\frac{x^{3}}{2}\right>_{\Omega^{-}},b_{1}=\left<f,\frac{x^{2}}{2}\right>,c_{1}=\left<f,\frac{2x^{2}}{2+x}\right>.$
(40)

It can be seen that $H_{1}(0)<0<H_{1}(2)$, therefore, a positive root of Eq.
(39) bounded in $(0,2)$ exists. As Eq. (39) is constructed by employing lower
order bounds on the logarithm, this root is called $\alpha_{\rm Lower}$,

$\alpha_{\rm
Lower}=\frac{-b_{1}+\sqrt{b_{1}^{2}-4a_{1}c_{1}}}{-2a_{1}}=\frac{2c_{1}}{b_{1}+\sqrt{b_{1}^{2}-4a_{1}c_{1}}}.$
(41)

To avoid numerical issues related to the precision loss while dealing with
small numbers, in the above expression we have multiplied the root with its
conjugate (Press  _et al._ , 1992).

Due to the nonnegative nature of the functions $G_{1},\,G_{2},\,G_{3}$ in
their respective domains [Eqs. (71), (74), (76)], and $\beta<1$, each term in
Eq. (38) is nonpositive, hence, $H_{1}^{(B)}\leq 0$. Therefore, from Eq. (37)
we see that the $H$ theorem is satisfied because $H_{1}(\alpha_{\rm
Lower})=0$, hence,

$\displaystyle\Delta H=H_{1}^{(B)}\leq 0.$ (42)

Upon expanding $\alpha_{\rm Lower}$ and ignoring higher order terms one
obtains

$\lim_{x_{i}\rightarrow 0}\alpha_{\rm
Lower}=2-\frac{\left<f,x^{3}\right>_{\Omega^{+}}}{\left<f,x^{2}\right>}+3\frac{\left<f,x^{3}\right>_{\Omega^{-}}}{\left<f,x^{2}\right>},$
(43)

which has the limiting value of $2$. Thus, for small departures from
equilibrium where $x_{i}\rightarrow 0$, the scheme reduces to the standard
LBM. It is also evident from Eq. (43) that $\alpha_{\rm Lower}<2$. This is
important as it is known that for ELBM the path length fluctuates around the
standard LBGK value of $\alpha=2$ (Karlin  _et al._ , 2015), a feature of ELBM
not mimicked by $\alpha_{\rm Lower}$. In the next section, we construct
another path length $\alpha_{\rm Higher}$ that fluctuates about the standard
LBGK value of $\alpha=2$.

###  III.2 Higher order solution

In this section, we derive the path length $\alpha$ by exploiting the sharper
bounds on the logarithms [Eqs. (72), (75), (77)]. Following the same
methodology as the previous section, we add and subtract terms from Eq. (67)
to obtain

$\displaystyle\begin{split}\Delta H&=H^{(B)}+\alpha\beta
H(\alpha),\end{split}$ (44)

where $H^{(B)}<0$ and

$\displaystyle H(\alpha)=-\alpha^{2}a+\alpha b-c.$ (45)

The coefficients $a,b,c$ are

$\displaystyle\begin{split}a&=\beta^{2}\left<f,\frac{x^{3}}{6}-\frac{h\beta
x^{4}}{12}+\frac{h^{2}\beta^{2}x^{5}}{20}-\frac{h^{3}\beta^{3}x^{6}}{5}\right>_{\Omega^{-}},c=\left<f,\frac{60x^{2}+60x^{3}+11x^{4}}{60+90x+36x^{2}+3x^{3}}\right>,\\\
b&=\bigg{<}f,\frac{x^{2}}{2}\bigg{>}-\bigg{<}f,\frac{2\alpha_{\rm
Lower}\beta^{2}x^{3}}{15}\bigg{(}\frac{2}{4+\alpha_{\rm
Lower}x}+\frac{1}{4+2\alpha_{\rm Lower}x}+\frac{2}{4+3\alpha_{\rm
Lower}x}\bigg{)}\bigg{>}_{\Omega^{+}},\\\ \end{split}$ (46)

The parameter $h$ in the above equation serves as an upper bound on the path
length and is found as the positive root of the quadratic equation

$\displaystyle\begin{split}H_{2}(\alpha)&=-\alpha^{2}a_{2}+\alpha
b-c,\end{split}$ (47)

where

$\displaystyle\begin{split}a_{2}=\beta^{2}\left<f,\frac{x^{3}}{6}\right>_{\Omega^{-}}.\end{split}$
(48)

Equation (45) has a positive root $\alpha_{\rm Higher}$ [as
$H(0)<0<H(\infty)$] which is the desired path length. It has the limit

$\lim_{x_{i}\rightarrow 0}\alpha_{\rm
Higher}=2+\left(\frac{4\beta^{2}}{3}-1\right)\frac{\left<f,x^{3}\right>}{\left<f,x^{2}\right>}.$
(49)

Unlike $\alpha_{\rm Lower}$, which was always less than 2, no such comment can
be made about $\alpha_{\rm Higher}$. Thus, $\alpha_{\rm Higher}$ mimics an
important feature of the ELBM where the path length fluctuates about the BGK
value of 2. A detailed derivation of $\alpha_{\rm Higher}$ is provided in
Appendix C. The details regarding the implementation of this exact solution
for the path length are given in Appendix D.

##  IV Comparison with ELBM and BGK

In this section, we compare the analytical solutions for the path length
($\alpha_{\rm Lower}$, $\alpha_{\rm Higher}$) with the BGK ($\alpha_{\rm
LBGK}=2$) and the iterative ELBM solution ($\alpha_{\rm ELBM}$). To this end,
we consider three canonical setups: the one-dimensional Sod shock tube, the
doubly periodic shear layer, and the lid-driven cavity. It is illustrated from
these examples that $\alpha_{\rm Lower}$ is more dissipative than $\alpha_{\rm
Higher}$ and hence is not the ideal choice for hydrodynamics. Nevertheless, it
is useful for the construction of $\alpha_{\rm Higher}$ as demonstrated in the
previous section. It is also demonstrated that there is an insignificant
difference between the path lengths $\alpha_{\rm Higher}$ and $\alpha_{\rm
ELBM}$.

###  IV.1 Sod shock tube

Figure 7: Density (left), velocity (middle) and entropy (right) plots from
$\alpha_{\rm LBGK}$, $\alpha_{\rm Lower}$, $\alpha_{\rm Higher}$, and
$\alpha_{\rm ELBM}$ at time $t=500$ for viscosity $\nu=1.0\times 10^{-5}$.

To compare the behaviour of $\alpha_{\rm Lower},\alpha_{\rm Higher}$ with
$\alpha_{\rm ELBM}$ and $\alpha_{\rm LBGK}$, we first simulate the one-
dimensional shock tube using the $D1Q3$ lattice. In this setup, a domain with
800 grid points is initialized with a step function for density as
$\rho\,(x\leq 400)=1.5$ and $\rho\,(x>400)=0.75$. The presence of a sharp
discontinuity in the initial condition at the center of the domain generates a
moving compressive shock front in the low-density region and a rarefaction
front in the high-density region. These two fronts give rise to a contact
region of uniform pressure and velocity in the center of the tube (Laney,
1998). The density, velocity, and entropy profiles shown in Figure 7
illustrate that the numerical oscillations are sharply reduced in the case of
$\alpha_{\rm Lower}$, thus pointing to its dissipative nature. It can also be
seen that the oscillations are prominent for $\alpha_{\rm LBGK}$ and that both
$\alpha_{\rm Higher}$ and $\alpha_{\rm ELBM}$ restore the $H$ theorem without
altering the fields.

Figure 8: Comparison between $\alpha_{\rm Lower}$ and $\alpha_{\rm ELBM}$ for
the Sod shock tube. Top: Snapshot of the path length. Bottom: Ratio of the
turbulent viscosity correction to the kinematic viscosity.

Figure 9: Comparison between $\alpha_{\rm Higher}$ and $\alpha_{\rm ELBM}$ for
the Sod shock tube. Top: Snapshot of the path length. Bottom: Ratio of the
turbulent viscosity correction to the kinematic viscosity.

Figure 8 (top) compares $\alpha_{\rm Lower}$ and $\alpha_{\rm ELBM}$. It is
evident that the path lengths show departure from $\alpha=2$ (BGK value) only
in the narrow regions of the compressive and the rarefaction fronts. It can
also be seen that the value of $\alpha_{\rm Lower}$ is always smaller than 2,
while that of $\alpha_{\rm ELBM}$ fluctuates about 2. Figure 8 (bottom) plots
the ratio of turbulent viscosity correction to kinematic viscosity
$\nu_{T}/\nu_{0}$ (more details in Sec. VI). From the figure, it is evident
that at the location of the shock front the $\alpha_{\rm Lower}$ is more than
twice the kinematic viscosity, while $\alpha_{\rm ELBM}$ is only $\sim 47\%$.
Similarly, Figure 9 (top) compares the path length from $\alpha_{\rm Higher}$
and it is seen that for this setup $\alpha_{\rm Higher}$ exhibits smaller
fluctuations than the $\alpha_{\rm ELBM}$. Figure 9 (bottom) shows that the
turbulent viscosity correction for $\alpha_{\rm ELBM}$ is $\sim 47\%$, whereas
for $\alpha_{\rm Higher}$ it is $\sim 42\%$. Hence, it can be concluded that
$\alpha_{\rm Higher}$ imposes the $H$ theorem (thus guaranteeing unconditional
numerical stability) with the least turbulent viscosity correction.

###  IV.2 Doubly periodic shear layer

In this section, we compare the behaviour of $\alpha_{\rm Lower},\alpha_{\rm
Higher}$ with $\alpha_{\rm LBGK}$ by considering the setup of doubly periodic
shear layer (Minion and Brown, 1997). The initial velocity field comprises of
two shear layers given by

$\displaystyle u_{x}(y)$
$\displaystyle=\begin{cases}U_{0}\tanh[(4y-1)/w],\qquad y\leq 1/2\\\
U_{0}\tanh[(3-4y)/w],\qquad y>1/2\end{cases}$ (50) $\displaystyle u_{y}(x)$
$\displaystyle=U_{0}\delta\sin[2\pi(x+1/4)],$ (51)

where $w=\delta=0.05,U_{0}=0.04$ and $x,y$ are nondimensionalized coordinates.
The viscosity is calculated from the Reynolds number which for the present
case is fixed at $3\times 10^{4}.$ It is known that at poor grid resolutions
for this setup, the numerical disturbances may lead to formation of spurious
vortices in the braids (Minion and Brown, 1997; Coreixas  _et al._ , 2017).

Figure 10 depicts the isovorticity contours for $\alpha_{\rm
Lower},\alpha_{\rm Higher}$ on a $256\times 256$ grid and for $\alpha_{\rm
LBGK}$ on $1024\times 1024$ grid obtained after one convection time. A
qualitative comparison of the three plots reveal that the vortex structure is
smudged for $\alpha_{\rm Lower}$, while the vortex structure of $\alpha_{\rm
Higher}$ on a $256\times 256$ grid is the same as that of BGK at $1024\times
1024$ grid. In Fig. 11 we show the magnitude of the path lengths $\alpha_{\rm
Lower},\alpha_{\rm Higher}$, from where it evident that while $\alpha_{\rm
Lower}$ always remains smaller than 2, $\alpha_{\rm Higher}$ fluctuates about
2, thus corroborating the dissipative nature of $\alpha_{\rm Lower}$. Finally,
a quantitative analysis of the flow is performed by measuring the change in
global enstrophy ($\Delta\Omega=\bar{\Omega}_{t}/\bar{\Omega}_{0}\times 100$),
where $\bar{\Omega}_{t}$ is the enstrophy at time $t$, and $\bar{\Omega}_{0}$
is the initial global enstrophy (defined as the square of the vorticity).
Figure 12 plots the time evolution of $\Delta\Omega$. It is evident that
$\alpha_{\rm Higher}$ on a $128\times 128$ grid behaves the same as the BGK on
a much larger $1024\times 1024$ grid, whereas $\alpha_{\rm Lower}$ exhibits
dissipation that manifests in the form of reduced enstrophy.

Figure 10: Nondimensional iso-vorticity contours for $\alpha_{\rm Lower}$
(left), $\alpha_{\rm Higher}$ (center) at grid size $256\times 256$ and for
BGK at $1024\times 1024$ (right) after one convection time.

Figure 11: Path length from $\alpha_{\rm Lower}$ (left) and $\alpha_{\rm
Higher}$ (right) after one convection time on a grid of size $256\times 256$.
Figure 12: Change in the global enstrophy $\Delta\Omega$ vs time for various
square grids. Here, $t^{*}$ is the nondimensional convection time.

###  IV.3 Lid-driven cavity

In this section, we consider the lid-driven cavity at a Reynolds number ($\rm
Re$) of 5000 where the motion of the top wall drives the flow in a 2D cavity.
We use the standard $D2Q9$ lattice and diffuse boundary condition Ansumali and
Karlin (2002c). For this setup, the LBGK ($\alpha=2$) is numerically unstable
at smaller grid sizes of $64\times 64,\,96\times 96$, and $128\times 128$,
however, it is stable at a larger grid of size $256\times 256$. The entropic
formulations $\alpha_{\rm Lower},\alpha_{\rm Higher},\alpha_{\rm ELBM}$ are
stable at all grid sizes.

Figure 13: Iso-vorticity contours for the lid-driven cavity at Reynolds number
of 5000 for various grid sizes: $64\times 64$ (left), $96\times 96$ (center),
$128\times 128$ (right).

Figure 13 depicts the iso-vorticity contours for various grid sizes obtained
using $\alpha_{\rm Higher}$. It is seen that even extremely under-resolved
grids remain numerically stable. However, at coarse resolutions like $64\times
64$ and $96\times 96$ the finer structures are distorted, which take the
expected form at a slightly higher grid size of $128\times 128$. It should be
repeated here that at grid size of $128\times 128$ the LBGK ($\alpha=2$) is
numerically unstable. In Fig. 14, we plot the velocities along vertical and
horizontal centerlines and observe a good match with Ghia  _et al._ (1982).

Figure 14: Velocity profiles for the lid-driven cavity at Reynolds number of
5000 and Mach number 0.05 for various grid sizes. Top: nondimensionalized
x-velocity along the vertical centerline. Bottom: nondimensionalized
y-velocity along the horizontal centerline.

Next, we establish that there is no appreciable difference between the path
lengths $\alpha_{\rm Higher}$ and $\alpha_{\rm ELBM}$. To this effect, we
compare the instantaneous value of $\alpha_{\rm Higher}$ and $\alpha_{\rm
ELBM}$ for three different grid resolutions. First, the simulation is
performed using $\alpha_{\rm Higher}$ for 100 convection times. On the
populations thus obtained, we evaluate $\alpha_{\rm Higher}$ and $\alpha_{\rm
ELBM}$ for the entire grid. The $L_{1},\,L_{2},\,L_{\infty}$ error norms of
$||\alpha_{\rm Higher}-\alpha_{\rm ELBM}||$ are tabulated in Table 1, whereas
the distribution of path lengths are given in Fig. 15. It is evident that
$\alpha_{\rm Higher}$ and $\alpha_{\rm ELBM}$ show insignificant deviation at
all grid sizes. From Fig. 15 and Table 2, it can also be seen that as the grid
size increases the distribution of the path lengths becomes narrower as the
region around the LBGK value of $\alpha=2$ where $90\%$ of the points lie
(inside solid vertical lines) becomes smaller.

Figure 15: Distribution of $\alpha_{\rm Higher}$ and $\alpha_{\rm ELBM}$ for lid-driven cavity at Reynolds number 5000 and Mach number 0.05. Grid sizes are $64\times 64$ (top), $96\times 96$ (middle), $128\times 128$ (bottom). The difference between the distribution of $\alpha_{\rm Higher}$ and $\alpha_{\rm ELBM}$ is seen to be insignificant. The solid black lines denote the region inside which $90\%$ of the points lie. The locations of the solid lines are tabulated in Table 2. | $64\times 64$ | $96\times 96$ | $128\times 128$
---|---|---|---
$L_{1}$ | $2.55\times 10^{-5}$ | $1.37\times 10^{-5}$ | $8.26\times 10^{-6}$
$L_{2}$ | $1.67\times 10^{-4}$ | $9.89\times 10^{-5}$ | $5.17\times 10^{-5}$
$L_{\infty}$ | $6.07\times 10^{-3}$ | $5.24\times 10^{-3}$ | $3.71\times 10^{-3}$
Table 1: Error norms for $||\alpha_{\rm Higher}-\alpha_{\rm ELBM}||$. | $64\times 64$ | $96\times 96$ | $128\times 128$
---|---|---|---
$\alpha_{\rm Higher}$ | $2\pm 1.79\times 10^{-3}$ | $2\pm 8.0\times 10^{-4}$ | $2\pm 3.9\times 10^{-4}$
$\alpha_{\rm ELBM}$ | $2\pm 1.77\times 10^{-3}$ | $2\pm 7.3\times 10^{-4}$ | $2\pm 2.9\times 10^{-4}$
Table 2: Region around the LBGK value of $\alpha=2$ where $90\%$ of the points
lie. It is seen that as the grid size increases the region becomes narrower.

We also briefly investigate the idea that the path length $\alpha_{\rm Lower}$
could be utilized as a good initial guess value for the iterative ELBM solver.
Typically, the iterative root solver converges in 4-5 iterations, however, it
is stipulated that the converged result should be obtained in a single
iteration when using $\alpha_{\rm Lower}$ as the initial guess value. We call
this first iterate $\alpha_{\rm Iterate1}$ and compare it with $\alpha_{\rm
ELBM}$. The $L_{1},L_{2},L_{\infty}$ error norms of $||\alpha_{\rm
Iterate1}-\alpha_{\rm ELBM}||$ are tabulated in Table 3 from where it can be
concluded that the difference is insignificant for all three grid sizes.

| $64\times 64$ | $96\times 96$ | $128\times 128$
---|---|---|---
$L_{1}$ | $2.27\times 10^{-7}$ | $7.48\times 10^{-8}$ | $2.56\times 10^{-8}$
$L_{2}$ | $3.02\times 10^{-6}$ | $1.37\times 10^{-6}$ | $8.26\times 10^{-7}$
$L_{\infty}$ | $1.10\times 10^{-4}$ | $9.00\times 10^{-5}$ | $7.00\times 10^{-5}$
Table 3: Error norms for $||\alpha_{\rm Iterate1}-\alpha_{\rm ELBM}||$.

##  V Exact solution to the entropic lattice ES–BGK model

The ES–BGK model proposed by Holway Jr (1965) overcomes the restriction on the
Prandtl number ($\rm Pr$) in BGK collision models without compromising the
conceptual simplicity. This model employs a quasi-equilibrium state $f^{\rm
QE}$ instead Maxwellian in the collision term. The quasi-equilibrium state is
an anisotropic Gaussian distribution that reduces to a Maxwellian at the
equilibrium. The continuous $H$ theorem for this model was proved by Andries
_et al._ (2000). In this section, we extend the discrete $H$ theorem to the
lattice ES–BGK model and derive the exact solution for the path length.

###  V.1 Lattice ES–BGK model

The collision term for the lattice ES–BGK collision model reads as (Meng  _et
al._ , 2013)

$\bm{f}^{*}({\bm{x}},t)=\bm{f}({\bm{x}},t)+\alpha\beta[\tilde{\bm{f}}^{\rm
QE}-\bm{f}({\bm{x}},t)],$ (52)

where $\beta=\Delta t/(2\tau_{1}+\Delta t),$ and the viscosity $\nu$ is
related to the relaxation time $\tau_{1}$ by $\nu=\tau_{1}\theta\,{\rm Pr}$
(Kolluru  _et al._ , 2022). In Eq. (52), $\alpha$ is the path length which is
equal to 2 in the standard case, and is found by solving Eq. (30) for the
entropic lattice ES–BGK model. The discrete quasi-equilibrium distribution
$\tilde{f}^{\rm QE}$ is found as the minimizer of the the discrete $H$
function under the constraints of mass and momentum being conserved with the
pressure tensor given by

$\displaystyle P_{\alpha\beta}\left(\tilde{\bm{f}}^{\rm QE}\right)$
$\displaystyle=\rho
u_{\alpha}u_{\beta}+\rho\theta\delta_{\alpha\beta}+\frac{1-1/{\rm
Pr}}{1+\Delta t/(2\tau_{1}{\rm Pr})}\sigma_{\alpha\beta}\left(\bm{f}\right).$
(53)

Solving the minimization problem, one obtains

$\displaystyle\tilde{f}_{i}^{\rm
QE}=\rho\exp\left(-\mu-{\zeta_{\kappa}c_{i\kappa}}-\gamma_{\alpha\beta}\sigma_{\alpha\beta}\right),$
(54)

where $\mu,\zeta_{\kappa},\gamma_{\alpha\beta}$ are the Lagrange multipliers
associated with the mass, momentum, and pressure tensor respectively. The
Lagrange multipliers are calculated by performing a perturbation expansion
around the equilibrium state as in Ref. Ansumali  _et al._ (2007).

###  V.2 Exact solution for the path-length

Following the procedure as detailed in Appendix A, the Eq. (30) for the
lattice ES–BGK model is rewritten as

$\displaystyle\begin{split}\Delta
H=\left<f,(1+\hat{z})\ln(1+\hat{z})\right>-\alpha\beta\left<f,z\ln(1+z)\right>\\\
+\frac{\alpha\beta}{{\rm Pr}}\frac{1+\delta t/(2\tau_{1})}{1+\delta
t/(2\tau_{1}{\rm Pr})}\gamma_{\alpha\beta}\sigma_{\alpha\beta}(f),\end{split}$
(55)

where $\hat{z}=\alpha\beta z,z=\tilde{f}^{\rm QE}/f-1$. The Lagrange
multipliers are evaluated numerically, however, using a series expansion it
can be shown that the last term in the above equation can be approximated as

$\gamma_{\alpha\beta}=-m\frac{\sigma_{\alpha\beta}}{\rho\theta^{2}},\quad
m=\frac{1}{2}\,{\rm for}\,\alpha=\beta,\quad m=1\,{\rm otherwise}.$ (56)

It is seen that the last term is negative definite hence it contributes only
to the entropy production. Thus, the analytical expression for the path length
remains the same with equivalent features as Sec. III.

###  V.3 Rayleigh-Bénard convection

Rayleigh-Bénard convection is a well-studied model of natural convection and
is considered a classical benchmark for thermal models (Shan, 1997). The
domain consists viscous fluid confined between two thermally well-conducting
parallel plates. The plates are kept at a distance $L$ with the bottom plate
maintained at higher temperature $\theta_{\rm bottom}$ and the top plate is
kept at a lower temperature $\theta_{\rm top}$. The flow is induced by the
unstable density gradients in the presence of a gravitational field (Atif  _et
al._ , 2018). The dynamics of the Rayleigh-Bènard convection is characterized
by two non-dimensional numbers: the Rayleigh number and the Prandtl number.
The Prandtl number is a property of the fluid (${\rm Pr}=\nu/\alpha_{T}$)
whereas the Rayleigh number (${\rm Ra}$) is defined as

${\rm Ra}=\frac{{g}\hat{\beta}\Delta\theta L^{3}}{\nu\alpha_{T}},$ (57)

where ${g}$ is the gravity, $\hat{\beta}=-1/\rho(\partial\rho/\partial T)_{P}$
is the thermal expansion coefficient, $\Delta\theta=\theta_{\rm
bottom}-\theta_{\rm top}$ is the temperature difference between the two walls,
$\nu$ is the kinematic viscosity, and $\alpha_{T}$ is the thermal diffusivity.

In this section, we simulate the turbulent Rayleigh-Bénard convection at ${\rm
Ra}=1.0\times 10^{7}$ and ${\rm Pr}=0.71$ on a grid of size $2N\times 2N\times
N$ with $N=112$ and $N=224$. The exact solution for the path length as derived
in the preceding section is used with Eq. (52) as the collision model. The
numerical simulations are performed using the 67 velocity crystallographic
lattice Atif  _et al._ (2018) with $\theta_{\rm bottom}=1.02\theta_{0}$ and
$\theta_{\rm top}=0.98\theta_{0}$. Constant temperature boundary conditions at
the top and the bottom walls were imposed and periodic boundary conditions
were applied in the horizontal directions. We calculate the Nusselt number and
time-averaged horizontal mean of nondimensional temperature
$T=(\theta-\theta_{\rm top})/\Delta\theta$. The calculated Nusselt number is
$13.4$ with $N=112$ and $15.3$ with $N=224$, whereas that reported by the
direct numerical simulation (DNS) of Ref. Togni  _et al._ (2015) is 15.59. In
Fig. 17 we compare the time-averaged mean horizontal temperature with the DNS
data and observe a good match. It can be seen that as expected the temperature
rises rapidly close to the wall and obtains a uniform profile in the bulk.
Hence, it can be concluded that the exact solution to the path length extends
the unconditional numerical stability to non-unity Prandtl number heat
transfer simulations too.

Figure 16: Iso-temperature contours (top) for Rayleigh-Bènard convection at
nondimensional temperatures 0.3, 0.7. The mid and bottom figures visualize the
temperature field at horizontal slices close to the two walls.

Figure 17: Time-averaged mean horizontal temperatures for Rayleigh-Bènard
convection at ${\rm Ra}=1\times 10^{7}$ and ${\rm Pr}=0.71$ compared with the
DNS profile from Ref.Togni  _et al._ (2015). Here, the nondimensional vertical
coordinate $z^{*}$ is $z\,{\rm Nu}/L$

.

##  VI Entropic route to modeling the subgrid viscosity

The entropic LBM has been interpreted as an implicit subgrid model of
turbulence Karlin  _et al._ (2003). The modification to the path length
$\alpha$ due to the compliance with the $H$ theorem can be understood as a
turbulent viscosity correction at the macroscopic scale. Several studies have
analyzed the form of the viscosity correction and found similarities to the
Smagorinsky’s model for viscosity correction Malaspinas  _et al._ (2008);
Buzzicotti and Tauzin (2021). In this section, we derive the subgrid model
corresponding to the exact path length. From the Chapman-Enskog expansion the
effective kinematic viscosity $\nu$ due to the entropic collision term is
found as

$\nu=\theta\tau=\theta\Delta t\left(\frac{1}{\alpha\beta}-\frac{1}{2}\right).$
(58)

The viscosity correction $\nu_{T}$ is defined as $\nu_{T}=\nu-\nu_{0}$, where
$\nu_{0}$ is the viscosity corresponding to the BGK path length $\alpha=2$,
and is obtained as

$\nu_{T}=\frac{\theta\Delta t}{2\alpha\beta}\left(2-{\alpha}\right).$ (59)

It is seen from the above expression that the path length $\alpha$ dictates
whether the viscosity correction is positive or negative. A path length
smaller than 2 implies an increment in the viscosity which in turn smoothens
the gradients, whereas, a path length larger than 2 corresponds to reduction
in the viscosity which sharpens the gradients (Karlin  _et al._ , 2015). Thus,
the entropic LBM permits backscatter of energy from the subgrid scales to the
resolved scales too.

We now evaluate the viscosity correction in terms of the macroscopic moments.
For this purpose we consider the path length from Eq. (49) in assuming small
departure from equilibrium, i.e., $x_{i}\rightarrow 0,\,\Delta t\gg\tau$, and
interpret $\left<\cdot\right>$ as a continuous integral. Substituting Eq. (49)
in Eq. (59) the turbulent viscosity correction $\nu_{T}$ is found as

$\nu_{T}=\frac{\theta\Delta t}{2}\frac{\Theta}{6+\Theta},\,{\rm
where}\,\quad\Theta=\frac{\left<f,x^{3}\right>}{\left<f,x^{2}\right>}.$ (60)

From Grad’s $13$ moment representation one can write the approximation
$f=f^{\rm MB}\left(1+\Omega\right)$, where

$\Omega=\frac{\sigma_{ij}\xi_{i}\xi_{j}}{2p\theta}-\frac{q_{k}\xi_{k}}{p\theta}\left(1-\frac{\xi^{2}}{5\theta}\right),$
(61)

$f^{\rm MB}$ is the Maxwell-Boltzmann distribution, $\xi_{i}$ is the peculiar
velocity, $p$ is the pressure, $\theta$ is the temperature, $\sigma_{ij}$ is
the traceless part of symmetric stress tensor and $q_{k}$ is the heat flux.
Thereafter, the leading terms of the two terms appearing in $\Theta$ are
evaluated as

$\displaystyle\begin{split}\left<f,x^{2}\right>&=\frac{1}{2p\theta}\sigma_{kl}\sigma_{lk}+{\cal
O}(\sigma^{3}),\\\
\left<f,x^{3}\right>&=-\frac{3}{p^{2}\theta}\sigma_{kl}\sigma_{lm}\sigma_{mk}+{\cal
O}(\sigma^{4}),\end{split}$ (62)

where assuming a small change in temperature ${\cal O}(q^{2})$ terms have been
ignored. Substituting $\sigma_{ij}=\rho\tau\theta S_{ij}$, $\bm{S}$ being the
strain rate tensor we find the viscosity correction as

$\nu_{T}=-{\tau\theta}\,\frac{\Delta
t}{2}\frac{S_{ij}S_{jk}S_{ki}}{S_{mn}S_{nm}-\tau S_{ab}S_{bc}S_{ca}}.$ (63)

It should be noted that for very fine grid resolutions ($\Delta t\rightarrow
0$) the viscosity correction vanishes. Similar expressions for the turbulent
viscosity have also been derived in Refs. Malaspinas  _et al._ (2008);
Buzzicotti and Tauzin (2021). The above expression for turbulent viscosity is
similar to Smagorinsky’s model where the turbulent viscosity $\nu_{T}$ is

$\nu_{T}=(C_{S}\Delta)^{2}\sqrt{S_{ij}S_{ji}},$ (64)

where $C_{S}$ is Smagorinsky’s constant, in that, both scale like the strain
rate tensor and is also distinct from it because of emergence of the third
invariant of the symmetrized strain rate tensor (Smagorinsky  _et al._ , 1965;
Deardorff, 1970).

##  VII Conclusion

In this paper, we present in detail the methodology to construct exact
solutions to the path length in the entropic lattice Boltzmann method. This
methodology can be extended to derive more accurate expressions, however, we
find that $\alpha_{\rm Higher}$ is sufficient for hydrodynamic applications.
The more dissipative solution $\alpha_{\rm Lower}$ could also be employed to
model viscous flows in the vicinity of walls and can also be used as a good
guess for the iterative solution. We have demonstrated that $\alpha_{\rm
Higher}$ shows no appreciable difference from the iterative solution by
studying the macroscopic behaviour of a few canonical setups. We have also
extended the exact solution to lattice ES-BGK model for nonlinear numerical
stability in non-unitary Prandtl heat transfer scenarios.

##  Appendix A Derivation of $\Delta H$

In this section, we derive the expression for $\Delta H=H[\bm{f}^{\rm
mirror}]-H[\bm{f}]$. We begin by using the form of $H$ [Eq. (3)] to obtain

$\displaystyle\begin{split}&H[\bm{f}^{\rm mirror}]-H[\bm{f}]=\left<f^{\rm
mirror},\log\frac{f^{\rm
mirror}}{w}\right>-\left<f,\log\frac{f}{w}\right>.\end{split}$ (65)

Substituting $\bm{f}^{\rm mirror}$ from Eq. (16) in the above equation yields

$\displaystyle\begin{split}&H[\bm{f}^{\rm mirror}]-H[\bm{f}]\\\
&=\left<f+\alpha(f^{\rm eq}-f),\log\frac{f+\alpha(f^{\rm
eq}-f)}{w}\right>-\left<f,\log\frac{f}{w}\right>.\end{split}$ (66)

Substituting $x$ from Eq. (19) in the above equation one obtains

$\displaystyle\begin{split}&H[\bm{f}^{\rm mirror}]-H[\bm{f}]\\\
&=\left<f(1+\alpha x),\log\frac{f(1+\alpha
x)}{w}\right>-\left<f,\log\frac{f}{w}\right>\\\ &=\left<f\left(1+\alpha
x\right),\log{\left(1+\alpha
x\right)}\right>-\alpha\left<fx,\log\frac{w}{f}\right>.\end{split}$ (67)

Now substituting $w_{i}$ from Eq. (8) one obtains

$\displaystyle\begin{split}&H[\bm{f}^{\rm mirror}]-H[\bm{f}]\\\
&=\left<f\left(1+\alpha x\right),\log{\left(1+\alpha x\right)}\right>\\\
&-\alpha\left<fx,\log\frac{f^{\rm
eq}\exp({\mu+\zeta_{\kappa}c_{i\kappa}+\gamma c_{i}^{2}})}{f}\right>\\\
&=\left<f\left(1+\alpha x\right),\log{\left(1+\alpha
x\right)}\right>-\alpha\left<fx,\log(1+x)\right>\\\
&-\underline{\alpha\lambda\left<f,x\right>}-\underline{\alpha\zeta_{\kappa}\left<f,xc_{\kappa}\right>}-\underline{\alpha\gamma\left<f,xc^{2}\right>},\end{split}$
(68)

where we have substituted $f_{i}^{\rm eq}/f_{i}=1+x_{i}$ and the underlined
terms are zero due to moments invariance, i.e.,

$\displaystyle\begin{split}\left<f,x\right>=\sum_{i}(f_{i}^{\rm
eq}-f_{i})=\rho-\rho=0,\\\ \left<f,xc_{\kappa}\right>=\sum_{i}(f_{i}^{\rm
eq}c_{i\kappa}-f_{i}c_{i\kappa})=\rho u_{\kappa}-\rho u_{\kappa}=0,\\\
\left<f,xc^{2}\right>=\sum_{i}(f_{i}^{\rm eq}c^{2}_{i}-f_{i}c^{2}_{i})=\rho
e-\rho e=0.\end{split}$ (69)

Thus, we obtain

$\displaystyle\begin{split}&H[\bm{f}^{\rm mirror}]-H[\bm{f}]\\\
&=\left<f\left(1+\alpha x\right),\log{\left(1+\alpha
x\right)}\right>-\alpha\left<fx,\log(1+x)\right>\\\ &=\left<f,\left(1+\alpha
x\right)\log{\left(1+\alpha
x\right)}\right>-\alpha\left<f,x\log(1+x)\right>.\end{split}$ (70)

##  Appendix B Bounds on the logarithm

In this section, we list a few positive definite functions along with their
domain of validity. In the interval $y\in(-1,0)$, using using the Taylor
series expansion of the logarithm we define

$\displaystyle G_{1}(y)$
$\displaystyle=(1+y)\left[-\log(1+y)+y-\frac{y^{2}}{2}\right]>0,$ (71)
$\displaystyle G_{4}(y)$
$\displaystyle=(1+y)\Bigg{[}-\log(1+y)+y-\frac{y^{2}}{2}+\frac{y^{3}}{3}-\frac{y^{4}}{4}$
$\displaystyle+\frac{y^{5}}{5}\Bigg{]}>0.$ (72)

Next, we exploit the integral definition of $\log(1+y)$, i.e.,

$\log(1+y)=\int_{0}^{y}\frac{1}{1+z}dz,$ (73)

and evaluate it using Gauss-Legendre and Newton-Cotes quadrature rules. As the
integrand is an $2n$-convex function, i.e., its even ($2n$) derivatives are
positive, the error due to the approximations are sign-definite, hence these
approximations can be used to construct upper and lower bounds on $\log(1+y)$.

Evaluating the integral in Eq. (73) via Gauss-Legendre quadratures, one
obtains

$\displaystyle\mathcal{I}_{\rm GL}^{(1)}(y)$ $\displaystyle=\frac{2y}{(2+y)},$
$\displaystyle\mathcal{I}_{\rm GL}^{(2)}(y)$
$\displaystyle=\frac{6y+3y^{2}}{6+6y+y^{2}},$ $\displaystyle\mathcal{I}_{\rm
GL}^{(3)}(y)$
$\displaystyle=\frac{60y+60y^{2}+11y^{3}}{60+90y+36y^{2}+3y^{3}},$

where $\mathcal{I}_{\rm GL}^{(n}$ is the intergral evaluated using $n^{\rm
th}$-order Gauss-Legendre quadrature. Similarly, evaluating the integral in
Eq. (73) via Newton-Cotes quadratures, one obtains Khattri (2009)

$\displaystyle\mathcal{I}_{\rm NC}^{(1)}(y)$
$\displaystyle=\frac{y}{2}\left[1+\frac{1}{1+y}\right],$
$\displaystyle\mathcal{I}_{\rm NC}^{(2)}(y)$
$\displaystyle=\frac{y}{6}\left[1+\frac{8}{2+y}+\frac{1}{1+y}\right],$
$\displaystyle\mathcal{I}_{\rm NC}^{(4)}(y)$
$\displaystyle=\frac{y}{90}\left[7+\frac{128}{4+y}+\frac{48}{4+2y}+\frac{128}{4+3y}+\frac{7}{1+y}\right],$

where $\mathcal{I}_{\rm GL}^{(n}$ is the intergral evaluated using $n^{\rm
th}$-order Newton-Cotes quadrature.

In the interval $y\in[0,\infty)$, exploiting the sign-definiteness of the
errors we define

$\displaystyle G_{2}(y)=(1+y)\left[-\log(1+y)+\mathcal{I}_{\rm
NC}^{(1)}\right]\geq 0,$ (74) $\displaystyle
G_{5}(y)=(1+y)\left[-\log(1+y)+\mathcal{I}_{\rm NC}^{(3)}\right]\geq 0,$ (75)

and in the interval $y\in(-1,\infty)$ we define

$\displaystyle G_{3}(y)=y\left[\log(1+y)-\mathcal{I}_{\rm GL}^{(1)}\right]\geq
0,$ (76) $\displaystyle G_{6}(y)=y\left[\log(1+y)-\mathcal{I}_{\rm
GL}^{(3)}\right]\geq 0.$ (77)

The functions $G_{1}(y),G_{2}(y),G_{3}(y)$ form loose bounds bounds on the
logarithm, whereas, $G_{4}(y),G_{5}(y),G_{6}(y)$ provide sharp bounds on it.

##  Appendix C Derivation of the higher-order solution

Following the same methodology as Section III.1, we add and subtract the same
terms from Eq. (67) to obtain

$\displaystyle\begin{split}\Delta
H&=H^{(B)}+\alpha\beta\hat{H}(\alpha),\end{split}$ (78)

where

$\displaystyle H^{(B)}=-\Big{<}f,{G_{6}(\alpha\beta
x)}\Big{>}_{\Omega^{-}}-\Big{<}f,{G_{7}(\alpha\beta
x)}\Big{>}_{\Omega^{+}}{-\alpha\beta G_{8}(x)}\leq 0,$ (79)

is nonpositive and contributes to the entropy production, and

$\displaystyle\begin{split}\hat{H}(\alpha)=-\left<f,\frac{\alpha^{2}\beta^{2}x^{3}}{6}-\frac{\alpha^{3}\beta^{3}x^{4}}{12}+\frac{\alpha^{4}\beta^{4}x^{5}}{20}-\frac{\alpha^{5}\beta^{5}x^{6}}{5}\right>_{\Omega^{-}}\\\
+\left<f,\frac{\alpha\beta
x^{2}}{2}\right>-\bigg{<}f,\frac{2\alpha^{2}\beta^{2}x^{3}}{15}\bigg{(}\frac{2}{4+\alpha\beta
x}+\frac{1}{4+2\alpha\beta x}\\\ +\frac{2}{4+3\alpha\beta
x}\bigg{)}\bigg{>}_{\Omega^{+}}-\left<f,\frac{60x^{2}+60x^{3}+11x^{4}}{60+90x+36x^{2}+3x^{3}}\right>.\end{split}$
(80)

The above equation has at least one positive root as
$\hat{H}(0)<0<\hat{H}(\infty)$ which can be found using any numerical method.
In order to preserve the computational efficiency of the method we solve the
above equation by converting it into a quadratic in $\alpha$.

###  C.1 Solving the higher degree polynomial

Figure 18: Behaviour of Eqs. (39),(80),(82),(45),(47) near the positive root.

In this section, we solve Eq. (80) by converting it to a quadratic. This
conversion to quadratic is performed by extracting negative terms from the Eq.
(80). The extracted terms then contribute to the entropy production $H^{(B)}$.
As stated earlier, the Eq. (80) has a positive root since
$\hat{H}(0)<0<\hat{H}(\infty)$. We assume that upper and lower bounds on the
root $\alpha$ exist. A suitable choice for the lower bound is $\alpha_{\rm
Lower}$, while the upper bound $h$ will be later evaluated. Therefore,
$\alpha_{\rm Lower}<\alpha<h$. Converting $\hat{H}(\alpha)$ to a quadratic is
a two step procedure and is explained in the following subsections.

####  C.1.1 Exploiting the lower bound

Using the lower bound $\alpha_{\rm Lower}$, in Eq. (80) we split the term

$\displaystyle\begin{split}-&\bigg{<}f,\frac{2\alpha^{2}\beta^{2}x^{3}}{15}\bigg{(}\frac{2}{4+\alpha\beta
x}+\frac{1}{4+2\alpha\beta x}+\frac{2}{4+3\alpha\beta
x}\bigg{)}\bigg{>}_{\Omega^{+}}\\\
&\equiv-\bigg{<}f,\frac{2\alpha\beta^{2}x^{3}}{15}\bigg{(}\frac{2}{\frac{4}{\alpha_{\rm
Lower}}+\beta x}+\frac{1}{\frac{4}{\alpha_{\rm Lower}}+2\beta
x}+\frac{2}{\frac{4}{\alpha_{\rm Lower}}+3\beta
x}\bigg{)}\bigg{>}_{\Omega^{+}}\\\
&-\bigg{<}f,\frac{2\alpha\beta^{2}x^{3}}{15}\bigg{(}\bigg{\\{}\frac{2}{\frac{4}{\alpha}+\beta
x}-\frac{2}{\frac{4}{\alpha_{\rm Lower}}+\beta
x}\bigg{\\}}+\bigg{\\{}\frac{1}{\frac{4}{\alpha}+2\beta
x}-\frac{1}{\frac{4}{\alpha_{\rm Lower}}+2\beta
x}\bigg{\\}}+\bigg{\\{}\frac{2}{\frac{4}{\alpha}+3\beta
x}-\frac{2}{\frac{4}{\alpha_{\rm Lower}}+3\beta
x}\bigg{\\}}\bigg{)}\bigg{>}_{\Omega^{+}},\end{split}$ (81)

where each term in curly braces is positive (as $\alpha_{\rm Lower}<\alpha$)
thereby making the second term negative. Here, recognizing that the negative
term contributes to the entropy production $H^{(B)}$, we obtain the quintic
polynomial $\tilde{H}(\alpha)$,

$\displaystyle\begin{split}\tilde{H}(\alpha)&=-\alpha^{2}\beta^{2}\left<f,\frac{x^{3}}{6}-\frac{\alpha\beta
x^{4}}{12}+\frac{\alpha^{2}\beta^{2}x^{5}}{20}-\frac{\alpha^{3}\beta^{3}x^{6}}{5}\right>_{\Omega^{-}}+\alpha\bigg{[}\bigg{<}f,\frac{x^{2}}{2}\bigg{>}-\bigg{<}f,\frac{2\alpha_{\rm
Lower}\beta^{2}x^{3}}{15}\bigg{(}\frac{2}{4+\alpha_{\rm Lower}x}\\\
&+\frac{1}{4+2\alpha_{\rm Lower}x}+\frac{2}{4+3\alpha_{\rm
Lower}x}\bigg{)}\bigg{>}_{\Omega^{+}}\bigg{]}-\left<f,\frac{60x^{2}+60x^{3}+11x^{4}}{60+90x+36x^{2}+3x^{3}}\right>.\end{split}$
(82)

Essentially, while converting $\hat{H}(\alpha)$ to $\tilde{H}(\alpha)$, we
have shifted the negative definite terms in Eq. (81) to the entropy
production, hence, the curve for $\tilde{H}(\alpha)$ lies above
$\hat{H}(\alpha)$ (see Fig. 18). It follows that an upper bound on the root of
$\hat{H}(\alpha)$ will also serve as the upper bound for the root of
$\tilde{H}(\alpha)$.

####  C.1.2 Exploiting the upper bound

Using the upper bound $h$, in Eq. (82) we split the term

$\displaystyle-\left<f,\frac{\alpha^{2}\beta^{2}x^{3}}{6}-\frac{\alpha^{3}\beta^{3}x^{4}}{12}+\frac{\alpha^{4}\beta^{4}x^{5}}{20}-\frac{\alpha^{5}\beta^{5}x^{6}}{5}\right>_{\Omega^{-}}\equiv-\alpha^{2}\beta^{2}\left<f,\frac{x^{3}}{6}-\frac{h\beta
x^{4}}{12}+\frac{h^{2}\beta^{2}x^{5}}{20}-\frac{h^{3}\beta^{3}x^{6}}{5}\right>_{\Omega^{-}}$
$\displaystyle-\alpha^{2}\beta^{2}\left<f,-\frac{(\alpha-h)\beta
x^{4}}{12}+\frac{(\alpha^{2}-h^{2})\beta^{2}x^{5}}{20}-\frac{(\alpha^{3}-h^{3})\beta^{3}x^{6}}{5}\right>_{\Omega^{-}},$
(83)

where the second term is negative, due to $x_{i}<0,x_{i}\in\Omega^{-}\,{\rm
and}\,\alpha<h$. Now, substituting Eq. (83) into Eq. (82) and again
recognizing that the negative terms contribute to the entropy production
$H^{(B)}$, we obtain the quadratic $H(\alpha)$.

It remains to specify the upper bound $h$. For this we consider the quadratic
equation $H_{2}(\alpha)=H(\alpha)|_{h=0}$,

$\displaystyle\begin{split}H_{2}(\alpha)&=-\alpha^{2}a_{2}+\alpha
b-c,\end{split}$ (84)
$\displaystyle\begin{split}a_{2}=\beta^{2}\left<f,\frac{x^{3}}{6}\right>_{\Omega^{-}},\end{split}$
(85)

whose positive root is $\alpha_{2}$. Therefore, $H_{2}(\alpha_{2})=0$ and

$\displaystyle H({\alpha_{2}})=\alpha_{2}^{2}\beta^{2}\left<f,\frac{h\beta
x^{4}}{12}-\frac{h^{2}\beta^{2}x^{5}}{20}+\frac{h^{3}\beta^{3}x^{6}}{5}\right>_{\Omega^{-}}$
$\displaystyle+H_{2}(\alpha_{2})>0.$ (86)

As $H(0)<0<H(\alpha_{2})$, a root of $H(\alpha)$ lies in the interval
$(0,\alpha_{2})$ (see Figure 18). Hence, a suitable choice for the upper bound
is $h=\alpha_{2}$.

##  Appendix D Implementing the analytical solution

The post-collisional populations are found via the routine

$f_{i}^{*}=f_{i}+\alpha\beta[f_{i}^{\rm eq}-f_{i}],$ (87)

where the path length $\alpha$ needs to be evaluated at each grid point. We
begin by calculating

$x_{i}=\frac{f_{i}^{\rm eq}}{f_{i}}-1,$ (88)

where $i=1\rightarrow N$ for a lattice with $N$ discrete velocities. To
evaluate a summation on one of the sub-divisions $\Omega^{-}$ or $\Omega^{+}$
we sum over the populations in the concerned subdivision. For instance, to
calculate

$a_{1}=\left<f,\frac{x^{3}}{2}\right>_{\Omega^{-}},\quad
b_{1}=\left<f,\frac{x^{2}}{2}\right>,$

the pseudo-code is:

1:$a_{1}=0,b_{1}=0$

2:for each integer $i$ in $1$ to $N$ do

3: if $x_{i}<0$ then

4: $a_{1}=a_{1}+f_{i}*x_{i}^{3}/2$

5: end if

6: $b_{1}=b_{1}+f_{i}*x_{i}^{2}/2$

7:end for

8:Return $a_{1}$

To find the path length $\alpha$ we execute the following steps:

1:Find $|x_{i}|^{\rm max}$, the $x_{i}$ with maximum magnitude.

2:if $|x_{i}|^{\rm max}<10^{-3}$ then

3: $\alpha=2$

4:else

5: Calculate $a_{1},b_{1},c_{1}$ from Eq. (40)

6: Calculate $\alpha_{\rm Lower}$ from Eq. (41)

7: Calculate $a_{2}$ from Eq. (85) and $b,c$ from Eq. (46)

8: Calculate $h$, the positive root of the Eq. (47)

9: Calculate $a,b,c$ from Eq. (46)

10: Find $\alpha_{\rm Higher}$, the positive root of Eq. (45)

11: $\alpha=\alpha_{\rm Higher}$

12:end if

Although, the exact solution to the path length is always found, we need to
ensure that the post collisional populations remain positive due to the
boundary conditions or in the case of of extremely under-resolved situations.
To this effect, an extra step might be required. We again stress that these
situations are extremely rare. The maximum permitted value of the path length
such that all the post collisional populations remain positive is $\alpha^{\rm
max}$. Therefore,

1:Find $x_{i}^{\rm min}$, the smallest $x_{i}$

2:Calculate $\alpha^{\rm max}=-1/(\beta x_{i}^{\rm min})$

3:if $\alpha>\alpha^{\rm max}$ then

4: $\alpha=(1+\alpha^{\rm max})/2$

5:end if

## References

  * Frisch  _et al._ (1986) U. Frisch, B. Hasslacher, and Y. Pomeau, Phys. Rev. Lett. 56, 1505 (1986).
  * Chen  _et al._ (1992) H. Chen, S. Chen, and W. H. Matthaeus, Phys. Rev. A 45, R5339 (1992).
  * Ansumali  _et al._ (2003) S. Ansumali, I. V. Karlin, and H. C. Öttinger, Europhys. Lett. 63, 798 (2003).
  * Yudistiawan  _et al._ (2010) W. P. Yudistiawan, S. K. Kwak, D. V. Patil, and S. Ansumali, Phys. Rev. E 82, 046701 (2010).
  * Adhikari  _et al._ (2005) R. Adhikari, K. Stratford, M. E. Cates, and A. J. Wagner, Europhys. Lett. 71, 473 (2005).
  * Mazloomi  _et al._ (2015) A. Mazloomi, S. S. Chikatamarla, and I. V. Karlin, Phys. Rev. Lett. 114, 174502 (2015).
  * Kolluru  _et al._ (2020a) P. K. Kolluru, M. Atif, and S. Ansumali, J. Comput. Sci. 45, 101179 (2020a).
  * Higuera  _et al._ (1989) F. J. Higuera, S. Succi, and R. Benzi, Europhys. Lett. 9, 345 (1989).
  * Qian  _et al._ (1992) Y. H. Qian, D. d’Humières, and P. Lallemand, Europhys. Lett. 17, 479 (1992).
  * Benzi  _et al._ (1992) R. Benzi, S. Succi, and M. Vergassola, Phys. Rep. 222, 145 (1992).
  * McNamara and Zanetti (1988) G. R. McNamara and G. Zanetti, Phys. Rev. Lett. 61, 2332 (1988).
  * Karlin  _et al._ (1999) I. V. Karlin, A. Ferrante, and H. C. Öttinger, Europhys. Lett. 47, 182 (1999).
  * Succi  _et al._ (2002) S. Succi, I. V. Karlin, and H. Chen, Rev. Mod. Phys. 74, 1203 (2002).
  * Boghosian  _et al._ (2001) B. M. Boghosian, J. Yepez, P. V. Coveney, and A. J. Wagner, Proc. R. Soc. London, Ser. A 457, 717 (2001).
  * Karlin  _et al._ (1998) I. V. Karlin, A. N. Gorban, S. Succi, and V. Boffi, Phys. Rev. Lett. 81, 6 (1998).
  * Wagner (1998) A. J. Wagner, Europhys. Lett. 44, 144 (1998).
  * Chen and Teixeira (2000) H. Chen and C. Teixeira, Comp. Phys. Commun. 129, 21 (2000).
  * Boghosian  _et al._ (2003) B. M. Boghosian, P. J. Love, P. V. Coveney, I. V. Karlin, S. Succi, and J. Yepez, Phys. Rev. E 68, 025103 (2003).
  * Ansumali  _et al._ (2006) S. Ansumali, I. Karlin, F. C.E., and K. Boulouchos, Physica A 359, 289 (2006).
  * Aidun and Clausen (2010) C. K. Aidun and J. R. Clausen, Annu. Rev. Fluid Mech. 42, 439 (2010).
  * Chikatamarla and Karlin (2013) S. Chikatamarla and I. Karlin, Physica A 392, 1925 (2013).
  * Atif  _et al._ (2017) M. Atif, P. K. Kolluru, C. Thantanapally, and S. Ansumali, Phys. Rev. Lett. 119, 240602 (2017).
  * Ansumali and Karlin (2000) S. Ansumali and I. V. Karlin, Phys. Rev. E 62, 7999 (2000).
  * Ansumali and Karlin (2002a) S. Ansumali and I. V. Karlin, J. Stat. Phys. 107, 291 (2002a).
  * Tosi  _et al._ (2006) F. Tosi, S. Ubertini, S. Succi, and I. V. Karlin, J. Sci. Comput. 30, 369 (2006).
  * Chikatamarla  _et al._ (2006) S. S. Chikatamarla, S. Ansumali, and I. V. Karlin, Phys. Rev. Lett. 97, 010201 (2006).
  * Brownlee  _et al._ (2007) R. A. Brownlee, A. N. Gorban, and J. Levesley, Phys. Rev. E 75, 036711 (2007).
  * Gorban and Packwood (2012) A. N. Gorban and D. Packwood, Phys. Rev. E 86, 025701 (2012).
  * Succi (2001) S. Succi,  _The Lattice Boltzmann Equation: for Fluid Dynamics and Beyond_ (Oxford University Press, Oxford, 2001).
  * Ansumali and Karlin (2005) S. Ansumali and I. V. Karlin, Phys. Rev. Lett. 95, 260605 (2005).
  * Bhatnagar  _et al._ (1954) P. L. Bhatnagar, E. P. Gross, and M. Krook, Phys. Rev. 94, 511 (1954).
  * Atif  _et al._ (2018) M. Atif, M. Namburi, and S. Ansumali, Phys. Rev. E 98, 053311 (2018).
  * Kolluru  _et al._ (2020b) P. K. Kolluru, M. Atif, M. Namburi, and S. Ansumali, Phys. Rev. E 101, 013309 (2020b).
  * Ansumali (2004) S. Ansumali,  _Minimal kinetic modeling of hydrodynamics_ , Ph.D. thesis, ETH Zurich (2004).
  * Gorban  _et al._ (1996) A. N. Gorban, I. V. Karlin, V. B. Zmievskii, and T. Nonnenmacher, Physica A 231, 648 (1996).
  * Mazloomi M.  _et al._ (2015) A. Mazloomi M., S. S. Chikatamarla, and I. V. Karlin, Phys. Rev. E 92, 023308 (2015).
  * Ansumali and Karlin (2002b) S. Ansumali and I. V. Karlin, Phys. Rev. E 65, 056312 (2002b).
  * Press  _et al._ (1992) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, The art of scientific computing 2, 1002 (1992).
  * Karlin  _et al._ (2015) I. Karlin, F. Bösch, S. Chikatamarla, and S. Succi, Entropy 17, 8099 (2015).
  * Laney (1998) C. B. Laney,  _Computational Gasdynamics_ (Cambridge University Press, 1998).
  * Minion and Brown (1997) M. L. Minion and D. L. Brown, J. Comput. Phys. 138, 734 (1997).
  * Coreixas  _et al._ (2017) C. Coreixas, G. Wissocq, G. Puigt, J. F. Boussuge, and P. Sagaut, Phys. Rev. E 96, 033306 (2017).
  * Ansumali and Karlin (2002c) S. Ansumali and I. V. Karlin, Phys. Rev. E 66, 026311 (2002c).
  * Ghia  _et al._ (1982) U. Ghia, K. Ghia, and C. Shin, J. Comput. Phys. 48, 387 (1982).
  * Holway Jr (1965) L. H. Holway Jr, Rarefied Gas Dyn. 1, 193 (1965).
  * Andries  _et al._ (2000) P. Andries, P. Le Tallec, J.-P. Perlat, and B. Perthame, Euro. J. Mech. B 19, 813 (2000).
  * Meng  _et al._ (2013) J. Meng, Y. Zhang, N. G. Hadjiconstantinou, G. A. Radtke, and X. Shan, J. Fluid Mech. 718, 347 (2013).
  * Kolluru  _et al._ (2022) P. K. Kolluru, M. Atif, and S. Ansumali, arXiv preprint arXiv:2201.05280 (2022).
  * Ansumali  _et al._ (2007) S. Ansumali, S. Arcidiacono, S. Chikatamarla, N. Prasianakis, A. Gorban, and I. Karlin, Eur. Phys. J. B 56, 135 (2007).
  * Shan (1997) X. Shan, Phys. Rev. E 55, 2780 (1997).
  * Togni  _et al._ (2015) R. Togni, A. Cimarelli, and E. De Angelis, J. Fluid Mech. 782, 380–404 (2015).
  * Karlin  _et al._ (2003) I. Karlin, S. Ansumali, E. De Angelis, H. Öttinger, and S. Succi, arXiv preprint cond-mat/0306003 (2003).
  * Malaspinas  _et al._ (2008) O. Malaspinas, M. Deville, and B. Chopard, Phys. Rev. E 78, 066705 (2008).
  * Buzzicotti and Tauzin (2021) M. Buzzicotti and G. Tauzin, Phys. Rev. E 104, 015302 (2021).
  * Smagorinsky  _et al._ (1965) J. Smagorinsky, S. Manabe, and J. L. Holloway, Mon. Weather Rev 93, 727 (1965).
  * Deardorff (1970) J. W. Deardorff, J. Fluid Mech. 41, 453 (1970).
  * Khattri (2009) S. Khattri, Teach. Math. 12, 7 (2009).