A proposed numerical method for matched-density Cahn–Hilliard–Navier–Stokes flows comes with an analytical guarantee of unconditional long-time stability under the theorem’s stated conditions. But one of the study’s own cavity tests supplied a warning: with theta = 0.51 and epsilon = 0, the computation eventually lost boundedness at all three displayed time-step sizes. Making the time step smaller did not restore robustness.
The scheme is a second-order linear IMEX method. It extrapolates nonlinear terms, uses an auxiliary-variable formulation for free energy, and adds temporal-curvature regularization controlled by epsilon. Each time step requires only linear solves.
A guarantee with a clear boundary
The stability result comes from a discrete energy estimate and is stated for theta in (1/2, 1] and epsilon >= 0. The paper’s unconditional claim is about long-time stability across that parameter range for the reported discretization.
That guarantee has an important boundary. The theorem assumes homogeneous velocity boundary conditions, while the moving lid in the cavity benchmark violates that assumption. Physical energy in that driven case is therefore not expected to decrease monotonically.
The temporal-accuracy study used two mixed finite-element configurations: P2 − P1 − P2 − P2 with degree k = 2, and P3 − P2 − P3 − P3 with degree k = 3. At Δt = 1.25 × 10−2, the P2 configuration produced reported convergence rates of 2.00 for velocity, 1.96 for pressure and 2.01 for the phase field. Those values are approximately second-order in the tested configuration.
The benchmarks tell a more mixed story
The principal cavity Cauchy tests varied theta across 0.51, 0.75 and 1.0, and epsilon across 0, 0.5ν, ν and 2ν. For most of those configurations, the measured interior Cauchy rates for velocity, phase and pressure were close to second order. Exploratory cases with epsilon = 5ν and 10ν instead produced irregular apparent rates that exceeded 3 or even 10; the study did not treat those numbers as genuine higher-order convergence.
The cavity sensitivity tests also separated moderate regularization from the unregularized case. With epsilon = 0.5ν and epsilon = ν, energy and enstrophy remained bounded over the tested interval. The numerical results also indicate that theta = 1 gives more robust behavior when epsilon = 0, while moderate positive epsilon can improve robustness near the weakly dissipative Crank–Nicolson limit.
In the spinodal decomposition calculation, phase mass stayed constant to numerical precision through T = 100. During the same run, energy continued to dissipate and the domains coarsened.
A separate square-droplet benchmark showed the interface approaching an approximately circular equilibrium by t = 1.0. The calculation remained stable and preserved the enclosed phase mass to numerical precision.
In a two-phase lid-driven cavity benchmark, the method remained stable while resolving a strongly deformed diffuse interface, with no visible spurious oscillations at the plotted resolution. That visual result is narrower than a general stability claim: the moving-lid boundary is outside the homogeneous-velocity condition used in the theorem.
The Rayleigh–Taylor benchmark added another interfacial-flow test. In both viscosity cases, the discretization remained numerically stable over the simulated interval and captured a qualitative shift from a strongly damped evolution to more intricate low-viscosity interfacial flow.
A candidate, not a verdict
Taken together, the calculations present the proposal as a second-order candidate with bounded behavior in the reported regularized and benchmark runs. The two clearest qualifications are that the analytical stability result does not cover the moving-lid boundary, and very large epsilon values generated irregular rates that should not be read as genuine higher-order accuracy.
The document is a preprint identified as arXiv:2608.26046v1 and dated 26 Aug 2026.
Paper data and sources
Original title: A family of second order, linear, unconditionally stable methods for the Cahn-Hilliard-Navier-Stokes equations
Authors: Daozhi Han, Nan Jiang, Jonah H. Nissan, Sayantan Sarkar
Journal/Repository: arXiv
Status: Preprint, not yet peer-reviewed
First online: 2026-08-26
DOI: Not available
Original paper · Full text