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

In 1986, Frisch, Hasslacher, and Pomeau dropped particles on a lattice and made them mimic Navier–Stokes. That was LGCA, the lattice gas cellular automaton. It died within five years — the integer occupancy of zero or one per direction generated noise that no amount of averaging could erase. In 1992, McNamara and Zanetti swapped integers for real numbers. They populated the lattice with particle distribution functions $f_i$ instead of bits. The noise vanished — but a new question rode in. How can nine real numbers stand in for the infinite-dimensional Maxwell–Boltzmann distribution? This post follows the answer: a second-order Hermite truncation gives LBM its polynomial equilibrium, and that same truncation locks it into low-Mach flow.

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

LGCA was elegant. Square lattice, six particle directions, one collision rule. Average it macroscopically and Navier–Stokes returns. But boolean particles (zero or one) amplified statistical noise. You could kill the noise with long time averages, but then the transient you came for was gone too.

McNamara and Zanetti placed real-valued distributions $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 lattice direction $\mathbf{c}_i$. The noise problem vanished by construction. One step splits into two halves — collision and streaming.

$\tau$ is the relaxation time and $f_i^{eq}$ is the equilibrium distribution. Collision shaves the non-equilibrium part on a $\tau$ timescale, streaming shifts the shaved distribution one cell over. The deep question hides in $f_i^{eq}$ — what makes nine real numbers enough?

Maxwell–Boltzmann carries a time bomb#

In continuum kinetic theory, the equilibrium distribution is Maxwell–Boltzmann.

$\mathbf{c}$ is the microscopic particle velocity, $\mathbf{u}$ is the macroscopic flow speed, $D$ is the spatial dimension, $T$ the temperature, $R$ the gas constant. The bomb is that $\mathbf{c}$ is continuous. Reducing it to nine or nineteen discrete directions requires you to skim the distribution intelligently.

Naive attempt — replace the Gaussian integral with a midpoint rule on a fixed grid. It fails. The moments that Navier–Stokes needs — $\rho$, $\rho \mathbf{u}$, the momentum tensor $\rho \mathbf{u}\mathbf{u}$, energy — disagree with the lattice sum. Once conservation breaks, LBM is dead in the water.

The right answer is series expansion. Expand Maxwell–Boltzmann in orthogonal polynomials, keep the orders that carry the needed moments, drop the rest. Those polynomials are Hermite's.

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 1D Maxwell–Boltzmann to be isothermal (rescale to $RT = 1$). The exponential generating function gives this identity straight away.

Each term grabs a slice of the universe. $n=0$ is the static Gaussian. $n=1$ is a small drift. $n=2$ is the mean-kinetic-energy correction. Higher $n$ are increasingly asymmetric tail effects. LBM's whole game is to cut this sum at a finite order.

Play with the simulation above. Push the drift $u$ to 0.4 and pick truncation order $N=2$ — the first three moments (density, momentum, energy) match Maxwell–Boltzmann exactly. Crank $u$ past 1.0 and the dashed curve peels off. That is exactly where LBM's Mach-number constraint $|\mathbf{u}|/c_s \lesssim 0.3$ comes from, drawn in one panel.

What second-order truncation buys you#

Which moments do you need for first-order viscous Navier–Stokes? Three orders, no more — $\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 second-order in $u$, the distribution function must carry a $u^2$ term.

So cut the Hermite series at $N=2$.

This is the continuous form of the LBM equilibrium. $N=3$ and beyond help with higher inviscid moments (heat flux) but are surplus to standard isothermal LBM. The dropped terms are responsible for LBM's two well-known limits — the isothermal assumption and the 0.3-Mach ceiling.

How D2Q9 was built#

Replacing the continuous integral $\int f^{eq}(c) \phi(c)\, dc$ with a lattice sum $\sum_i w_i f_i^{eq} \phi(\mathbf{c}_i)$ demands quadrature. The key fact — quadrature accurate to order $N$ requires only that it integrate polynomials of that order exactly.

In 2D, accuracy through the $u^2$ term means quadrature exact for polynomials up to fifth order under the Gaussian weight (the momentum tensor produces a $c^2 \times u^2$ term). Three-point Gauss–Hermite does precisely that job.

In 1D, three nodes $c \in \{-\sqrt{3}, 0, +\sqrt{3}\}$ with weights $\{1/6, 2/3, 1/6\}$. The tensor product in 2D gives $3 \times 3 = 9$ nodes — that is D2Q9. Rescale to lattice units and you get $c_s = 1/\sqrt{3}$ with lattice velocity vectors at integer coordinates $(0,0), (\pm 1, 0), (0, \pm 1), (\pm 1, \pm 1)$. The lattice weights $w_i$ absorb the Gaussian factor into the Gauss–Hermite weights.

Plugging $c_s^2 = 1/3$ into the continuous form gives the on-lattice equilibrium.

That one line is the heart of every isothermal LBM code. A simple polynomial — but the first three moments of Maxwell–Boltzmann are sealed inside it.

Python — one BGK collision step#

Run a single collision step on one lattice 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}")

About twenty steps and the residual drops below $10^{-6}$ — the nine $f_i$ settle onto the equilibrium polynomial. What survives each step? Density $\rho$ and momentum $\rho \mathbf{u}$. Collision only chips away non-equilibrium moments.

Pull the $\tau$ slider down toward 0.55 and the residual vanishes in one or two steps. Viscosity is $\nu = c_s^2 (\tau - 1/2)$ — as $\tau \to 1/2$, viscosity disappears and the lattice goes unstable. Push $\tau$ up and viscosity grows but relaxation slows. That trade-off is the dial every LBM practitioner spends time on.

Three lines to remember#

• LBM's equilibrium $f_i^{eq}$ is Maxwell–Boltzmann expanded in Hermite polynomials and chopped at second order — that cut preserves the three moments $\rho$, $\rho \mathbf{u}$, and the momentum flux exactly.
• D2Q9 is derived as the 2D tensor product of three-point Gauss–Hermite quadrature, the smallest rule that integrates the second-order equilibrium exactly. Nine is not an arbitrary number; it is an artifact of orthogonal integration.
• The truncation order sets every limit. Isothermal assumption, $|\mathbf{u}|/c_s \lesssim 0.3$ ceiling — both follow automatically from the single line "we cut the series at $N=2$".