The double-distribution energy backend

The double-distribution-function (DDF) energy backend transports the conserved energy variable on its own lattice distribution \(g_i\), as an alternative to the finite-difference (FD) backend of the energy equation page.

Both backends solve the same energy model - the same conserved variable, the same coupling to the fluid through the equation of state and buoyancy. They differ only in how the energy is transported in space: the FD backend integrates the enthalpy material derivative by finite differences; the DDF backend streams a second lattice population. The DDF route is conservative by construction and shares the fluid’s discrete operator, which removes the FD\(\leftrightarrow\)LBM coupling instability that the FD backend controls with an added pressure filter.

What must be recovered

The target is the variable-density low-Mach energy equation the FD backend already integrates, Eq. (7) of the energy-equation page: the material derivative of specific enthalpy balances conduction and a volumetric heat release,

(1)\[ \rho\,\frac{\mathrm{D}h}{\mathrm{D}t} = \partial_\alpha\!\left(\lambda\,\partial_\alpha T\right) + Q . \]

With a constant specific heat \(c_p\) and the caloric law \(h = c_p\,(T - T_{\text{ref}})\) (datum \(h=0\) at \(T_{\text{ref}}\)), so that \(\partial_\alpha T = \partial_\alpha h / c_p\), and the conductivity \(\lambda = c_p\,\mu(T)/\mathrm{Pr}\) at a fixed Prandtl number, the conduction term reduces to a diffusion of enthalpy with coefficient \(D = \mu(T)/\mathrm{Pr}\),

(2)\[ \partial_\alpha\!\left(\lambda\,\partial_\alpha T\right) = \partial_\alpha\!\left(\frac{\mu(T)}{\mathrm{Pr}}\,\partial_\alpha h\right) = \partial_\alpha\!\left(D\,\partial_\alpha h\right), \qquad D = \frac{\mu(T)}{\mathrm{Pr}} . \]

This is the same reduction the FD backend uses (Eq. (9)); the DDF backend must reproduce it.

The continuity identity: the material derivative is exactly conservative

The single fact that makes an LBM energy transport natural is that, given continuity, the material derivative of enthalpy equals the conservative (divergence) form exactly, with no approximation. Writing continuity as \(\partial_t\rho + \partial_\alpha(\rho u_\alpha) = 0\),

(3)\[ \partial_t(\rho h) + \partial_\alpha(\rho h\,u_\alpha) = \rho\,\frac{\mathrm{D}h}{\mathrm{D}t} + h\,\underbrace{\left[\partial_t\rho + \partial_\alpha(\rho u_\alpha)\right]}_{=\,0\ \text{(continuity)}} , \]

so that

(4)\[ \rho\,\frac{\mathrm{D}h}{\mathrm{D}t} = \partial_t(\rho h) + \partial_\alpha(\rho h\,u_\alpha) . \]

The target equation (1) is therefore identical to the conservative energy equation

(5)\[ \partial_t(\rho h) + \partial_\alpha(\rho h\,u_\alpha) = \partial_\alpha\!\left(D\,\partial_\alpha h\right) + Q, \qquad D = \frac{\mu(T)}{\mathrm{Pr}} . \]

The conservative and non-conservative forms target the same PDE

Equation (4) is an identity, not a modelling choice: the FD backend’s non-conservative form \(u_\alpha\partial_\alpha h\) and the DDF backend’s conservative form \(\partial_\alpha(\rho h\,u_\alpha)\) are the same continuous equation. The FD backend does not solve a different physics. The difference is one of discretisation fidelity at finite dilatation \(\partial_\alpha u_\alpha\): a first-order upwind discretisation of \(u_\alpha\partial_\alpha h\) can mis-represent the dilatation contribution where gradients are under-resolved Wissocq and Sagaut[1], whereas the DDF carries \(\partial_\alpha(\rho h\,u_\alpha)\) in conservative flux form and represents it faithfully. In the smooth, acoustically filtered low-Mach regime the two coincide; the conservative form earns its keep as gradients sharpen.

The conserved variable and conservation by construction

The DDF carries the volumetric energy density as the zeroth moment of a second distribution \(g_i\),

(6)\[ \sum_i g_i = E \equiv \rho\,h = \rho\,c_p\,(T - T_{\text{ref}}) . \]

Lattice-Boltzmann collision conserves the zeroth moment locally and streaming only relabels populations, so the domain sum \(\sum_\Omega E\) changes only through boundary flux and the integrated source \(Q\). Energy is conserved to machine precision by construction - the property the FD upwind advection structurally lacks, and the primary correctness oracle for this backend. The conserved primary variable is the energy, not the temperature.

Per node the backend stores only the conserved zeroth moment \(E\) and the rank-1 non-equilibrium flux \(q_\alpha^{neq}\), reconstructing the full \(g_i\) in registers at collision rather than holding the population set in memory - the same macroscopic-only strategy the fluid uses. Carrying the thermal flux as a first-order moment and storing macroscopic variables plus that flux in place of the populations is an established low-memory DDF route that recovers the energy equation at roughly half the memory of a population-storing DDF Mochizuki and Tagawa[2].

The datum couples the energy sum to mass

Because the datum sets \(h=0\) at \(T_{\text{ref}}\), the density \(E = \rho\,c_p\,(T - T_{\text{ref}})\) can be negative, and “conservation of \(\sum E\)” is the conservation of a datum-shifted quantity that involves \(T_{\text{ref}}\sum\rho\). Under variable density with open boundaries this couples the energy bookkeeping to the mass flux. The bitwise conservation oracle is therefore stated on a closed or periodic domain, where there is no boundary mass or energy flux and the statement is unambiguous; open-domain runs are checked against the energy-budget closure (advective plus diffusive boundary flux plus integrated source) rather than against bitwise invariance.

The equilibrium: an advection-diffusion distribution, not a second fluid

The energy density is transported by the prescribed fluid velocity \(u_\alpha\); \(g_i\) is not a second self-advecting fluid. The correct equilibrium is therefore the advection-diffusion (linear, first-order in \(u\)) form He et al.[3], Krüger et al.[4],

(7)\[ g_i^{\text{eq}} = w_i\,E\left(1 + \frac{c_{i\alpha}u_\alpha}{c_{s,g}^2}\right), \]

whose velocity moments are

(8)\[ \sum_i g_i^{\text{eq}} = E, \qquad \sum_i g_i^{\text{eq}} c_{i\alpha} = E\,u_\alpha, \qquad \sum_i g_i^{\text{eq}} c_{i\alpha}c_{i\beta} = E\,c_{s,g}^2\,\delta_{\alpha\beta} . \]

The first moment \(E\,u_\alpha\) reproduces the conservative advective flux \(\partial_\alpha(\rho h\,u_\alpha)\) of (5). The lattice sound speed \(c_{s,g}^2\) of the energy distribution is a free scaling that sets the stable diffusivity range; the default reuses the fluid stencil, \(c_{s,g}^2 = 1/3\).

Considered and rejected: the fluid-type second moment

A natural-looking but incorrect choice is the full hydrodynamic second moment \(\sum_i g_i^{\text{eq}}c_{i\alpha}c_{i\beta} = E\,u_\alpha u_\beta + E\,c_{s,g}^2\delta_{\alpha\beta}\), copied from the fluid equilibrium. That moment belongs to a distribution that advects itself; for a scalar transported by a given \(u\) it injects a spurious \(E\,u_\alpha u_\beta\) stress that is precisely the source of the deviation term in the next section. The energy DDF must use the advection-diffusion moment (8), with no \(E\,u_\alpha u_\beta\) term.

The deviation term: the term that must not be dropped

Carrying the Chapman-Enskog expansion of the BGK advection-diffusion DDF to second order does not give a clean diffusion equation. The diffusive flux picks up a deviation term from the time derivative of the first-order non-equilibrium,

(9)\[ \partial_t E + \partial_\alpha(E\,u_\alpha) = \partial_\alpha\!\left(D_g\,\partial_\alpha E\right) \;-\;\partial_\alpha\!\left[\left(\tau_g - \tfrac12\right)\partial_t(E\,u_\alpha)\right] \;+\;\mathcal{O}(u^2), \qquad D_g = c_{s,g}^2\left(\tau_g - \tfrac12\right) . \]

The extra flux \(-\partial_\alpha[(\tau_g-\tfrac12)\,\partial_t(E\,u_\alpha)]\) is the well-known spurious term of advection-diffusion LBM Krüger et al.[4]. It is \(\mathcal{O}(\mathrm{Ma}\,\partial_t)\) and is not negligible by fiat. The backend handles it explicitly, by the default route of adding the standard correction source that cancels it,

(10)\[ S_i^{\text{corr}} = w_i\,\frac{c_{i\alpha}}{c_{s,g}^2}\,\partial_t(E\,u_\alpha), \]

so that the recovered equation is the exact \(\partial_t E + \partial_\alpha(E u_\alpha) = \partial_\alpha(D_g\,\partial_\alpha E)\) with no regime restriction.

A common pitfall: the deviation term must not be dropped silently

A quick derivation drops the deviation term silently. The correct closure cancels it with the source (10), so the recovered equation matches the target (5) rather than carrying a spurious \(\mathcal{O}(\mathrm{Ma}\,\partial_t)\) flux.

Conduction under variable density: the explicit correction

With the deviation term handled, the natural diffusion of the zeroth moment is \(\partial_\alpha(D_g\,\partial_\alpha E)\). But the physics (5) needs \(\partial_\alpha(D\,\partial_\alpha h)\) with \(D = \mu/\mathrm{Pr}\), a diffusion of the intensive enthalpy \(h\), not of the density \(E = \rho h\). Under variable density these differ, because

(11)\[ \partial_\alpha E = \partial_\alpha(\rho h) = \rho\,\partial_\alpha h + h\,\partial_\alpha\rho , \]

so diffusing \(E\) carries a spurious \(h\,\partial_\alpha\rho\) contribution. Setting the per-node relaxation so that the energy diffusivity equals the physical thermal diffusivity,

(12)\[ D_g \equiv \alpha_{\text{eff}} = \frac{\mu(T)}{\rho\,\mathrm{Pr}} \quad\Longleftrightarrow\quad \tau_g = \frac12 + \frac{\alpha_{\text{eff}}}{c_{s,g}^2}, \]

and adding a correction term that implements \(-\partial_\alpha(\alpha_{\text{eff}}\,h\,\partial_\alpha\rho)\), the net diffusion collapses exactly onto the target:

(13)\[ \partial_\alpha\!\left(\alpha_{\text{eff}}\,\partial_\alpha E\right) - \partial_\alpha\!\left(\alpha_{\text{eff}}\,h\,\partial_\alpha\rho\right) = \partial_\alpha\!\left(\alpha_{\text{eff}}\left[\partial_\alpha(\rho h) - h\,\partial_\alpha\rho\right]\right) = \partial_\alpha\!\left(\alpha_{\text{eff}}\,\rho\,\partial_\alpha h\right) = \partial_\alpha\!\left(\frac{\mu}{\mathrm{Pr}}\,\partial_\alpha h\right) = \partial_\alpha\!\left(\lambda\,\partial_\alpha T\right). \]

Equivalently, the correction turns “diffuse the conserved density \(E\)” into “diffuse the intensive enthalpy \(h\) with coefficient \(\mu/\mathrm{Pr}\)” - which is precisely the conduction operator (2) the FD backend already uses. The DDF conduction and the FD conduction are the same operator; the DDF differs only by carrying it conservatively. The relaxation \(\tau_g\) is built per node from the local equation-of-state density and the viscosity law \(\mu(T)\), the identical machinery to the FD backend’s per-node \(\alpha_{\text{eff}}\) (Eq. (10)) and to the fluid’s per-node LES relaxation.

The correction is \(\mathcal{O}(\partial_\alpha\rho)\): it vanishes at constant density and is non-zero only in the variable-density regime. This is the variable-density term that must never be silently dropped, and the DDF makes it an explicit, named flux rather than an implicit approximation.

Implementation constraint: the correction must be a conservative face flux

The correction \(-\partial_\alpha(\alpha_{\text{eff}}\,h\,\partial_\alpha\rho)\) must be discretised in conservative face-flux (divergence) form, the same way the FD conduction is (Eq. (11)). If it is added instead as a node-local collision source, its discrete sum over the domain is not identically zero at finite stencil, so it acts as a spurious energy source or sink wherever \(\partial_\alpha\rho\) is large - at density gradients and, worse, at multiblock coarse/fine interfaces - which would break the machine-precision conservation oracle exactly where it matters most. Conservative face-flux form is a design constraint, not an implementation detail.

Reductions: the two limits the backend must collapse onto

Two limits check the closure against code that is already trusted.

• Constant density (\(\partial_\alpha\rho = 0\)): the variable-density correction vanishes, and \(\partial_\alpha(\alpha_{\text{eff}}\,\partial_\alpha E) = \partial_\alpha((\mu/(\rho_0\mathrm{Pr}))\,\rho_0\,\partial_\alpha h) = \partial_\alpha((\mu/\mathrm{Pr})\,\partial_\alpha h)\), the FD backend’s conduction (2).

• Constant density and divergence-free flow (\(\partial_\alpha u_\alpha = 0\)): (5) becomes \(\rho_0(\partial_t h + u_\alpha\partial_\alpha h) = \partial_\alpha(D\,\partial_\alpha h)\), the plain advection-diffusion equation Eq. (3) that the scalar-transport thermal route already integrates, with diffusivity \(\kappa = \nu/\mathrm{Pr}\).

So the DDF energy reduces, in form, to both the FD conduction and the Boussinesq scalar route in the appropriate limits. These are reduction oracles, distinct from the conservation and diffusion-fidelity oracles.

Why energy needs its own distribution: the stencil-order argument

A second-order polynomial equilibrium (\(N=2\)) needs the even-parity lattice velocity tensors isotropic up to rank \(2N = 4\) Philippi et al.[5]. The advection-diffusion energy equilibrium (7) is \(N=2\), so it needs rank-4 isotropy, which both D3Q19 and D3Q27 provide (verified in tests/mod/t_LBM/testLatticeMomentIsotropy.py: the rank-2 stress and the rank-4 lattice tensor are exact on both lattices). The minimal D3Q7 scalar stencil supports a rank-2 moment but its rank-4 tensor is not isotropic, so it is adequate for pure diffusion but not for an advection-diffusion energy with a correct flux.

This is why energy is a DDF, not a higher moment of the fluid distribution

A full compressible thermal equilibrium, carrying energy as a higher moment of the single fluid distribution \(f\), is a fourth-order equilibrium (\(N=4\), the energy-flux and hyperflux moments) and needs rank-8 isotropy Philippi et al.[5]. The production D3Q19 and D3Q27 stencils do not provide rank-8 isotropy. That is the reason energy is carried on its own second distribution \(g_i\) - a double-distribution function - rather than as a higher moment of \(f\): the energy-only distribution needs only the rank-4 isotropy the production stencils already have. The default energy stencil reuses the fluid \(c_i\) (D3Q27); D3Q19 is also viable.

The alternative the literature takes for genuinely compressible, high-Mach thermal flow is to keep energy on the second distribution but give it the full high-order equilibrium on an extended stencil: Latt et al.[6], for instance, match the first thirteen Maxwell-Boltzmann moments through a numerically guided equilibrium on a D3Q39 lattice to reach supersonic, polyatomic flow. That is the compressible endgame of the DDF idea; the low-Mach energy backend here deliberately needs none of it, because in the low-Mach, low-Eckert regime the energy reduces to the advection-diffusion transport (5) that a second-order equilibrium on D3Q27 carries exactly.

Terms kept and terms dropped, stated explicitly

For the record, against the full compressible energy equation (2):

• Viscous dissipation \(\Phi\) (\(\sim\mathrm{Ec}\)): dropped, low Eckert number. Same as the FD backend.

• Pressure work \(\mathrm{D}p/\mathrm{D}t\) (\(\sim\mathrm{Ma}^2\)): dropped, low Mach number. Same as the FD backend. The conservative DDF form admits pressure work as a divergence-form term; the genuinely compressible route adds it through the high-order guided equilibrium of Latt et al.[6].

• Dilatation in the advection: kept, in conservative flux form (the continuity identity, (4)).

• The advection-diffusion deviation term (9): kept and cancelled by the correction source (10) - not silently dropped.

• The variable-density conduction correction (13): kept as an explicit conservative face flux.

Coupling to the fluid: unchanged from the FD backend

The DDF changes only how the energy variable is transported. The coupling to the fluid is the same seam the FD backend uses: \(g_i \to h \to T = T_{\text{ref}} + h/c_p \to \rho = P/(rT) \to\) buoyancy \((\rho - \rho_\infty)\,G_\alpha \to\) momentum, with the flow advecting \(g_i\) in return. Nothing in the equation-of-state, reduced-pressure or buoyancy core changes when the transport backend is switched; that is what makes FD and DDF two backends of one energy model rather than two thermal solvers.

Free choices that leave the recovered equation unchanged

Two settings are free to choose without altering the recovered equation (5): the value of the energy lattice sound speed \(c_{s,g}^2\), and whether the energy stencil is D3Q19 or D3Q27. The closure - the conserved variable (6), the advection-diffusion equilibrium (7), the deviation-cancelling correction (10), the per-node relaxation (12), and the conservative variable-density correction (13) - is the same regardless.

---

[1]

Gauthier Wissocq and Pierre Sagaut. Hydrodynamic limits and numerical errors of isothermal lattice boltzmann schemes. Journal of Computational Physics, 450:110858, 2022. doi:10.1016/j.jcp.2021.110858.

[2]

Y. Mochizuki and T. Tagawa. Low-memory-consumption system with the lattice Boltzmann method based on optional moment definitions including gradient information. Physics of Fluids, 37(10):107161, 2025. doi:10.1063/5.0290448.

[3]

Xiaoyi He, Shiyi Chen, and Gary D. Doolen. A novel thermal model for the lattice Boltzmann method in incompressible limit. Journal of Computational Physics, 146(1):282–300, 1998. doi:10.1006/jcph.1998.6057.

[4] (1,2)

Timm Krüger, Halim Kusumaatmaja, Alexandr Kuzmin, Orest Shardt, Goncalo Silva, and Erlend Magnus Viggen. The lattice Boltzmann method. Graduate Texts in Physics. Springer International Publishing, Basel, Switzerland, 1 edition, November 2016.

[5] (1,2)

Paulo C. Philippi, Luiz A. Hegele Jr., Luís O. E. dos Santos, and Rodrigo Surmas. From the continuous to the lattice Boltzmann equation: the discretization problem and thermal models. Physical Review E, 73(5):056702, 2006. doi:10.1103/PhysRevE.73.056702.

[6] (1,2)

Jonas Latt, Christophe Coreixas, Joel Beny, and Andrea Parmigiani. Efficient supersonic flows through high-order guided equilibrium with lattice Boltzmann. arXiv preprint arXiv:1910.13515, 2020. arXiv:1910.13515, doi:10.48550/arXiv.1910.13515.