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

Calling entropy "lost information" does not explain why this run stays up. In the lattice Boltzmann method (LBM), lowering viscosity means pushing $\tau$ onto $1/2$, and BGK walks past equilibrium. Populations go negative. Karlin-style entropic LBM (ELBM) pins the discrete H to its pre-collision value right there. The price is one nonlinear root per node.

Drag the non-equilibrium perturbation slider to the right and the yellow dot leaves the dashed line at $\alpha=2$, so the second root of $G(\alpha)=0$ shifts left.

With a small perturbation the dot sits on 2. The gap between that dashed line and the yellow dot is the overshoot BGK refuses to drop.

When τ approaches 0.5, BGK falls apart#

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

$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 $\nu = c_s^2(\tau - 1/2)$. $c_s^2 = 1/3$ is the squared lattice sound speed. High-Reynolds flow, meaning low viscosity, forces $\tau \to 1/2$, so $\omega \to 2$.

As $\omega$ nears 2, the collision lands on the far side of equilibrium. A population that goes negative swings harder on the next step. The lattice blows up. BGK has nothing that stops this over-relaxation.

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

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

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

H is negative entropy. A monotone drop in H is the same statement as a monotone rise in entropy. Irreversibility. The arrow of time. The inequality is a stability certificate. A scheme in which H cannot grow cannot let the amplitude run away, because an unbounded amplitude would have to drag H up with it.

The discrete H-function and entropic equilibrium#

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

$w_i$ are the lattice weights: center $4/9$, faces $1/9$, diagonals $1/36$ on D2Q9. The log of the weight-normalized ratio.

Entropic equilibrium is the population that minimizes this $H$ under the conserved moments. Mass $\sum_i f_i = \rho$ and momentum $\sum_i \mathbf{c}_i f_i = \rho\mathbf{u}$ stay fixed while H is driven down. A Lagrange-multiplier problem. D2Q9 has a closed product form.

$u_\alpha$ is a velocity component and $c_{i\alpha}$ is a lattice-velocity component (−1, 0, 1). The second-order truncation in standard LBM is a polynomial equilibrium cut from a Hermite expansion, and that polynomial is not the minimizer of H. The product above is.

The entropic stabilizer — solving for α each step#

ELBM's collision has the same shape as BGK. The rate is split in two.

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

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

The update travels toward equilibrium only far enough to leave H identical to its pre-collision value. There is no room for H to grow. $\alpha=0$ is the trivial root. The root we need sits past equilibrium. Near equilibrium $\alpha \approx 2$, which is BGK again. Only a large departure pulls $\alpha$ off 2 and shaves the overshoot.

Numbers make the viscosity formula less abstract. With $c_s^2 = 1/3$ and $\nu = c_s^2(\tau - 1/2)$, $\tau = 0.51$ leaves a gap of $0.01$, so $\nu = 1/300 \approx 0.00333$. The rate is $\omega = 1/0.51 \approx 1.961$. A departure of length 1 is carried 1.961 of the way and lands 0.961 past equilibrium. At $\tau = 0.501$, $\omega \approx 1.996$. The leftover gap to the mirror value $\omega = 2$ is 0.004.

The listing's $\beta = 0.95$ is the same identity. $\tfrac{1}{2\beta}-\tfrac{1}{2} = \tfrac{1}{1.9}-\tfrac{1}{2} = \tfrac{1}{38}$, so $\nu = 1/114 \approx 0.00877$. BGK at that viscosity has $\tau = 1/(2\beta) \approx 0.526$. Near equilibrium, $\alpha \approx 2$ gives an effective rate $\alpha\beta = 1.9$, and the viscosity matches that BGK step. Once a large departure pulls $\alpha$ below 2, $\alpha\beta$ drops under 1.9 and the overshoot $\alpha\beta-1$ shrinks below 0.9. Viscosity stays where $\beta$ put it. Only the overshoot is cut.

Python — solving the entropic α by bisection#

Product-form equilibrium, the H-function, and one bisection for $\alpha$.

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 printed $\alpha$ is close to 2 and not exactly 2. That gap is the correction that holds H fixed.

The two perturbation lines are a moment shift. Adding $0.03$ to $f_1$ and subtracting $0.03$ from $f_3$ leaves the mass sum unchanged. Lattice speeds on those links are $+1$ and $-1$, so $x$-momentum rises by $0.03-(-0.03)=0.06$. From $\rho=1$, $u_x=0.05$, the state just after the kick has $u_x=0.11$. The product factor is $(2u+\sqrt{1+3u^2})/(1-u)$: about 1.162 at $u=0.05$, about 1.391 at $u=0.11$, and about 2.467 at nondimensional speed 0.3, where $1-u$ is only 0.7. The band $|u| \lesssim 0.1$ is where that factor is still mild.

The bisection cap is a sign constraint. A component with negative $d_i$ hits zero at $\alpha=-f_i/d_i$, and the smallest of those values is $a_{\max}$. The search stops at $\min(0.999\,a_{\max},\,4)$. A root past 4, or one that flips a population's sign, is not a candidate.

Interactive — the second root of G(α)#

The curve above is this $G$. The yellow dot is the non-trivial root $\alpha^*$. A larger perturbation peels it off 2. BGK nails that slot to 2. Stack the difference for enough steps and the lattice diverges.

Interactive — blocking the moment H grows#

Feed the same non-equilibrium state to both collision rules and only the populations that went negative turn into red bars.

The entropic side stays on $H_{ent} = H_{pre}$. That is the construction. On the BGK side, once the departure grows, $H_{BGK}$ crosses the pre-collision H. The red warning lights there. That rise in H is where a low-viscosity lattice starts to rupture.

Three things I hit while implementing#

A nonlinear solve at every node costs more than the viscosity term is worth. Near equilibrium, when $\|f - f^{eq}\|$ is small, set $\alpha \approx 2$ and turn bisection on only where the departure crosses a threshold. The cost falls to about a tenth.

The product equilibrium is defined only for $|u| < 1$. The denominator carries $1 - u$. Past a nondimensional speed of 0.3 the pressure tensor starts to skew by direction. Even ELBM is safer inside $|u| \lesssim 0.1$.

Keep the second-order polynomial as the equilibrium and the H-minimizer guarantee is gone. The stabilizer pays off fully only when the equilibrium itself is the entropic product. Bolting $\alpha$ onto a polynomial gets half of that.

Locking α at 2 while pushing β undoes the H bound#

Skipping to $\alpha \approx 2$ near equilibrium is a cost cut for small $\|f - f^{eq}\|$. Push $\beta$ toward 1 and freeze that shortcut on every node, and the effective rate $\alpha\beta$ sits on 2 again. The overshoot out to the mirror comes back, and discrete H can exceed its pre-collision value. Viscosity dropped. The root was never solved.

A collision can hold H flat and still be late if streaming imports a population that is already negative from a neighbor. Whether a sign broken in advection still needs a positivity-preserving flux limiter, or whether the entropic root is enough on its own, is the next cut.