Nine Arrows of Maxwell–Boltzmann — Why LBM's Equilibrium Is a Polynomial

Lattice gas died in 1986. Occupancy was zero or one, and the noise would not average out. In 1992 those integers became a real distribution $f_i$. The noise went away. The question stayed. How do nine real numbers stand in for infinite-dimensional Maxwell–Boltzmann?

Move drift $u$ and truncation order $N$ and the solid Maxwell–Boltzmann curve either sits on the dashed polynomial or pulls away from it.

At $N=2$ and $u=0.4$, density, momentum, and energy lie on top of each other. Past $u=1$ only the dashed curve leaves. The Mach ceiling $|\mathbf{u}|/c_s \lesssim 0.3$ is not that peel-off. The numbers come later.

From integers to reals — LGCA died in 1986 and LBM was born#

LGCA, the lattice gas cellular automaton, was a square lattice, six directions, one collision rule. The macroscopic average came back as Navier–Stokes. Boolean particles (zero or one) grew statistical noise. A longer time average kills the noise and the transient with it.

McNamara and Zanetti put a real distribution $f_i(\mathbf{x}, t)$ on every node. $f_i$ is the expected count of particles at position $\mathbf{x}$ and time $t$ moving along $\mathbf{c}_i$. One step splits in two. Collision, then streaming.

$\tau$ is the relaxation time. $f_i^{eq}$ is the equilibrium distribution. Collision shaves the non-equilibrium part on a $\tau$ scale, and streaming shifts that distribution by one cell. Nine numbers are enough only because of the shape of $f_i^{eq}$.

Maxwell–Boltzmann carries a time bomb#

In continuum kinetic theory the equilibrium is Maxwell–Boltzmann.

$\mathbf{c}$ is the molecular velocity, $\mathbf{u}$ the macroscopic speed, $D$ the dimension, $T$ the temperature, $R$ the gas constant. The bomb is that $\mathbf{c}$ is continuous. Keeping nine or nineteen directions means thinning this distribution.

A rectangle-rule stand-in for the Gaussian integral fails. The moments Navier–Stokes asks for — $\rho$, $\rho\mathbf{u}$, the momentum tensor $\rho\mathbf{u}\mathbf{u}$, energy — miss the lattice sum. Break conservation and LBM stops there.

The way through is a series. Expand Maxwell–Boltzmann in orthogonal polynomials and keep only the moments you need. Those polynomials are Hermite.

Hermite polynomials chop the distribution into pieces#

The probabilist Hermite polynomials $He_n(c)$ are orthogonal under the standard Gaussian weight $w(c)=(2\pi)^{-1/2}e^{-c^2/2}$.

The first few are $He_0=1$, $He_1=c$, $He_2=c^2-1$, $He_3=c^3-3c$. Set the one-dimensional equilibrium to isothermal, $RT=1$, and the generating function gives this.

$n=0$ is a Gaussian at rest. $n=1$ is a small drift. $n=2$ corrects the mean kinetic energy. Higher $n$ are asymmetric tails. LBM cuts the sum at a finite order.

In that isothermal expansion $RT=1$, so the thermal speed is 1. The $u$ on the figure is already divided by that speed. The quoted ceiling, 0.3, is $|u|$ on this scale. $u=0.4$ sits 0.1 above the ceiling, and the first three moments at $N=2$ still match. The peel-off at $|u|>1$ is past Mach 1. The 0.3 and the 1 are not the same failure.

The $n=3$ term is $(u^3/6)\,He_3(c)$, with $He_3=c^3-3c$. The coefficient alone is $0.3^3/6=0.027/6=0.0045$ at $u=0.3$. At $u=0.4$ it is $0.4^3/6=0.064/6\approx 0.0107$, about 2.37 times the ceiling value. At $u=1$ it is $1/6\approx 0.1667$, about 37 times 0.0045. That ratio is $1^3/0.3^3=1/0.027$. $N=2$ throws the whole term away. So the density, momentum, and energy can still match at $u=0.4$ while the discarded cubic is already more than twice its size at the ceiling. The peel-off at $u=1$ is where that coefficient has climbed to about 0.17. The Mach ceiling at 0.3 is not the point where the curve visibly leaves. It is a line drawn in advance, at the moment the cubic coefficient passes 0.0045.

What second-order truncation buys you#

First-order viscous Navier–Stokes needs three layers. $\rho$, $\rho\mathbf{u}$, and the momentum flux $\Pi_{\alpha\beta}=\rho u_\alpha u_\beta+p\delta_{\alpha\beta}$. Because $\rho\mathbf{u}\mathbf{u}$ is $u^2$, the distribution has to carry a $u^2$ term.

Cut the series at $N=2$.

That is the continuous equilibrium. $N=3$ and above feed the heat flux. Standard isothermal LBM does not keep them. The discarded terms are what leave both the isothermal assumption and the Mach 0.3 ceiling. Put the cubic back and the isothermal assumption breaks first. Heat flux is the piece that term carries, so a code stopped at second order cannot take a problem where temperature flows and still use this polynomial. The Mach ceiling and the isothermal assumption are two faces of one cut. Even with speed held under 0.173, a hot wall and a cold wall on the same mesh get no heat transport from this equilibrium.

How D2Q9 was built#

Turning the integral $\int f^{eq}(c)\,\phi(c)\,dc$ into the lattice sum $\sum_i w_i f_i^{eq}\phi(\mathbf{c}_i)$ needs a quadrature. Matching through order $N$ means integrating polynomials of that order exactly.

In two dimensions, accuracy through $u^2$ means exactness for fifth-order polynomials under the Gaussian weight. The momentum tensor produces a $c^2\times u^2$ term. Three-point Gauss–Hermite does that integral.

The one-dimensional nodes are $c\in\{-\sqrt{3},0,+\sqrt{3}\}$, with weights $\{1/6,2/3,1/6\}$. Why those three roots sit where they do is counted in the Hermite quadrature that folded a sixth moment down onto nine speeds. The two-dimensional tensor product is $3\times 3=9$. That is D2Q9. In lattice units $c_s=1/\sqrt{3}$, and the velocities land on integer coordinates $(0,0),(\pm 1,0),(0,\pm 1),(\pm 1,\pm 1)$. The weights $w_i$ are the Gauss–Hermite weights with the Gaussian factor absorbed.

On-lattice equilibrium is the continuous form with $c_s^2=1/3$ substituted.

$1/c_s^2=3$, $1/(2c_s^4)=9/2$, $1/(2c_s^2)=3/2$. The 3, the $9/2$, and the $3/2$ in the polynomial are those three numbers. At $\mathbf{u}=0$ the bracket is 1 and $f_i^{eq}=w_i\rho$. The weights sum to $4/9+4\cdot(1/9)+4\cdot(1/36)=1$. Four ninths of a fluid at rest sits on direction 0. The three table weights are products of the one-dimensional set $\{1/6,\ 2/3,\ 1/6\}$. Rest is $(2/3)\times(2/3)=4/9$. An axis direction is $(2/3)\times(1/6)=1/9$. A diagonal is $(1/6)\times(1/6)=1/36$. Four axes and four diagonals close the sum again, $4/9+4\cdot(1/9)+4\cdot(1/36)=1$. The count of nine and the three denominators come out of the same tensor product.

That line is the equilibrium in an isothermal LBM code. A polynomial, with the first three moments of Maxwell–Boltzmann locked inside it.

Python — one BGK collision step#

One collision step, one node.

import numpy as np
 
cx = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
cy = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w  = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
 
def equilibrium(rho, ux, uy):
    """D2Q9 equilibrium — nine polynomials"""
    ciu = cx * ux + cy * uy            # dot product
    usq = ux * ux + uy * uy
    return w * rho * (1 + 3 * ciu + 4.5 * ciu**2 - 1.5 * usq)
 
def bgk_collision(f, tau):
    rho = f.sum()
    ux  = (cx * f).sum() / rho
    uy  = (cy * f).sum() / rho
    feq = equilibrium(rho, ux, uy)
    return f - (f - feq) / tau, rho, ux, uy
 
# off-equilibrium initial state
f = w * 1.0
f[1] += 0.05; f[3] -= 0.04          # slight eastward bias
 
for step in range(20):
    f, rho, ux, uy = bgk_collision(f, tau=0.8)
    err = np.linalg.norm(f - equilibrium(rho, ux, uy))
    print(f"step {step:2d}: rho={rho:.4f}, u=({ux:+.4f},{uy:+.4f}), ||delta||={err:.2e}")

The initial state starts from $f=w$, adds 0.05 to $f_1$, and subtracts 0.04 from $f_3$. Mass is $1+0.05-0.04=1.01$. The $x$-momentum is $1\cdot 0.05+(-1)\cdot(-0.04)=0.09$. So $u_x=0.09/1.01\approx 0.0891$ and $u_y=0$. Divide by $c_s=1/\sqrt{3}$ and the Mach number is about 0.154. The lattice speed at the 0.3 ceiling is $0.3/\sqrt{3}\approx 0.173$, and 0.0891 is under it. The residual the run drives below $10^{-6}$ in about twenty steps is the approach to equilibrium on this conserved pair $(1.01,\ 0.0891)$. Each collision keeps only $\rho$ and $\rho\mathbf{u}$.

Drop $\tau$ toward 0.55 and the residual is gone in a step or two. Viscosity is $\nu=c_s^2(\tau-1/2)$. At the code's $\tau=0.8$, $\nu=(1/3)\cdot 0.3=0.1$. At $\tau=0.55$, $\nu=(1/3)\cdot 0.05=1/60\approx 0.0167$, one sixth of the code value. Push $\tau$ all the way to $1/2$ and viscosity hits zero. The lattice goes unstable. Raise $\tau$ and viscosity grows while the pull toward equilibrium slows down. One step multiplies the non-equilibrium part by $1-1/\tau$. At the code's $\tau=0.8$, $1/\tau=1.25$, so the factor is $-0.25$. The sign flips. That is over-relaxation. Twenty factors of $0.25$ give $0.25^{20}\approx 9.1\times 10^{-13}$. A deviation of order 1 is already far under $10^{-6}$ after twenty steps. Taking twenty steps as the mark for $10^{-6}$ is generous on this factor alone. At $\tau=1$ the factor is 0 and the state lands on equilibrium in one step. The code's 0.8 stands one notch off that single-step equilibrium.

Past Mach 0.3, do not trust this polynomial#

$N=2$ matches $\rho$, $\rho\mathbf{u}$, and the momentum flux. Everything above that was dropped. Once the lattice speed passes about 0.173, the dropped terms are no longer small. Lowering $\tau$ so viscosity falls from $0.1$ to $1/60$ does not change the equilibrium polynomial. Raising the speed through the ceiling and lowering the viscosity are different failures.

Whether a viscosity tuned in a rest frame survives a shift of that frame is not answered here. The discarded cubic term holds that answer. Whether the missing third moment changes viscosity once the frame moves still sits outside this polynomial.