Incompressible Navier-Stokes-Fourier System

• The incompressible Navier-Stokes–Fourier system is a coupled PDE model describing a divergence-free velocity field, pressure as a Lagrange multiplier, and a thermal variable with Fourier diffusion.
• It encompasses various formulations including deterministic, stochastic, and electromagnetic-coupled models, with kinetic derivations and scaling playing a crucial role.
• The framework supports rigorous convergence analysis, handling of boundary effects, and regime transitions, which are essential for modeling fluid dynamics in diverse geometries.

The incompressible Navier–Stokes–Fourier (NSF) system denotes, in the cited literature, the coupled evolution of a divergence-free velocity field, a pressure, and a scalar thermal variable, with viscous momentum diffusion and Fourier-type heat transport. In its canonical low-Mach kinetic-limit form on $\mathbb{R}^3$, it is written as

$$\partial_t u + (u\cdot\nabla_x)u + \nabla_x p = \nu\,\Delta_x u,\qquad \nabla_x\cdot u = 0,$$

$$\partial_t \theta + u\cdot\nabla_x\theta = \kappa\,\Delta_x \theta,$$

with transport coefficients determined by the underlying collision operator (Jiang et al., 2014). Across current work, however, the same label also covers normalized deterministic formulations with viscous heating, Poisson- or Maxwell-coupled variants, stochastic thermal-noise models, and generalized constitutive laws, so the exact meaning of “incompressible NSF” is scaling- and model-dependent (Gess et al., 27 Mar 2026).

1. Core equations and constitutive structure

A common structural feature is the incompressibility constraint

$\nabla\cdot u = 0,$

which makes the pressure a Lagrange multiplier rather than an independent thermodynamic variable. In kinetic low-Mach derivations, the temperature equation is typically advective–diffusive, with no viscous heating term surviving at leading order under the chosen incompressible scaling (Jiang et al., 2014). By contrast, in the normalized deterministic system on the three-dimensional torus $T^3$ with constant mass density $\rho=1$, shear viscosity $\eta=1/2$ and thermal conductivity $\kappa=1$, the equations are

$$\partial_t u = \nabla\cdot(\nabla_{\mathrm{sym}}u)-\Pi\nabla\cdot(u\otimes u),\qquad \nabla\cdot u=0,$$

$$\partial_t \theta = \Delta \theta - \nabla\cdot(u\theta) + |\nabla_{\mathrm{sym}}u|^2,$$

where $$\partial_t u + (u\cdot\nabla_x)u + \nabla_x p = \nu\,\Delta_x u,\qquad \nabla_x\cdot u = 0,$$0 and $$\partial_t u + (u\cdot\nabla_x)u + \nabla_x p = \nu\,\Delta_x u,\qquad \nabla_x\cdot u = 0,$$1 is the Helmholtz projection (Gess et al., 27 Mar 2026).

The thermal variable itself is not fixed once and for all. Depending on the formulation, it appears as $$\partial_t u + (u\cdot\nabla_x)u + \nabla_x p = \nu\,\Delta_x u,\qquad \nabla_x\cdot u = 0,$$2, as an internal-energy variable, or as a transformed variable such as

$$\partial_t u + (u\cdot\nabla_x)u + \nabla_x p = \nu\,\Delta_x u,\qquad \nabla_x\cdot u = 0,$$3

introduced to restore coercivity and strict convexity in the stochastic thermal-noise setting (Gess et al., 27 Mar 2026). Likewise, the constitutive assumptions range from constant coefficients to temperature-dependent viscosity and conductivity, and further to nonlinear constitutive maps on the symmetric gradient in generalized NSF-like systems (Maringová et al., 2017).

2. Principal deterministic and stochastic variants

The plain incompressible NSF system is only one member of a broader family of incompressible velocity–temperature models. In the one-species Vlasov–Poisson–Boltzmann diffusive limit, the incompressible Navier–Stokes–Fourier–Poisson system takes the form

$$\partial_t u + (u\cdot\nabla_x)u + \nabla_x p = \nu\,\Delta_x u,\qquad \nabla_x\cdot u = 0,$$4

$$\partial_t u + (u\cdot\nabla_x)u + \nabla_x p = \nu\,\Delta_x u,\qquad \nabla_x\cdot u = 0,$$5

$$\partial_t u + (u\cdot\nabla_x)u + \nabla_x p = \nu\,\Delta_x u,\qquad \nabla_x\cdot u = 0,$$6

$$\partial_t u + (u\cdot\nabla_x)u + \nabla_x p = \nu\,\Delta_x u,\qquad \nabla_x\cdot u = 0,$$7

so that electrostatics is encoded by an elliptic Boussinesq-type relation rather than by a direct Poisson force in the final projected momentum equation (Guo et al., 2020). When the electrostatic coupling is removed, this reduces to the standard incompressible NSF system with the same kinetic transport coefficients (Guo et al., 2020).

In the presence of general external forcing on the torus, the limiting incompressible system acquires both a body-force contribution and a Boussinesq relation,

$$\partial_t u + (u\cdot\nabla_x)u + \nabla_x p = \nu\,\Delta_x u,\qquad \nabla_x\cdot u = 0,$$8

$$\partial_t u + (u\cdot\nabla_x)u + \nabla_x p = \nu\,\Delta_x u,\qquad \nabla_x\cdot u = 0,$$9

$$\partial_t \theta + u\cdot\nabla_x\theta = \kappa\,\Delta_x \theta,$$0

with $$\partial_t \theta + u\cdot\nabla_x\theta = \kappa\,\Delta_x \theta,$$1 not required to be small, decaying, or of gradient form (Briant et al., 2019).

For two-species kinetic plasma models, the incompressible NSF component is part of a larger incompressible Navier–Stokes–Fourier–Maxwell or Navier–Stokes–Fourier–Poisson system. A representative one-fluid momentum–temperature subsystem is

$$\partial_t \theta + u\cdot\nabla_x\theta = \kappa\,\Delta_x \theta,$$2

$$\partial_t \theta + u\cdot\nabla_x\theta = \kappa\,\Delta_x \theta,$$3

supplemented by Maxwell equations and Ohm’s law

$$\partial_t \theta + u\cdot\nabla_x\theta = \kappa\,\Delta_x \theta,$$4

or closely related closures, depending on the kinetic source model and scaling (Lei et al., 2023).

A qualitatively different variant is the incompressible NSF system with thermal noise derived through the GENERIC framework. In Stratonovich form, the equations acquire nonlinear gradient noise terms

$$\partial_t \theta + u\cdot\nabla_x\theta = \kappa\,\Delta_x \theta,$$5

which couple velocity and temperature even though the deterministic incompressible equations decouple at the momentum level (Gess et al., 27 Mar 2026). The transformed variable $$\partial_t \theta + u\cdot\nabla_x\theta = \kappa\,\Delta_x \theta,$$6 produces a strictly convex energy

$$\partial_t \theta + u\cdot\nabla_x\theta = \kappa\,\Delta_x \theta,$$7

and a positive Itô-sign structure in the gradient nonlinearity (Gess et al., 27 Mar 2026).

3. Kinetic origin and transport coefficients

A major part of the modern theory treats incompressible NSF as a hydrodynamic limit of kinetic equations. In the classical Boltzmann setting, one expands around the global Maxwellian

$$\partial_t \theta + u\cdot\nabla_x\theta = \kappa\,\Delta_x \theta,$$8

using either

$$\partial_t \theta + u\cdot\nabla_x\theta = \kappa\,\Delta_x \theta,$$9

or

$\nabla\cdot u = 0,$0

depending on the normalization (Jiang et al., 2014). The linearized collision operator has null space generated by collision invariants, and the macro–micro decomposition separates the hydrodynamic part from the dissipative kinetic remainder. In the classical Boltzmann case, the null space is

$\nabla\cdot u = 0,$1

or, equivalently under a different normalization,

$\nabla\cdot u = 0,$2

and the incompressible limit resides precisely in this macroscopic subspace (Jiang et al., 2014).

The transport coefficients are given by linear response problems for the linearized operator. In the VPB-to-NSFP derivation, one defines

$\nabla\cdot u = 0,$3

solves

$\nabla\cdot u = 0,$4

and obtains

$\nabla\cdot u = 0,$5

These coefficients are strictly positive and depend on the collision kernel (Guo et al., 2020). In the classical Boltzmann derivation on $\nabla\cdot u = 0,$6, the same mechanism appears through auxiliary problems for the traceless stress source

$\nabla\cdot u = 0,$7

and the heat-flux source

$\nabla\cdot u = 0,$8

again leading to kinetic formulas for $\nabla\cdot u = 0,$9 and $T^3$0 (Jiang et al., 2014).

Quantum statistics modifies both the linearized operator and the nonlinear structure. For the Boltzmann–Fermi–Dirac equation, the equilibrium weight becomes $T^3$1, the null space remains generated by $T^3$2, $T^3$3, and $T^3$4, but the Pauli factors generate a trilinear operator $T^3$5 absent in the classical Maxwellian case. The limiting incompressible NSF–Fourier system is still

$T^3$6

with $T^3$7 and $T^3$8 characterized by the Fermi–Dirac-weighted equations

$T^3$9

and the associated weighted inner products (Jiang et al., 2021).

Charged multispecies models further enlarge the macroscopic closure. For two-species VPB, VMB, or VML systems, the limit carries incompressible velocity and temperature together with charge density, current, electric field, magnetic field, and Ohm’s law. In these settings, conductivity $\rho=1$0 is again a kinetic transport coefficient, defined through an auxiliary equation for the charge mode, and the incompressible NSF subsystem becomes one block in a larger NSFP or NSFM structure (Fang et al., 2021).

4. Solution theories and convergence frameworks

On the deterministic kinetic side, a recurring theorem pattern is uniform-in-$\rho=1$1 global existence near equilibrium, followed by rigorous convergence of moments to the incompressible NSF or NSFP equations. For the scaled Boltzmann equation with general collision kernels, global classical solutions $\rho=1$2 are obtained for small initial data, with an $\rho=1$3-uniform estimate

$\rho=1$4

and the limit $\rho=1$5 yields a global classical solution of incompressible NSF in the whole space (Jiang et al., 2014). The one-species VPB analysis proves the analogous uniform estimate

$\rho=1$6

and then identifies the NSFP limit with strong convergence of $\rho=1$7 and $\rho=1$8 in $\rho=1$9 (Guo et al., 2020).

The stochastic thermal-noise model on $\eta=1/2$0 has a different solution theory. There the $\eta=1/2$1 transformation leads to global-in-time weak martingale solutions for $\eta=1/2$2 initial data, local-in-time strong solutions for regular strictly positive data, and weak–strong uniqueness. The temperature equation is relaxed by a nonnegative random measure $\eta=1/2$3 satisfying

$\eta=1/2$4

and the global admissibility condition is the pathwise energy inequality

$\eta=1/2$5

A relative-energy method yields weak–strong uniqueness, with all stochastic integral contributions canceling in the Itô expansion of the relative energy (Gess et al., 27 Mar 2026).

For generalized incompressible NSF systems with temperature-dependent non-Newtonian stresses and inhomogeneous Dirichlet temperature data, a different analytical issue is central: entropy equality at the weak-solution level. In that setting one has

$\eta=1/2$6

and the entropy variable $\eta=1/2$7 satisfies the weak equality

$\eta=1/2$8

not merely an inequality (Abbatiello et al., 2022). This is established by Galerkin approximation, monotonicity identification of the stress, and a delicate strong-convergence argument for the entropy production term (Abbatiello et al., 2022).

5. Domain, boundary conditions, and singular-limit effects

The incompressible NSF system is highly sensitive to geometry and boundary scaling. On $\eta=1/2$9, kinetic derivations avoid boundary layers and produce whole-space incompressible NSF or NSFP limits (Jiang et al., 2014). On $\kappa=1$0, periodicity eliminates boundary terms and simplifies Helmholtz projection, Korn identities, and Fourier-based compactness arguments (Gess et al., 27 Mar 2026). In bounded domains, by contrast, the boundary layer structure becomes part of the hydrodynamic limit itself.

For the Boltzmann equation in a bounded $\kappa=1$1 domain with Maxwell reflection, the incompressible NSF limit depends on the ratio between the accommodation coefficient $\kappa=1$2 and the Knudsen number $\kappa=1$3. If $\kappa=1$4, the macroscopic limit has Dirichlet boundary conditions

$\kappa=1$5

and the convergence is strong. If $\kappa=1$6, the limit boundary condition is Navier-slip for velocity and Robin for temperature. The decisive new mechanism is immediate damping of acoustic waves in a boundary layer in time, caused by coupled viscous and kinetic boundary layers in space (Jiang et al., 2015).

External forcing changes not only the macroscopic momentum equation but also the kinetic construction used to justify the limit. On the torus, a general non-potential force can be handled by linearizing around a time-dependent Maxwellian rather than a static global equilibrium. The resulting incompressible limit contains the force $\kappa=1$7, while the kinetic theory remains uniform in $\kappa=1$8 on arbitrary finite time intervals (Briant et al., 2019).

Free boundaries introduce yet another layer of structure. A recent generalized incompressible NSF traveling-wave theory treats a viscous, incompressible, heat-conducting fluid in a horizontally infinite strip with flat rigid bottom and moving upper surface, with gravity, external body force, external boundary stress, external heat source, temperature-dependent viscosity, and temperature-dependent capillary coefficients producing Marangoni stresses. In the moving frame, the unknowns solve a steady free-boundary system with kinematic, dynamic, and thermal boundary conditions, and a small-data Sobolev well-posedness theory is proved (Choi et al., 22 Mar 2026).

At the opposite asymptotic extreme, the full compressible Navier–Stokes–Fourier system with small Mach number and large Reynolds and Peclet numbers converges, on $\kappa=1$9, to an incompressible Euler–Boussinesq-type limit rather than to incompressible NSF. This shows that incompressible NSF occupies a specific place inside a broader singular-limit hierarchy governed by the relative scaling of Mach, Reynolds, and Peclet numbers (Feireisl et al., 2012).

6. Generalizations, regime transitions, and conceptual clarifications

Several recent directions extend the incompressible NSF framework beyond Newtonian constant-coefficient models. One class replaces the linear stress law by a temperature-parametrized maximal monotone $$\partial_t u = \nabla\cdot(\nabla_{\mathrm{sym}}u)-\Pi\nabla\cdot(u\otimes u),\qquad \nabla\cdot u=0,$$0-graph: $$\partial_t u = \nabla\cdot(\nabla_{\mathrm{sym}}u)-\Pi\nabla\cdot(u\otimes u),\qquad \nabla\cdot u=0,$$1 with $$\partial_t u = \nabla\cdot(\nabla_{\mathrm{sym}}u)-\Pi\nabla\cdot(u\otimes u),\qquad \nabla\cdot u=0,$$2. This admits temperature-driven transitions between Navier–Stokes, Bingham, and Euler-with-activation regimes, together with a boundary threshold-slip graph depending continuously on temperature (Maringová et al., 2017). The corresponding weak-solution theory is large-data and long-time, but it belongs to an NSF-like generalized constitutive class rather than to the canonical linear incompressible NSF model (Maringová et al., 2017).

Another direction uses specially designed microstructure terms in a Boltzmann model with periodic velocity to derive an incompressible NSF-type system with added forcing and “rubbish” terms. In that setting, a global-in-time weak solution with bounded enstrophy is obtained for arbitrary initial data, but the limiting equations are explicitly not the classical NSF system; the result depends on periodic velocity and on a nonstandard collision operator with a cubic cancellation term (Gu et al., 2023).

Two recurrent misconceptions are therefore corrected by the cited work. First, incompressible NSF is not a single universal PDE. Under the diffusive kinetic scalings of the Boltzmann and VPB derivations, the temperature equation is purely advective–diffusive, whereas in the normalized deterministic and stochastic torus model it contains viscous heating and, under GENERIC, nonlinear gradient noise (Jiang et al., 2014). Second, additional electrostatic or electromagnetic physics does not merely append a lower-order forcing term; it can change the closure itself through elliptic Boussinesq-type constraints, Maxwell evolution, and Ohm’s law (Guo et al., 2020).

A further plausible implication, explicitly noted in the strict-temperature-positivity study for compressible NSF, is that comparison-principle arguments based on the transform

$$\partial_t u = \nabla\cdot(\nabla_{\mathrm{sym}}u)-\Pi\nabla\cdot(u\otimes u),\qquad \nabla\cdot u=0,$$3

should be simpler in the incompressible case because the pressure–dilatation term is absent. That paper states that the positivity mechanism carries over to incompressible NSF with inhomogeneous Dirichlet temperature, and may even permit time-independent lower bounds in idealized settings (Feireisl, 2023).

Taken together, these results present the incompressible Navier–Stokes–Fourier system as both a concrete PDE model and a limit class. Its exact form is fixed not only by incompressibility and Fourier diffusion, but also by the kinetic scaling, constitutive law, stochastic structure, electromagnetic or electrostatic coupling, and domain geometry under consideration.