Entropy Stops the Blow-Up — Karlin's Entropic Lattice Boltzmann

Entropy usually arrives labeled as "lost information" — the direction disorder grows, a measure of irreversible loss. Yet in the Lattice Boltzmann Method (LBM), that very quantity is what catches a runaway simulation and holds it still. This post is about how Karlin-style Entropic LBM (ELBM) enforces the H-theorem at the discretization level to kill blow-up in the low-viscosity regime. By the end you'll see why a single line of the collision operator is worth solving a nonlinear equation at every step.

When τ approaches 0.5, BGK falls apart#

Standard LBM uses the BGK collision (a single-relaxation-time approximation). It pulls the distribution $f_i$ toward equilibrium $f_i^{eq}$ by a factor $\omega$.

Here $f_i$ is the population in lattice direction $i$, $f_i^{eq}$ is the equilibrium, and $\omega = 1/\tau$ is the relaxation rate.

Viscosity is set by $\nu = c_s^2(\tau - 1/2)$, where $c_s^2 = 1/3$ is the squared lattice sound speed. To resolve high-Reynolds (low-viscosity) flow, you push $\tau \to 1/2$, i.e. $\omega \to 2$.

The trouble is that as $\omega$ nears 2, the collision overshoots equilibrium and lands on the far side. Populations dip negative, and once an $f_i$ goes negative it oscillates harder on the next step. The lattice blows up. BGK has nothing to stop this over-relaxation.

Boltzmann's H-theorem — the arrow of time#

In 1872 Boltzmann proved that the following quantity never grows under collisions alone.

$f$ is the velocity distribution and $\mathbf{c}$ is the molecular velocity. We have $dH/dt \le 0$, with equality only at equilibrium (the Maxwell-Boltzmann distribution).

H is negative entropy. A monotone decrease of H is the same statement as a monotone increase of entropy. That is irreversibility — the arrow of time. The key point is that this inequality is a stability certificate. A numerical scheme in which H never grows cannot diverge, because unbounded amplitude growth would require H to increase.

The discrete H-function and entropic equilibrium#

LBM slices the continuous velocity into 9 directions (D2Q9). H is discretized to match.

$w_i$ are the lattice weights (in D2Q9: center $4/9$, faces $1/9$, diagonals $1/36$). It is the log of the weight-normalized ratio.

The entropic equilibrium is defined as the distribution that minimizes this $H$ under the conserved-moment constraints. Holding mass $\sum_i f_i = \rho$ and momentum $\sum_i \mathbf{c}_i f_i = \rho\mathbf{u}$ fixed, you minimize H (a Lagrange-multiplier problem). In D2Q9 a closed product form exists.

$u_\alpha$ is a velocity component and $c_{i\alpha}$ is a lattice-velocity component (−1, 0, 1). Unlike the second-order truncated polynomial equilibrium of standard LBM, this one is the true minimizer of H.

The entropic stabilizer — solving for α each step#

ELBM's collision has the same shape as BGK but splits the relaxation rate in two.

$\beta \in (0,1)$ sets the viscosity ($\nu = c_s^2(\tfrac{1}{2\beta} - \tfrac{1}{2})$, with $\beta\to1$ being inviscid), and $\alpha$ sets the stability.

$\alpha$ is not a constant. At every step, at every lattice node, it is solved as the non-trivial root of

The meaning is plain: relax toward equilibrium by exactly the amount that leaves H identical to its pre-collision value. Any increase of H is structurally forbidden. $\alpha=0$ is the trivial root; what we want is the second root beyond equilibrium. Near equilibrium $\alpha \approx 2$, which is exactly BGK. Only when the non-equilibrium part is large does $\alpha$ depart from 2 and shave off the over-relaxation.

Python — solving the entropic α by bisection#

Here is the product-form equilibrium, the H-function, and one step that solves $\alpha$ by bisection.

import numpy as np
 
# D2Q9 lattice weights & velocities
W  = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
CX = np.array([0, 1, 0, -1,  0, 1, -1, -1,  1])
CY = np.array([0, 0, 1,  0, -1, 1,  1, -1, -1])
 
def h_function(f):
    # discrete H-function: H = Σ f_i ln(f_i / w_i)  (lattice version of Boltzmann's H)
    return np.sum(f * np.log(f / W))
 
def entropic_equilibrium(rho, ux, uy):
    # exact Ansumali–Karlin product-form entropic equilibrium
    def factor(u):
        s = np.sqrt(1 + 3 * u * u)
        return (2 - s), (2 * u + s) / (1 - u)
    bx, rx = factor(ux)
    by, ry = factor(uy)
    return rho * W * bx * by * (rx ** CX) * (ry ** CY)
 
def entropic_alpha(f, feq):
    # non-trivial root α* of H(f + α(feq − f)) = H(f), via bisection
    d = feq - f
    H0 = h_function(f)
    G = lambda a: h_function(f + a * d) - H0
    amax = np.min(np.where(d < 0, -f / d, np.inf))   # positivity bound
    hi = min(0.999 * amax, 4.0) if np.isfinite(amax) else 4.0
    if G(hi) <= 0:
        return 2.0                                   # near-equilibrium safe value
    lo = 1.0
    while G(lo) > 0 and lo > 1e-3:
        lo *= 0.5
    for _ in range(60):
        mid = 0.5 * (lo + hi)
        if G(mid) > 0: hi = mid
        else:          lo = mid
    return 0.5 * (lo + hi)
 
# --- one node, one step ---
feq = entropic_equilibrium(1.0, 0.05, 0.0)
f = feq.copy()
f[1] += 0.03; f[3] -= 0.03            # non-equilibrium carried in by streaming
rho = f.sum(); ux = (CX * f).sum() / rho; uy = (CY * f).sum() / rho
feq = entropic_equilibrium(rho, ux, uy)
 
alpha = entropic_alpha(f, feq)
beta  = 0.95                          # ν = c_s²(1/(2β) − 1/2), β→1 is low viscosity
f_post = f + alpha * beta * (feq - f)
 
print(f"entropic α  = {alpha:.4f}   (BGK fixed value 2.0)")
print(f"H before    = {h_function(f):.8f}")
print(f"H after     = {h_function(f_post):.8f}")
print(f"min pop     = {f_post.min():.3e}  (> 0)")

The output shows $\alpha$ close to but not exactly 2. That tiny gap is the correction that conserves H.

Interactive — the second root of G(α)#

Play with the simulation below. Increase the non-equilibrium perturbation and watch the curve bend.

The yellow dot is the $\alpha^*$ we want. For small perturbations it sits right on the dashed line ($\alpha=2$). Push the slider right and $\alpha^*$ peels away from 2. BGK ignores this shift and always uses 2. When the gap accumulates, the lattice diverges.

Interactive — blocking the moment H grows#

The same non-equilibrium state, processed by two collision rules side by side. Red bars are populations that went negative.

The entropic side always has $H_{ent} = H_{pre}$ — that is how it is built. On the BGK side, once the non-equilibrium part grows, a moment arrives where $H_{BGK}$ exceeds the pre-collision value. That is when the red warning lights up. That increase in H is exactly the seed that ruptures the lattice at low viscosity.

Three things I hit while implementing#

First, a nonlinear solve at every node is expensive. So in practice you set $\alpha \approx 2$ near equilibrium (when $\|f - f^{eq}\|$ is small) and turn on bisection only at nodes where the non-equilibrium part crosses a threshold. The cost drops to roughly a tenth.

Second, the product-form equilibrium is only defined for $|u| < 1$ — there is a $1 - u$ in the denominator. Beyond a non-dimensional speed of 0.3 the pressure tensor begins to skew anisotropically, so even ELBM is safest in the $|u| \lesssim 0.1$ regime.

Third, if you keep the second-order truncated polynomial equilibrium, the H-minimization guarantee breaks. To get the full benefit of the entropic stabilizer the equilibrium itself must be the entropic product form. Bolting $\alpha$ onto a polynomial equilibrium gets you only half the effect.

Three-line takeaway#

• BGK blows up at low viscosity because over-relaxation increases the discrete H. ELBM prevents this by solving, each step, the $\alpha$ that makes H after collision equal to H before.
• $\alpha$ is the non-trivial root of $H(f + \alpha(f^{eq}-f)) = H(f)$, and near equilibrium $\alpha \approx 2$ automatically recovers BGK.
• Viscosity is carried by $\beta$, stability by $\alpha$. The cost is a per-step nonlinear solve, but the near-equilibrium shortcut keeps it affordable in practice.