Inhomogeneous Navier-Stokes Equations

• Inhomogeneous Navier-Stokes Equations are models for incompressible flows where variable density is actively transported and coupled with momentum dynamics.
• They employ advanced analytical frameworks like Besov and Sobolev spaces to address variable coefficients and achieve global well-posedness in diverse settings.
• The theory extends to free-boundary, stationary, and computational models, providing insights into interface regularity and boundary effects.

The inhomogeneous Navier–Stokes equations govern incompressible viscous flows with non-constant density, and in many formulations also with density-dependent viscosity. In their standard incompressible form they couple transport of density to a variable-coefficient momentum balance, so the system lies between the homogeneous Navier–Stokes equations and more general compressible models: incompressibility imposes $\operatorname{div}u=0$, but the density remains an active unknown transported by the flow. Modern analysis treats this class across whole-space, half-space, bounded-domain, free-boundary, perforated-domain, stationary, and numerical settings, with critical Besov, Sobolev, multiplier, maximal-regularity, and Lagrangian frameworks all playing central roles (Danchin et al., 2013, Danchin et al., 2011, Huang et al., 2024, Mucha et al., 8 Dec 2025).

1. Governing equations and model variants

A standard incompressible inhomogeneous Navier–Stokes system on a domain $\Omega\subset\mathbb R^d$ is

$$\begin{cases} \partial_t \rho + \operatorname{div}(\rho u)=0,\ \partial_t(\rho u)+\operatorname{div}(\rho u\otimes u)-\mu\Delta u+\nabla \Pi=0,\ \operatorname{div}u=0, \end{cases}$$

with either no-slip, stress, or free-boundary conditions depending on the geometry. Because $\operatorname{div}u=0$, the continuity equation is equivalently

$\partial_t \rho + u\cdot\nabla \rho =0,$

so density is advected by a volume-preserving flow rather than diffused or compressed (Danchin et al., 2013, Danchin et al., 2016).

A frequently used reformulation assumes the density stays away from zero and introduces

$a=\frac1\rho-1.$

Then $\rho=(1+a)^{-1}$, and the system becomes

$$\begin{cases} \partial_t a + u\cdot\nabla a =0,\ \partial_t u + u\cdot\nabla u + (1+a)(\nabla\Pi-\mu\Delta u)=0,\ \operatorname{div}u=0. \end{cases}$$

This form isolates the pure transport structure of the density fluctuation and makes the inhomogeneity appear as a variable coefficient multiplying the Stokes operator (Danchin et al., 2013).

A broader class replaces constant viscosity by $\mu(\rho)$. In that setting the momentum equation is written as

$$(\rho u)_t+\operatorname{div}(\rho u\otimes u)+\nabla P-\operatorname{div}\bigl(2\mu(\rho)d\bigr)=0, \qquad d=\frac12(\nabla u+(\nabla u)^\top),$$

which substantially complicates elliptic estimates because the viscous operator itself becomes density-dependent (Huang et al., 2015, Huang et al., 2024, Huang et al., 2012).

Free-boundary formulations add geometric evolution. For a time-dependent domain $\Omega\subset\mathbb R^d$0, the Eulerian system is supplemented by the stress condition

$\Omega\subset\mathbb R^d$1

and the kinematic condition

$\Omega\subset\mathbb R^d$2

so the interface is transported by the fluid itself (Mucha et al., 8 Dec 2025).

Stationary variants replace time derivatives by zero. In three dimensions, a stationary incompressible inhomogeneous system in the whole space is

$\Omega\subset\mathbb R^d$3

while in two dimensions, density-dependent viscosity admits a stream-function reduction to a fourth-order nonlinear elliptic equation (Ding et al., 7 Jan 2025, He et al., 2020).

A persistent misconception is that “inhomogeneous” should be identified with “compressible.” The cited literature treats the opposite regime: the flow is incompressible, but density is variable and transported. The analytic difficulty therefore comes from variable coefficients and transport, not from acoustic modes or volumetric compression.

2. Critical scaling and functional settings

The basic parabolic scaling is

$\Omega\subset\mathbb R^d$4

with initial data

$\Omega\subset\mathbb R^d$5

Under this scaling, density has order zero, while velocity has order one. Consequently, scale-invariant spaces for the velocity include homogeneous Besov classes $\Omega\subset\mathbb R^d$6, and for the density fluctuation one expects either $\Omega\subset\mathbb R^d$7, critical Besov spaces, or multiplier spaces acting on the velocity class (Danchin et al., 2013, Danchin et al., 2016).

Several critical frameworks coexist in the literature. In the whole-space Lagrangian theory, the velocity is taken in $\Omega\subset\mathbb R^d$8, while the density belongs to the multiplier space

$\Omega\subset\mathbb R^d$9

so that multiplication by $$\begin{cases} \partial_t \rho + \operatorname{div}(\rho u)=0,\ \partial_t(\rho u)+\operatorname{div}(\rho u\otimes u)-\mu\Delta u+\nabla \Pi=0,\ \operatorname{div}u=0, \end{cases}$$0 preserves the velocity regularity required by the momentum equation (Danchin et al., 2011). In the half-space theory with merely bounded density, the critical pair is

$$\begin{cases} \partial_t \rho + \operatorname{div}(\rho u)=0,\ \partial_t(\rho u)+\operatorname{div}(\rho u\otimes u)-\mu\Delta u+\nabla \Pi=0,\ \operatorname{div}u=0, \end{cases}$$1

again consistent with scaling (Danchin et al., 2013). In the three-dimensional strong-solution theory with vacuum, the critical smallness condition is imposed instead on the Sobolev norm

$$\begin{cases} \partial_t \rho + \operatorname{div}(\rho u)=0,\ \partial_t(\rho u)+\operatorname{div}(\rho u\otimes u)-\mu\Delta u+\nabla \Pi=0,\ \operatorname{div}u=0, \end{cases}$$2

which is scale invariant in dimension three (Craig et al., 2013).

The rough-density literature makes essential use of multiplier spaces. A key fact is that for a uniformly $$\begin{cases} \partial_t \rho + \operatorname{div}(\rho u)=0,\ \partial_t(\rho u)+\operatorname{div}(\rho u\otimes u)-\mu\Delta u+\nabla \Pi=0,\ \operatorname{div}u=0, \end{cases}$$3-bounded domain $$\begin{cases} \partial_t \rho + \operatorname{div}(\rho u)=0,\ \partial_t(\rho u)+\operatorname{div}(\rho u\otimes u)-\mu\Delta u+\nabla \Pi=0,\ \operatorname{div}u=0, \end{cases}$$4, the characteristic function $$\begin{cases} \partial_t \rho + \operatorname{div}(\rho u)=0,\ \partial_t(\rho u)+\operatorname{div}(\rho u\otimes u)-\mu\Delta u+\nabla \Pi=0,\ \operatorname{div}u=0, \end{cases}$$5 belongs to $$\begin{cases} \partial_t \rho + \operatorname{div}(\rho u)=0,\ \partial_t(\rho u)+\operatorname{div}(\rho u\otimes u)-\mu\Delta u+\nabla \Pi=0,\ \operatorname{div}u=0, \end{cases}$$6 whenever

$$\begin{cases} \partial_t \rho + \operatorname{div}(\rho u)=0,\ \partial_t(\rho u)+\operatorname{div}(\rho u\otimes u)-\mu\Delta u+\nabla \Pi=0,\ \operatorname{div}u=0, \end{cases}$$7

which is precisely why piecewise constant densities can be admissible in critical theories (Danchin et al., 2016). This multiplier perspective is one of the main mechanisms by which density patches enter incompressible variable-density analysis (Danchin et al., 2011).

Anisotropic criticality appears when horizontal and vertical directions are treated differently. For the three-dimensional anisotropic theory, the spaces

$$\begin{cases} \partial_t \rho + \operatorname{div}(\rho u)=0,\ \partial_t(\rho u)+\operatorname{div}(\rho u\otimes u)-\mu\Delta u+\nabla \Pi=0,\ \operatorname{div}u=0, \end{cases}$$8

are built from horizontal and vertical Littlewood–Paley decompositions, and shorthand classes such as $$\begin{cases} \partial_t \rho + \operatorname{div}(\rho u)=0,\ \partial_t(\rho u)+\operatorname{div}(\rho u\otimes u)-\mu\Delta u+\nabla \Pi=0,\ \operatorname{div}u=0, \end{cases}$$9 and $\operatorname{div}u=0$0 encode scale-invariant control of $\operatorname{div}u=0$1 and $\operatorname{div}u=0$2 separately (Chemin et al., 2013). This suggests that criticality in the inhomogeneous problem is not tied to a single canonical Banach space, but to a family of scale-compatible spaces adapted to transport, Stokes regularity, and geometric structure.

3. Well-posedness regimes and global solvability

The contemporary theory splits into several regimes: small-data critical well-posedness, strong solutions with vacuum, anisotropic large-data results under structural smallness, and density-dependent viscosity problems in bounded domains. A common pattern is that density regularity can often be weaker than velocity regularity, because density is transported while velocity is diffused (Danchin et al., 2011, Danchin et al., 2013, Craig et al., 2013, Huang et al., 2015, Huang et al., 2024).

Setting	Hypotheses	Representative conclusion
Whole space, Lagrangian critical theory	$\operatorname{div}u=0$3, $\operatorname{div}u=0$4, small	global unique solution; existence and uniqueness under the same smallness condition (Danchin et al., 2011)
Half-space, bounded density	$\operatorname{div}u=0$5 small, $\operatorname{div}u=0$6 in critical Besov spaces	global weak solution; uniqueness with additional $\operatorname{div}u=0$7 control (Danchin et al., 2013)
3D strong solutions with vacuum	$\operatorname{div}u=0$8, compatibility, $\operatorname{div}u=0$9 small	unique global strong solution; vacuum allowed (Craig et al., 2013)
3D variable viscosity in bounded domains	either $\partial_t \rho + u\cdot\nabla \rho =0,$0 small, or $\partial_t \rho + u\cdot\nabla \rho =0,$1 with $\partial_t \rho + u\cdot\nabla \rho =0,$2 and density sufficiently large	unique global strong solutions in regimes not covered by constant-viscosity critical theory (Huang et al., 2015, Huang et al., 2024)

The whole-space Lagrangian approach is distinctive because it freezes density along trajectories and converts the Eulerian system into a variable-coefficient Stokes problem. The resulting contraction mapping is performed directly in the critical space $\partial_t \rho + u\cdot\nabla \rho =0,$3, which yields Lipschitz dependence on initial data and removes the usual discrepancy between existence and uniqueness assumptions (Danchin et al., 2011).

In the half-space, the central result establishes global existence with density only bounded and close to a positive constant, no derivative regularity of density, and velocity in critical Besov spaces. The smallness condition is anisotropic: it couples $\partial_t \rho + u\cdot\nabla \rho =0,$4 and the horizontal component $\partial_t \rho + u\cdot\nabla \rho =0,$5, while the vertical component $\partial_t \rho + u\cdot\nabla \rho =0,$6 may be large and appears only through an exponential factor. Uniqueness is proved under additional regularity yielding

$\partial_t \rho + u\cdot\nabla \rho =0,$7

for a specific exponent $\partial_t \rho + u\cdot\nabla \rho =0,$8 (Danchin et al., 2013).

The three-dimensional strong theory with vacuum takes a different route. Here density may vanish, but one assumes a compatibility condition and a small critical Sobolev norm $\partial_t \rho + u\cdot\nabla \rho =0,$9. The result is a unique global strong solution together with decay

$a=\frac1\rho-1.$0

in the setting of either $a=\frac1\rho-1.$1 or bounded domains (Craig et al., 2013).

Variable-viscosity problems exhibit two distinct bounded-domain mechanisms. One result proves global-in-time unique strong solutions with arbitrary large initial density and possible vacuum when $a=\frac1\rho-1.$2 is suitably small and $a=\frac1\rho-1.$3 remains positive and bounded (Huang et al., 2015). Another treats

$a=\frac1\rho-1.$4

and shows that sufficiently large lower density stabilizes the three-dimensional system; the abstract explicitly describes this as the first result concerning the existence of large strong solution for the inhomogeneous Navier–Stokes equations in three dimensions (Huang et al., 2024).

Two further large-data mechanisms are structural rather than isotropic. In two dimensions with variable viscosity, global well-posedness is obtained in a critical Besov framework under a non-linear smallness condition requiring the fluctuation of the initial density to be doubly exponential small compared with the size of the initial velocity; in the patch setting this removes the velocity smallness condition from an earlier theorem (Huang et al., 2012). In three dimensions, anisotropic critical Besov theory allows a large vertical initial velocity, provided the density fluctuation and the horizontal component are exponentially small relative to it, and the same anisotropic strategy also yields global solutions for data slowly varying in one direction (Chemin et al., 2013).

4. Rough density, density patches, and interface regularity

One of the defining developments in the subject is the relaxation of density regularity from smooth or Besov classes to merely bounded functions, multipliers, and piecewise constants. In the half-space theory, the density fluctuation $a=\frac1\rho-1.$5 is only assumed to belong to $a=\frac1\rho-1.$6 and to be small, yet the transport structure preserves

$a=\frac1\rho-1.$7

for all time (Danchin et al., 2013). In the whole-space Lagrangian theory, the multiplier space $a=\frac1\rho-1.$8 is broad enough to include small jumps across $a=\frac1\rho-1.$9 interfaces, so piecewise constant densities become admissible data (Danchin et al., 2011).

The archetypal density patch is

$\rho=(1+a)^{-1}$0

with $\rho=(1+a)^{-1}$1 bounded and $\rho=(1+a)^{-1}$2. Since density is transported,

$\rho=(1+a)^{-1}$3

where $\rho=(1+a)^{-1}$4 is the flow map of $\rho=(1+a)^{-1}$5. The nontrivial question is not preservation of the patch form, but preservation of interface regularity (Danchin et al., 2016).

The patch regularity theory answers this using transported vector fields. If $\rho=(1+a)^{-1}$6 is tangent to the initial interface, then the evolved field $\rho=(1+a)^{-1}$7 solves

$\rho=(1+a)^{-1}$8

The analysis of para-vector fields $\rho=(1+a)^{-1}$9, commutators such as $$\begin{cases} \partial_t a + u\cdot\nabla a =0,\ \partial_t u + u\cdot\nabla u + (1+a)(\nabla\Pi-\mu\Delta u)=0,\ \operatorname{div}u=0. \end{cases}$$0, and multiplier estimates for discontinuous densities shows that $$\begin{cases} \partial_t a + u\cdot\nabla a =0,\ \partial_t u + u\cdot\nabla u + (1+a)(\nabla\Pi-\mu\Delta u)=0,\ \operatorname{div}u=0. \end{cases}$$1 regularity of the interface propagates globally in time for small density contrast and small critical velocity, in both two and three dimensions (Danchin et al., 2016). The key point is that regularity is propagated along selected directions, not isotropically; this is the striated-regularity mechanism inherited from vortex-patch theory.

In the multiplier setting, piecewise constant density interacts naturally with critical velocity spaces because characteristic functions of $$\begin{cases} \partial_t a + u\cdot\nabla a =0,\ \partial_t u + u\cdot\nabla u + (1+a)(\nabla\Pi-\mu\Delta u)=0,\ \operatorname{div}u=0. \end{cases}$$2 domains act as multipliers on the relevant Besov classes (Danchin et al., 2011). In the two-dimensional variable-viscosity theory, the same perspective yields a global patch result for

$$\begin{cases} \partial_t a + u\cdot\nabla a =0,\ \partial_t u + u\cdot\nabla u + (1+a)(\nabla\Pi-\mu\Delta u)=0,\ \operatorname{div}u=0. \end{cases}$$3

with $$\begin{cases} \partial_t a + u\cdot\nabla a =0,\ \partial_t u + u\cdot\nabla u + (1+a)(\nabla\Pi-\mu\Delta u)=0,\ \operatorname{div}u=0. \end{cases}$$4 a bounded $$\begin{cases} \partial_t a + u\cdot\nabla a =0,\ \partial_t u + u\cdot\nabla u + (1+a)(\nabla\Pi-\mu\Delta u)=0,\ \operatorname{div}u=0. \end{cases}$$5 domain and $$\begin{cases} \partial_t a + u\cdot\nabla a =0,\ \partial_t u + u\cdot\nabla u + (1+a)(\nabla\Pi-\mu\Delta u)=0,\ \operatorname{div}u=0. \end{cases}$$6 small; the solution preserves patch form

$$\begin{cases} \partial_t a + u\cdot\nabla a =0,\ \partial_t u + u\cdot\nabla u + (1+a)(\nabla\Pi-\mu\Delta u)=0,\ \operatorname{div}u=0. \end{cases}$$7

and the corollary explicitly removes the smallness condition for the initial velocity from the corresponding earlier theorem (Huang et al., 2012).

A common misconception is that interface regularity follows automatically from the transport identity $$\begin{cases} \partial_t a + u\cdot\nabla a =0,\ \partial_t u + u\cdot\nabla u + (1+a)(\nabla\Pi-\mu\Delta u)=0,\ \operatorname{div}u=0. \end{cases}$$8. The patch literature shows that this is false in critical regularity classes: preserving a $$\begin{cases} \partial_t a + u\cdot\nabla a =0,\ \partial_t u + u\cdot\nabla u + (1+a)(\nabla\Pi-\mu\Delta u)=0,\ \operatorname{div}u=0. \end{cases}$$9 interface requires quantitative control of transported tangent fields, commutator structure, and, ultimately, integrability of $\mu(\rho)$0 in $\mu(\rho)$1 (Danchin et al., 2016).

5. Geometry, boundaries, and domain effects

Boundary geometry changes both the linear theory and the admissible smallness mechanisms. In the half-space, explicit Stokes representations based on horizontal/vertical decomposition, Riesz transforms, and reflection operators yield maximal $\mu(\rho)$2 regularity with Dirichlet boundary condition. This makes it possible to prove global small-data existence with bounded density, and it also reveals a genuinely anisotropic feature: the horizontal component of the initial velocity can be estimated from horizontal data alone, which is why the smallness condition distinguishes $\mu(\rho)$3 from $\mu(\rho)$4 (Danchin et al., 2013).

For bounded domains, the same paper gives a partial extension through fractional domains $\mu(\rho)$5 of the Stokes operator. The bounded-domain theorem is isotropic rather than anisotropic, and uniqueness is left open there because the Lagrangian reformulation with nonzero divergence and boundary terms becomes technically harder (Danchin et al., 2013).

Geometry	Main analytical device	Representative outcome
Half-space $\mu(\rho)$6	explicit Stokes formulas and maximal regularity	global small-data existence with $\mu(\rho)$7 (Danchin et al., 2013)
Bounded perforated domains with holes of size $\mu(\rho)$8, $\mu(\rho)$9	zero extension, Bogovskii operator, Aubin–Lions compactness	homogenized limit is the unchanged inhomogeneous Navier–Stokes system (Lu et al., 10 Jan 2025)
Free boundary as perturbation of a half-space	Lagrangian coordinates, maximal $$(\rho u)_t+\operatorname{div}(\rho u\otimes u)+\nabla P-\operatorname{div}\bigl(2\mu(\rho)d\bigr)=0, \qquad d=\frac12(\nabla u+(\nabla u)^\top),$$0 regularity, complex interpolation	global small-data well-posedness in critical $$(\rho u)_t+\operatorname{div}(\rho u\otimes u)+\nabla P-\operatorname{div}\bigl(2\mu(\rho)d\bigr)=0, \qquad d=\frac12(\nabla u+(\nabla u)^\top),$$1-$$(\rho u)_t+\operatorname{div}(\rho u\otimes u)+\nabla P-\operatorname{div}\bigl(2\mu(\rho)d\bigr)=0, \qquad d=\frac12(\nabla u+(\nabla u)^\top),$$2 and Lorentz frameworks (Mucha et al., 8 Dec 2025)

Homogenization in perforated domains provides a different geometric effect. For three-dimensional inhomogeneous incompressible flow in a bounded domain punctured by very tiny holes with diameter $$(\rho u)_t+\operatorname{div}(\rho u\otimes u)+\nabla P-\operatorname{div}\bigl(2\mu(\rho)d\bigr)=0, \qquad d=\frac12(\nabla u+(\nabla u)^\top),$$3 and $$(\rho u)_t+\operatorname{div}(\rho u\otimes u)+\nabla P-\operatorname{div}\bigl(2\mu(\rho)d\bigr)=0, \qquad d=\frac12(\nabla u+(\nabla u)^\top),$$4, weak solutions can be extended by zero to the full domain, and the combination of time-derivative estimates, Leray projection, and Aubin–Lions compactness yields strong convergence of density and momentum. In the subcritical regime of very tiny holes, no Brinkman or Darcy correction appears: the homogenized limit is exactly the original inhomogeneous incompressible Navier–Stokes system on the unperforated domain (Lu et al., 10 Jan 2025).

Free-boundary theory adds a second layer of geometry. A recent framework treats the free surface as a perturbation of the half-space and uses the natural Lagrangian change of variables so that density becomes frozen,

$$(\rho u)_t+\operatorname{div}(\rho u\otimes u)+\nabla P-\operatorname{div}\bigl(2\mu(\rho)d\bigr)=0, \qquad d=\frac12(\nabla u+(\nabla u)^\top),$$5

and the momentum equation differs from the homogeneous one only by the perturbative term $$(\rho u)_t+\operatorname{div}(\rho u\otimes u)+\nabla P-\operatorname{div}\bigl(2\mu(\rho)d\bigr)=0, \qquad d=\frac12(\nabla u+(\nabla u)^\top),$$6. Maximal regularity is then developed in $$(\rho u)_t+\operatorname{div}(\rho u\otimes u)+\nabla P-\operatorname{div}\bigl(2\mu(\rho)d\bigr)=0, \qquad d=\frac12(\nabla u+(\nabla u)^\top),$$7, in time-weighted variants, and in Lorentz classes $$(\rho u)_t+\operatorname{div}(\rho u\otimes u)+\nabla P-\operatorname{div}\bigl(2\mu(\rho)d\bigr)=0, \qquad d=\frac12(\nabla u+(\nabla u)^\top),$$8; complex interpolation is used to control nonlinear boundary terms in fractional Sobolev spaces (Mucha et al., 8 Dec 2025). This suggests a path beyond Besov-multiplier techniques for variable density in moving-boundary problems.

6. Stationary theory, large-time behavior, numerical approximation, and open problems

The time-dependent theory is complemented by precise asymptotic and stationary results. In two dimensions, optimal decay estimates have been established in a Besov framework: $$(\rho u)_t+\operatorname{div}(\rho u\otimes u)+\nabla P-\operatorname{div}\bigl(2\mu(\rho)d\bigr)=0, \qquad d=\frac12(\nabla u+(\nabla u)^\top),$$9 for any $\Omega\subset\mathbb R^d$00 and $\Omega\subset\mathbb R^d$01, provided the initial momentum satisfies

$\Omega\subset\mathbb R^d$02

The result is stated to be optimal even for the classical homogeneous Navier–Stokes equations (Liu, 2024).

At the stationary level, a three-dimensional Liouville-type theorem shows that if a smooth stationary solution satisfies

$\Omega\subset\mathbb R^d$03

together with the low-frequency condition

$\Omega\subset\mathbb R^d$04

then necessarily

$\Omega\subset\mathbb R^d$05

The proof localizes the Dirichlet energy near the origin in frequency space through two successive localizations, first for the interaction between $\Omega\subset\mathbb R^d$06 and $\Omega\subset\mathbb R^d$07, then for the interaction between $\Omega\subset\mathbb R^d$08 and $\Omega\subset\mathbb R^d$09 (Ding et al., 7 Jan 2025). This does not resolve the finite-Dirichlet-energy Liouville problem without extra low-frequency control; rather, it identifies a sufficient rigidity condition in the inhomogeneous setting.

Two-dimensional stationary variable-viscosity theory follows a different route. Using a stream function $\Omega\subset\mathbb R^d$10 with $\Omega\subset\mathbb R^d$11, the system is reduced to a fourth-order nonlinear elliptic equation for $\Omega\subset\mathbb R^d$12. This yields existence of weak solutions, $\Omega\subset\mathbb R^d$13 and higher regularity under smoother data, and explicit parallel, concentric, and radial flows. Piecewise-constant viscosity profiles can be handled explicitly, and the analysis shows that while $\Omega\subset\mathbb R^d$14 may lie in $\Omega\subset\mathbb R^d$15 for all finite $\Omega\subset\mathbb R^d$16, quantities like $\Omega\subset\mathbb R^d$17 or $\Omega\subset\mathbb R^d$18 may fail to belong to $\Omega\subset\mathbb R^d$19 across sharp interfaces (He et al., 2020).

On the computational side, an elementary fully discrete finite-difference method combines a Lax–Friedrichs-type explicit scheme for the density transport equation with a Ladyzhenskaya-type implicit scheme for the momentum equation. Under the assumption that the initial density profile is strictly away from $\Omega\subset\mathbb R^d$20, the scheme is proved to converge strongly, up to a subsequence, to a weak solution on an arbitrary time interval. The analysis introduces a new Aubin–Lions–Simon type compactness argument that interpolates between strong norms of the velocity and a weak norm of the momentum $\Omega\subset\mathbb R^d$21 (Soga, 2023).

Several open directions recur across the literature. The patch-regularity theory explicitly states that extending critical global well-posedness to large density variations with critical velocity fields is totally open (Danchin et al., 2016). The half-space work leaves bounded-domain uniqueness unresolved and identifies less regular boundaries as a serious obstacle because the necessary Stokes maximal-regularity theory becomes delicate or unavailable (Danchin et al., 2013). The large-density stabilization result for $\Omega\subset\mathbb R^d$22 with $\Omega\subset\mathbb R^d$23 is restricted to bounded domains, and the whole-space case remains open in that regime (Huang et al., 2024). More broadly, the coexistence of transport-dominated density dynamics, rough coefficients, and boundary geometry continues to make the inhomogeneous Navier–Stokes equations a testing ground for harmonic analysis, maximal regularity, geometric transport, and compactness methods.