July 10, 2026
The skew-gradient embedding (SGE) framework [1] reformulates a thermodynamically consistent system as a generalized gradient flow by embedding its zero-energy contribution in a skew-symmetric operator. In a time-discrete scheme, the profiles defining this operator may be evaluated at previous time levels. The resulting operator remains skew-symmetric, so its contribution to the discrete energy balance vanishes; this explicit treatment often decouples multiphysics systems. We show that this operator is not unique: the admissible gauges form an affine space, and we call the resulting family generalized skew-gradient embeddings (GSGE). For any positive definite metric, least squares selects a unique minimum-Hilbert–Schmidt gauge, and the native metric recovers SGE. This construction also gives regularized approximations, corrections of non-neutral residuals, and gauges that preserve prescribed invariants. For rank-two gauges, we use a necessary and sufficient Jacobi criterion. Applying this criterion to a compatible MAC discretization of the incompressible Navier–Stokes equations gives a finite-dimensional rank-two Poisson–GENERIC formulation at the semi-discrete level; the implicit midpoint rule preserves this rank-two GENERIC structure at the fully discrete level and satisfies the exact discrete energy law. For the Cahn–Hilliard–Navier–Stokes system, the regularized GSGE–BDF2 scheme preserves mass, dissipates the discrete energy unconditionally, and admits a decoupled implementation.
generalized Onsager principle ,skew-gradient embedding ,gauge freedom ,Poisson structure ,GENERIC ,structure-preserving discretization
Thermodynamically consistent models in classical electrodynamics, fluid and solid mechanics, quantum mechanics, complex fluids, phase-field hydrodynamics, and statistical physics often arise from conservation laws coupled with constitutive relations. Their coupled reversible and irreversible parts admit several structural descriptions, such as the port-Hamiltonian, metriplectic, and GENERIC (General Equation for the Non-Equilibrium Reversible–Irreversible Coupling) formulations [2]–[5]. After a compatible spatial discretization, or at the formal PDE level used to identify the structural identities, such systems can be written in the following generalized form \[\label{eq:intro-onsager-zec} \partial_t\Phi=M(\Phi)G(\Phi)+J(\Phi), \quad G(\Phi):=\nabla F(\Phi).\tag{1}\] Here \(\Phi\) denotes the state variable in a real inner-product space \(H\), \(F(\Phi)\) is the free energy, \(G(\Phi) = \nabla F(\Phi)\) is the thermodynamic force, and \(M(\Phi)\) is a symmetric negative semidefinite mobility operator [6]–[8]. The generalized Onsager principle provides a systematic route to such models and has been applied to complex fluids, liquid crystals, active matter, and interacting particle systems [9]–[13]. The vector field \(J(\Phi)\) satisfying \(\left\langle G(\Phi),J(\Phi)\right\rangle=0\) is the zero-energy contribution (ZEC) [14], [15]; it collects reversible mechanisms such as convection, transport, rotation, and electromagnetic coupling. This orthogonality condition gives the energy dissipation law \[\label{eq:intro-energy-law} \tfrac{\mathrm d}{\mathrm dt}F(\Phi(t)) =\left\langle G(\Phi),\partial_t\Phi\right\rangle =\left\langle G(\Phi),M(\Phi)G(\Phi)\right\rangle \le0 .\tag{2}\]
The skew-gradient embedding (SGE) framework [1] exploits the zero-energy property of \(J\) by embedding it into a skew-symmetric operator. Specifically, 1 is rewritten as \[\label{eq:intro-sge-choice} \partial_t\Phi=\bigl(M(\Phi)+S_*(\Phi)\bigr)G(\Phi), \quad S_*(\Phi)=\tfrac{G(\Phi)\wedge J(\Phi)}{\left\|G(\Phi)\right\|^2}.\tag{3}\] Since \(S_*(\Phi)\) is skew-symmetric, the reversible contribution may be treated explicitly, for instance through frozen profiles from previous time levels, without destroying the discrete energy law. This treatment often decouples multiphysics systems and yields efficient first- and higher-order energy-stable schemes [1].
The rank-two skew operator \(S_*\) in 3 is useful, but it is not the only way to embed the ZEC term. It is enough to find a skew two-form \(\omega\) such that \[\label{eq:intro-gsge} J(\Phi)=(\iota_{G(\Phi)}\,\omega)^\sharp .\tag{4}\] Equivalently, \(\omega(G(\Phi),v)=\left\langle J(\Phi),v\right\rangle\) for every \(v\in H\). Thus the contraction of \(\omega\) with the force equals the covector associated with the reversible field. Since \(\omega(G(\Phi),G(\Phi))=0\), the energy law is unchanged. We call the resulting reformulations generalized skew-gradient embeddings (GSGE). The admissible two-forms form an affine space: the model fixes only the contraction of \(\omega\) along \(G(\Phi)\), while all components annihilated by \(\iota_{G(\Phi)}\) do not affect the energy balance. In this sense, \(S_*\) is one representative of an algebraic gauge freedom. The rest of the paper develops three selection principles for this affine gauge space: least-squares optimality, invariant preservation, and Poisson geometry.
The first principle selects gauges by least-squares optimality. For any symmetric positive definite gauge map \(A\), the operator-weighted two-form \[\label{eq:intro-A-gauge} \omega_A =\tfrac{(A G(\Phi))\wedge J(\Phi)^\flat} {\left\langle A G(\Phi),G(\Phi)\right\rangle}\tag{5}\] is the unique minimum-Hilbert–Schmidt representative in the \(A\)-metric. The original SGE operator corresponds to the native metric. Thus changing \(A\) changes the metric in which the admissible gauge is selected, while preserving the same energy law. The same least-squares principle also gives regularized approximations and projections of non-neutral residuals. It further justifies the regularized differential gauge used later in the Cahn–Hilliard–Navier–Stokes (CHNS) scheme, while the energy estimate remains in the discrete \(L^2\) pairing.
The second principle enforces invariant preservation. If 1 conserves additional quantities, such as mass, momentum, or linear constraints, the gauge can be chosen to annihilate the gradients of these invariants. Then the embedded reversible block makes no contribution to the corresponding discrete balances. An \(A\)-orthogonal projection of the force profile gives an explicit construction, while leaving the remaining gauge freedom available.
The third principle selects gauges with Poisson geometry. In this setting, a rank-two gauge is written as a bivector \(L=Y\wedge J\), acting on covectors by \[L\alpha=\alpha(Y)J-\alpha(J)Y, \quad \alpha\in H' .\] Thus \(L\,\mathrm dF=J\) whenever \(\mathrm dF(Y)=1\) and \(\mathrm dF(J)=0\). The second condition is exactly the ZEC condition. In finite dimensions, this rank-two operator defines a Poisson structure precisely when \([Y,J]\in\operatorname{span}\{Y,J\}\). Under mild nondegeneracy assumptions, such gauges exist locally and can include prescribed invariants as degeneracies of the bracket. They give low-rank, Hamiltonian-specific reversible operators satisfying the GENERIC degeneracy conditions in the corresponding rank-two Poisson sense [16]–[18]. Since the systems considered here are isothermal, this connection is stated in the free-energy convention. For incompressible Navier–Stokes, a formal radial field \(Y=\tfrac{\boldsymbol{u}}{\left\|\boldsymbol{u}\right\|^2}\) motivates the construction, while the rigorous finite-dimensional rank-two Poisson–GENERIC formulation is obtained after a compatible marker-and-cell (MAC) discretization with skew convection.
The main contributions are summarized as follows:
We introduce the GSGE formulation 4 , characterize its affine gauge space, and identify rank-two decomposable gauges and operator-weighted gauges as computationally useful representatives.
We establish a unified least-squares principle for operator-weighted gauges, including regularized gauges, residual-projection corrections, and constructive invariant-preserving gauges.
We use a finite-dimensional necessary and sufficient condition for a rank-two gauge to satisfy the Jacobi identity, prove a local existence result with prescribed degeneracies, and formulate the corresponding GENERIC interpretation in the rank-two Poisson sense.
For the incompressible Navier–Stokes equations, we use a formal rank-two gauge as motivation and then construct a compatible MAC discretization. The resulting semi-discrete system admits a rank-two Poisson–GENERIC formulation distinct from the classical Lie–Poisson structure. The implicit midpoint rule retains this representation at the midpoint and satisfies the discrete energy law.
For the CHNS system, we propose a GSGE–BDF2 scheme based on a three-parameter regularized differential gauge. The scheme preserves mass, dissipates the discrete energy unconditionally, and admits a decoupled implementation. Two parameter limits recover the SGE gauge and a gradient-weighted mass-preserving gauge, respectively.
The remainder of the paper is organized as follows. 2 introduces the finite-dimensional setting, reviews SGE, and formulates GSGE. 3 develops gauge-selection principles based on least-squares optimality, invariant preservation, and Poisson geometry. 4 applies the theory to the Navier–Stokes and CHNS systems. 5 contains concluding remarks.
Let \((H,\left\langle \cdot,\cdot\right\rangle)\) be a finite-dimensional real inner-product space with dual space \(H'\). The structural theory is finite-dimensional. For PDE models, \(H\) denotes the spatially discrete state space, and the continuum formulas in 4 only indicate the identities to be preserved by compatible discretizations. Denote the Riesz isomorphism and its inverse by the musical maps \[\flat:H\to H',\quad v\mapsto v^\flat=\left\langle v,\cdot\right\rangle, \quad \sharp=\flat^{-1}:H'\to H.\] We use \(\left\langle \cdot,\cdot\right\rangle\) for both the inner product on \(H\) and the induced duality pairing between \(H'\) and \(H\). For the free energy \(F:H\to\mathbb{R}\), its differential \(\mathrm dF(\Phi)\in H'\) is a covector, and throughout the paper we write its gradient as \[G(\Phi):=\nabla F(\Phi)=\bigl(\mathrm dF(\Phi)\bigr)^\sharp\in H.\] Let \(\Lambda^2H'\) denote the space of skew-symmetric bilinear forms \(\omega:H\times H\to\mathbb{R}\). For \(\alpha,\beta\in H'\), define their wedge by \[\label{eq:wedge-def} (\alpha\wedge\beta)(u,v)=\alpha(u)\beta(v)-\alpha(v)\beta(u),\tag{6}\] Define the contraction by \(\iota_v\omega=\omega(v,\cdot)\in H'\). In particular, \(\iota_u(\alpha\wedge\beta)=\alpha(u)\beta-\beta(u)\alpha\).
Each \(\omega\in\Lambda^2H'\) defines a skew-symmetric operator \(S_\omega:H\to H\) by \[\label{eq:omega-operator} S_\omega u=(\iota_u\omega)^\sharp, \quad \left\langle S_\omega u,v\right\rangle=\omega(u,v)=-\left\langle u,S_\omega v\right\rangle.\tag{7}\] For a linear map \(\mathcal{S}:H\to H'\), skew-symmetric means \(\left\langle \mathcal{S}u,v\right\rangle=-\left\langle \mathcal{S}v,u\right\rangle\) for all \(u,v\in H\). For \(a,b\in H\), we abbreviate \(a\wedge b:=S_{a^\flat\wedge b^\flat}\), so that \[\label{eq:vector-wedge} (a\wedge b)\,u=\left\langle a,u\right\rangle\,b-\left\langle b,u\right\rangle\,a .\tag{8}\] For the geometric statements, \(\mathcal{M}\) is an open subset of a finite-dimensional real inner-product space \(H\). Affine constraints are handled by identifying the affine space with its tangent space, endowed with the induced inner product. Hence \(T_\Phi\mathcal{M}\simeq H\) and \(T_\Phi^*\mathcal{M}\simeq H'\). The Lie bracket of vector fields is \[[X,Z]=DZ[X]-DX[Z].\] For \(\mathcal{Q}\in C^\infty(\mathcal{M})\), write \(X(\mathcal{Q}):=\mathrm d\mathcal{Q}(X)\); then \([X,Z](\mathcal{Q})=X(Z(\mathcal{Q}))-Z(X(\mathcal{Q}))\). A bivector field \(L\) assigns to each point a skew-symmetric map \(T^*\mathcal{M}\to T\mathcal{M}\). In particular, a decomposable bivector \(Y\wedge J\) acts on covectors by \[\label{eq:bivector-action} (Y\wedge J)\,\alpha=\alpha(Y)\,J-\alpha(J)\,Y .\tag{9}\] It induces the bracket \(\{\mathcal{A},\mathcal{B}\}_L=\left\langle \mathrm d\mathcal{A},L\,\mathrm d\mathcal{B}\right\rangle\) and defines a Poisson structure if for all \(\mathcal{A},\mathcal{B},\mathcal{C}\), \[\{\mathcal{A},\{\mathcal{B},\mathcal{C}\}_L\}_L+\{\mathcal{B},\{\mathcal{C},\mathcal{A}\}_L\}_L+\{\mathcal{C},\{\mathcal{A},\mathcal{B}\}_L\}_L=0.\] Equivalently, in local coordinates, for all \(i,j,k\), \[\label{eq:jacobi-coordinates} \sum\limits_{\ell} L^{i\ell}\partial_\ell L^{jk} +L^{j\ell}\partial_\ell L^{ki} +L^{k\ell}\partial_\ell L^{ij} = 0.\tag{10}\]
We recall the SGE reformulation of 1 .
Theorem 1 ([1]). Assume \(G(\Phi)\ne0\). Then 1 can be written as \[\partial_t\Phi=\bigl(M(\Phi)+S_*(\Phi)\bigr)G(\Phi), \quad S_*(\Phi)=\tfrac{G(\Phi)\wedge J(\Phi)}{\left\|G(\Phi)\right\|^2},\] where \(S_*(\Phi)\) is skew-symmetric and satisfies \[S_*(\Phi)G(\Phi)=J(\Phi).\]
The skew form separates two tasks in time discretization. The dissipative part is controlled by the discretization of the thermodynamic force \(G=\nabla F\), where one may use discrete-gradient methods, stabilization, convex splitting, EQ/SAV-type approaches, averaged-vector-field methods, or supplementary-variable formulations [19]–[34]. The skew part \(S_*G\), which represents the ZEC term, may be treated explicitly by freezing the profiles defining \(S_*\). This explicit treatment preserves the discrete energy cancellation and often decouples multiphysics variables, for example the velocity and phase-field unknowns in CHNS-type systems. Furthermore, \(S_*\) need not be assembled explicitly. By 8 , its action is given by \[\label{eq:matrix-free-S} S_*v=\tfrac{\left\langle G(\Phi),v\right\rangle\,J(\Phi)-\left\langle J(\Phi),v\right\rangle\,G(\Phi)}{\left\|G(\Phi)\right\|^2},\tag{11}\] which requires only two inner products and two vector updates. All gauges constructed below retain this matrix-free structure.
The SGE construction is a special case of the following definition.
Definition 1. A two-form \(\omega\in\Lambda^2H'\) is an admissible ZEC gauge* for 1 at the state \(\Phi\) if \[\label{eq:compatibility} (\iota_{G(\Phi)}\,\omega)(v)=\left\langle J(\Phi),v\right\rangle, \quad \forall v\in H,\tag{12}\] or equivalently, \[S_\omega G(\Phi)=J(\Phi).\] The corresponding GSGE is \[\label{eq:gsge-system} \partial_t\Phi=\bigl(M(\Phi)+S_\omega(\Phi)\bigr)G(\Phi).\tag{13}\] *
For any admissible gauge, \[\tfrac{\mathrm d}{\mathrm dt}F(\Phi) =\left\langle G(\Phi),M(\Phi)G(\Phi)\right\rangle +\omega(G(\Phi),G(\Phi)) =\left\langle G(\Phi),M(\Phi)G(\Phi)\right\rangle \le0,\] since every two-form is skew-symmetric. At states with \(G(\Phi)\ne0\), admissible gauges exist by 1. The remaining task is to characterize the entire admissible family.
Proposition 2. Suppose \(G\ne0\) and \((\iota_G\omega_0)^\sharp=J\). Then every admissible gauge is of the form \[\omega=\omega_0+\zeta, \quad \iota_G\zeta=0,\]
Proof. The compatibility condition is affine in \(\omega\), so two solutions differ by a two-form annihilated by \(\iota_G\). ◻
Thus compatibility with the energy law leaves an affine gauge freedom.
Selecting a representative from this family requires an additional criterion. We consider three: least-squares minimality, preservation of prescribed invariants, and compatibility with Poisson geometry.
Let \(A:H\to H'\) be symmetric positive definite, with norms \[\left\|v\right\|_A^2=\left\langle Av,v\right\rangle\quad\text{on }H, \quad \left\|r\right\|_{A^{-1}}^2=\left\langle r,A^{-1}r\right\rangle\quad\text{on }H'.\] For \(\omega,\eta\in\Lambda^2H'\), define the \(A\)-weighted Frobenius inner product and norm by \[\label{eq:two-form-Frobenius} \left\langle \omega,\eta\right\rangle_{\mathrm F,A} :=\sum_{i,j=1}^{\dim H}\omega(e_i,e_j)\eta(e_i,e_j), \quad \left\|\omega\right\|_{\mathrm F,A}^2:=\left\langle \omega,\omega\right\rangle_{\mathrm F,A},\tag{14}\] where \(\{e_i\}\) is any \(A\)-orthonormal basis. This is the Frobenius structure induced by the \(A\)-inner product on covariant two-tensors, restricted to \(\Lambda^2H'\). The inner-product properties, and hence the norm properties, are immediate once basis independence is established. We verify the latter briefly. Let \(\{f_k\}\) be another \(A\)-orthonormal basis and write \(f_k=\sum_i c_{ik}e_i\). Since \(\sum_k c_{ik}c_{pk}=\delta_{ip}\), bilinearity gives \[\begin{align} \sum_{k,l}\omega(f_k,f_l)\eta(f_k,f_l) &=\sum_{i,j,p,q}\omega(e_i,e_j)\eta(e_p,e_q)\, \delta_{ip}\delta_{jq} \\ &=\sum_{i,j}\omega(e_i,e_j)\eta(e_i,e_j). \end{align}\] Thus 14 is well-defined.
Theorem 3 (Unified least-squares gauge). Let \(X\ne0\) and let \(r\in H'\) satisfy \(r(X)=0\). Among all two-forms \(\omega\in\Lambda^2H'\) with \(\iota_X\omega=r\), the gauge \[\label{eq:A-minimal-form} \omega_{A,r}=\tfrac{AX\wedge r}{\left\langle AX,X\right\rangle}\qquad{(1)}\] uniquely minimizes \(\left\|\cdot\right\|_{\mathrm F,A}\), and \[\left\|\omega_{A,r}\right\|_{\mathrm F,A}^2 =\tfrac{2\left\|r\right\|_{A^{-1}}^2}{\left\|X\right\|_A^2}.\] For \(X=G(\Phi)\), \(r=J(\Phi)^\flat\), and \(A=\flat\), ?? is the two-form corresponding to the SGE gauge in 1. Thus the embedding of [1] is the least-squares gauge in the native metric.
Proof. Let \(e_1=X/\left\|X\right\|_A\) and extend it to an \(A\)-orthonormal basis. Since \(r(X)=0\), contraction of ?? with \(X\) gives \(\iota_X\omega_{A,r}=r\). By 2, every other admissible gauge has the form \(\omega=\omega_{A,r}+\eta\) with \(\iota_X\eta=0\), so \(\eta(e_1,e_j)=0\) for every \(j\). On the other hand, \(\omega_{A,r}(e_i,e_j)=0\) whenever \(i,j\ge2\). Hence \(\left\langle \omega_{A,r},\eta\right\rangle_{\mathrm F,A}=0\), and \[\left\|\omega\right\|_{\mathrm F,A}^2 =\left\|\omega_{A,r}\right\|_{\mathrm F,A}^2 +\left\|\eta\right\|_{\mathrm F,A}^2.\] This proves uniqueness and minimality. The only nonzero components of \(\omega_{A,r}\) have one index equal to \(1\), and \(r(e_1)=0\); therefore 14 gives the stated norm. ◻
Theorem 4 (Regularized gauge). Let \(X\in H\) and let \(r\in H'\) satisfy \(r(X)=0\). For \(\sigma>0\), the problem \[\min_{\omega\in\Lambda^2H'} \left\{ \tfrac12\left\|\iota_X\omega-r\right\|_{A^{-1}}^2 +\tfrac{\sigma}{4}\left\|\omega\right\|_{\mathrm F,A}^2 \right\}\] has the unique solution \[\label{eq:regularized-solution} \omega^{\sigma} =\tfrac{AX\wedge r}{\left\langle AX,X\right\rangle+\sigma}.\qquad{(2)}\] If \(X\ne0\), then \(\omega^\sigma\) converges to the minimum-norm gauge ?? as \(\sigma\downarrow0\).
Proof. If \(X=0\), the unique minimizer is \(\omega^\sigma=0\). Suppose \(X\ne0\). Set \(q=\iota_X\omega\). Then \(q(X)=0\), and 3 reduces the problem to \[\tfrac12\left\|q-r\right\|_{A^{-1}}^2 +\tfrac{\sigma}{2\left\|X\right\|_A^2}\left\|q\right\|_{A^{-1}}^2\] over \(q(X)=0\). Its unique minimizer is \(q=\tfrac{\left\|X\right\|_A^2}{\left\|X\right\|_A^2+\sigma}r\). Substitution into ?? gives ?? . ◻
Remark 5. Contraction of ?? gives \[\iota_X\omega^\sigma =\tfrac{\left\langle AX,X\right\rangle}{\left\langle AX,X\right\rangle+\sigma}\,r, \quad \left\|\iota_X\omega^\sigma-r\right\|_{A^{-1}} \le \tfrac{\sigma}{c}\left\|r\right\|_{A^{-1}}\] whenever \(\left\langle AX,X\right\rangle\ge c>0\). Thus, for \(X=G(\Phi)\) and \(r=J(\Phi)^\flat\), \[S_{\omega^\sigma}G(\Phi) =J(\Phi)+O(\sigma), \quad \left\langle G(\Phi),S_{\omega^\sigma}G(\Phi)\right\rangle=0.\] Choosing \(\sigma=O(\rho^p)\) for a method of order \(p\) therefore preserves its formal order while retaining exact energy cancellation. The regularized formula is defined at \(X=0\), where it reproduces \(r\) only if \(r=0\).
The preceding remark regularizes a compatible profile satisfying \(r(X)=0\). If a computed profile does not satisfy this condition exactly, the same wedge construction gives its nearest compatible correction.
Proposition 6 (Optimal ZEC correction). For \(X\ne0\) and \(r\in H'\), define \[\Pi_X^A r :=\iota_X\!\left(\tfrac{AX\wedge r}{\left\langle AX,X\right\rangle}\right) =r-\tfrac{r(X)}{\left\langle AX,X\right\rangle}AX.\] Then \(\Pi_X^A r\) uniquely minimizes \(\left\|\xi-r\right\|_{A^{-1}}^2\) over all \(\xi\in H'\) satisfying \(\xi(X)=0\).
Proof. The contraction formula gives \((\Pi_X^A r)(X)=0\). Moreover, \(r-\Pi_X^A r\) is a multiple of \(AX\), which is \(A^{-1}\)-orthogonal to every \(\zeta\in H'\) satisfying \(\zeta(X)=0\), since \(\left\langle AX,\zeta\right\rangle_{A^{-1}}=\zeta(X)\). Thus \(\Pi_X^A r\) is the stated orthogonal projection. ◻
If 1 has additional invariants \(\mathcal{C}_1,\ldots,\mathcal{C}_m\), we choose \(\omega\) so that the embedded reversible term preserves each of them. Assume that their gradients are linearly independent at the current state, and write \[G_{\mathcal{C}_\alpha}(\Phi):=\nabla\mathcal{C}_\alpha(\Phi), \quad \alpha=1,\ldots,m.\] Assume that the original ZEC field preserves these invariants, \[\left\langle J(\Phi)^\flat,G_{\mathcal{C}_\alpha}(\Phi)\right\rangle=0, \quad \alpha=1,\ldots,m.\] To construct the force leg, form the Gram system \[B_{\alpha\beta}=\left\langle A G_{\mathcal{C}_\beta}(\Phi),G_{\mathcal{C}_\alpha}(\Phi)\right\rangle, \quad b_\alpha=\left\langle A G(\Phi),G_{\mathcal{C}_\alpha}(\Phi)\right\rangle,\] solve \(B\lambda=b\), and set \[\label{eq:projected-force-leg} \widetilde{G}(\Phi)=G(\Phi)-\sum_{\alpha=1}^m\lambda_\alpha G_{\mathcal{C}_\alpha}(\Phi).\tag{15}\]
Proposition 7. Under the assumptions above, suppose also that \(\left\langle J(\Phi)^\flat,G(\Phi)\right\rangle=0\). If \(\widetilde{G}(\Phi)\ne0\), then the two-form \[\label{eq:projected-invariant-gauge} \omega_{A,\mathcal{C}} =\tfrac{A\widetilde{G}(\Phi)\wedge J(\Phi)^\flat} {\left\langle A\widetilde{G}(\Phi),G(\Phi)\right\rangle}\qquad{(3)}\] satisfies \[\iota_{G(\Phi)}\omega_{A,\mathcal{C}}=J(\Phi)^\flat, \quad \iota_{G_{\mathcal{C}_\alpha}(\Phi)}\omega_{A,\mathcal{C}}=0, \quad \alpha=1,\ldots,m.\]
Proof. The definition gives \(\left\langle A\widetilde{G},G_{\mathcal{C}_\alpha}\right\rangle=0\) and \[\left\langle A\widetilde{G},G\right\rangle=\left\langle A\widetilde{G},\widetilde{G}\right\rangle>0.\] The two contraction identities now follow directly from ?? , the ZEC condition, and \(\left\langle J(\Phi)^\flat,G_{\mathcal{C}_\alpha}(\Phi)\right\rangle=0\). ◻
Remark 8. Suppose in addition that \[\left\langle M(\Phi)G(\Phi),G_{\mathcal{C}_\alpha}(\Phi)\right\rangle=0, \quad \alpha=1,\ldots,m.\] Then each \(\mathcal{C}_\alpha\) is an invariant of the full system, since \[\tfrac{\mathrm d}{\mathrm dt}\mathcal{C}_\alpha(\Phi) =\left\langle M(\Phi)G(\Phi),G_{\mathcal{C}_\alpha}(\Phi)\right\rangle +\left\langle J(\Phi)^\flat,G_{\mathcal{C}_\alpha}(\Phi)\right\rangle=0.\] The same balance holds discretely whenever the temporal discretization retains the corresponding discrete chain rule and orthogonality conditions.
The compatibility condition 12 alone does not imply the Jacobi identity. We now determine when a rank-two gauge defines a Poisson structure.
Consider the rank-two gauge \(L_2=Y\wedge J\). By 9 , \(L_2\,\mathrm dF=J\) if \(\mathrm dF(Y)=1\) and \(\mathrm dF(J)=0\). The Jacobi identity adds a differential condition on \((Y,J)\).
Proposition 9 (Rank-two Jacobi criterion). Let \(Y\) and \(J\) be smooth vector fields. Then \(L_2=Y\wedge J\) defines a Poisson structure if and only if \[\label{eq:frobenius-condition} Y\wedge J\wedge[Y,J]=0,\qquad{(4)}\] Equivalently, \([Y,J]\in\operatorname{span}\{Y,J\}\) wherever \(Y\) and \(J\) are linearly independent; at points where \(Y\wedge J=0\), ?? holds automatically. In particular, \([Y,J]=0\) suffices [35].
Proposition 10 (Local existence). Let \(J(z_0)\ne0\), and suppose that \(\mathrm dF(J)=0\) and \(\mathrm d\mathcal{C}_a(J)=0\), \(a=1,\ldots,m\), near \(z_0\). If \(\mathrm dF,\mathrm d\mathcal{C}_1,\ldots,\mathrm d\mathcal{C}_m\) are linearly independent at \(z_0\), then locally there is a vector field \(Y\) with \[[Y,J]=0, \quad \mathrm dF(Y)=1, \quad \mathrm d\mathcal{C}_a(Y)=0,\quad a=1,\ldots,m,\] and \(L_2=Y\wedge J\) defines a Poisson structure with \[\label{eq:L-casimir} L_2\,\mathrm dF=J, \quad L_2\,\mathrm d\mathcal{C}_a=0,\quad a=1,\ldots,m .\qquad{(5)}\]
Proof. Let \(\Psi_t\) be the local flow of \(J\) and \(\Sigma\) a section transverse to \(J\) at \(z_0\). Independence gives a field \(Y_\Sigma\) on \(\Sigma\) with \(\mathrm dF(Y_\Sigma)=1\) and \(\mathrm d\mathcal{C}_a(Y_\Sigma)=0\). Extend it by \(Y_{\Psi_t(z)}=(\Psi_t)_*Y_\Sigma(z)\), so \([Y,J]=0\). For any \(\mathcal{Q}\) satisfying \(J(\mathcal{Q})=0\), \[J(Y(\mathcal{Q}))=Y(J(\mathcal{Q}))+[J,Y](\mathcal{Q})=0,\] so the identities on \(\Sigma\) propagate along the flow. The conclusion follows from 9 and 9 . ◻
The \(\mathcal{C}_a\) are Casimirs of the constructed bracket. This local construction realizes the selected field \(J=L_2\,\mathrm dF\) and need not coincide with the physical Poisson bracket.
The GENERIC formalism [16]–[18] writes nonequilibrium dynamics as \[\label{eq:generic-form} \dot{z}=L(z)\,\mathrm dE(z)+\mathsf M(z)\,\mathrm dS_{\rm ent}(z), \quad L\,\mathrm dS_{\rm ent}=0, \quad \mathsf M\,\mathrm dE=0, \quad \mathsf M\ge0,\tag{16}\] where \(L\) is Poisson. The two degeneracy conditions state that reversible motion preserves entropy and irreversible motion preserves energy. For the isothermal free-energy convention 1 , ?? gives the corresponding reversible degeneracy. This rank-two bracket reproduces the selected field \(J\) but need not encode the physical Lie–Poisson structure. The Navier–Stokes example below carries the construction to the fully discrete level.
Let \(\Omega\) be a periodic box and \(P_\sigma\) the Helmholtz–Leray projection. The incompressible Navier–Stokes equations are [36] \[\label{eq:NS} \partial_t\boldsymbol{u}+P_\sigma(\boldsymbol{u}\cdot\nabla\boldsymbol{u})=\nu P_\sigma\Delta\boldsymbol{u}.\tag{17}\] On the divergence-free space, set \(F(\boldsymbol{u})=\tfrac12\left\|\boldsymbol{u}\right\|^2\), \(M\nabla F=\nu P_\sigma\Delta\boldsymbol{u}\), and \(J_E(\boldsymbol{u})=-P_\sigma(\boldsymbol{u}\cdot\nabla\boldsymbol{u})\). Then \((J_E(\boldsymbol{u}),\boldsymbol{u})=0\), and the SGE gauge is \(S_*=\tfrac{\boldsymbol{u}\wedge J_E}{\left\|\boldsymbol{u}\right\|^2}\). The Euler part also carries the noncanonical Lie–Poisson bracket [37], [38] \[\label{eq:euler-lie-poisson} \{F,G\}_{\rm Eul}(\boldsymbol{u}) =-\int_\Omega\boldsymbol{u}\cdot \Bigl[\tfrac{\delta F}{\delta\boldsymbol{u}},\tfrac{\delta G}{\delta\boldsymbol{u}}\Bigr]\mathrm dx ,\tag{18}\] where \([\cdot,\cdot]\) is the Lie bracket of divergence-free vector fields. The following continuum calculation motivates the MAC construction below.
Proposition 11 (Formal continuum rank-two gauge). On \(\{\boldsymbol{u}\ne0\}\), the bivector \(L_2^{NS}=Y\wedge J_E\), with \(Y(\boldsymbol{u})=\tfrac{\boldsymbol{u}}{\left\|\boldsymbol{u}\right\|^2}\), formally satisfies the Jacobi identity and \(L_2^{NS}\,\mathrm dF=J_E\).
Proof. The identity \(L_2^{NS}\mathrm dF=J_E\) follows from \(\mathrm dF(Y)=1\), \(\mathrm dF(J_E)=0\), and 9 . Since \(J_E\) is quadratic and energy neutral, \[DJ_E[Y]=\tfrac{2J_E}{\left\|\boldsymbol{u}\right\|^2}, \quad DY[J_E]=\tfrac{J_E}{\left\|\boldsymbol{u}\right\|^2} -\tfrac{2(\boldsymbol{u},J_E)\boldsymbol{u}}{\left\|\boldsymbol{u}\right\|^4} =\tfrac{J_E}{\left\|\boldsymbol{u}\right\|^2}.\] Thus \([Y,J_E]=\tfrac{J_E}{\left\|\boldsymbol{u}\right\|^2}\), and 9 applies formally. ◻
This rank-two bracket generates \(J_E\) from \(F\) but is not the physical bracket 18 . The MAC construction retains the bilinearity and energy neutrality used above.
On a uniform two-dimensional MAC grid [39], [40], store \(p\) at cell centers and \((u,v)\) at the corresponding face centers. With the mesh-weighted inner products, define \[\begin{align} (D_h\boldsymbol{u})_{i,j} & =\tfrac{u_{i+1/2,j}-u_{i-1/2,j}}{h} +\tfrac{v_{i,j+1/2}-v_{i,j-1/2}}{h}, \\ (G_hp)^x_{i+1/2,j} & =\tfrac{p_{i+1,j}-p_{i,j}}{h}, \quad (G_hp)^y_{i,j+1/2} =\tfrac{p_{i,j+1}-p_{i,j}}{h}. \end{align}\] Periodic summation by parts gives \(G_h=-D_h^{\top}\), equivalently \((D_h\boldsymbol{v},q)_h=-(\boldsymbol{v},G_hq)_h\) [41], [42]. Set \(V_h=\ker D_h\), let \(P_h\) be the orthogonal projection onto \(V_h\), and define the componentwise Laplacian by \[\begin{align} (\Delta_h u)_{i+1/2,j} & =\tfrac{u_{i+3/2,j}-2u_{i+1/2,j}+u_{i-1/2,j}}{h^2} +\tfrac{u_{i+1/2,j+1}-2u_{i+1/2,j}+u_{i+1/2,j-1}}{h^2}, \\ (\Delta_h v)_{i,j+1/2} & =\tfrac{v_{i+1,j+1/2}-2v_{i,j+1/2}+v_{i-1,j+1/2}}{h^2} +\tfrac{v_{i,j+3/2}-2v_{i,j+1/2}+v_{i,j-1/2}}{h^2}. \end{align}\]
For face fields \(\boldsymbol{w}\) and \(\boldsymbol{z}\), let \(N_h(\boldsymbol{w})\boldsymbol{z}\) and \(K_h(\boldsymbol{w})\boldsymbol{z}\) be the centered advective and conservative MAC approximations, with arithmetic averages placing products on the required faces [40]. Define \[\label{eq:Ch-def} C_h(\boldsymbol{w})\boldsymbol{z} =\tfrac12\bigl(N_h(\boldsymbol{w})\boldsymbol{z}+K_h(\boldsymbol{w})\boldsymbol{z}\bigr) \approx\tfrac12\bigl((\boldsymbol{w}\cdot\nabla)\boldsymbol{z}+\operatorname{div}(\boldsymbol{z}\otimes\boldsymbol{w})\bigr).\tag{19}\] Periodic summation by parts gives \(K_h(\boldsymbol{w})=-N_h(\boldsymbol{w})^{\top}\), hence \(C_h(\boldsymbol{w})=\tfrac12\bigl(N_h(\boldsymbol{w})-N_h(\boldsymbol{w})^{\top}\bigr)\) and \(C_h(\boldsymbol{w})^{\top}=-C_h(\boldsymbol{w})\). If \(D_h\boldsymbol{w}=0\), 19 is second-order consistent with \((\boldsymbol{w}\cdot\nabla)\boldsymbol{z}\). The semi-discrete scheme is \[\label{eq:NS-semidiscrete} \dot{\boldsymbol{u}}=\nu P_h\Delta_h\boldsymbol{u}+J_h(\boldsymbol{u}) \quad\text{on }V_h,\tag{20}\] where \(J_h(\boldsymbol{u})=-P_hC_h(\boldsymbol{u})\boldsymbol{u}\) and \(F_h(\boldsymbol{u})=\tfrac12\left\|\boldsymbol{u}\right\|_h^2\). The operator \(\Delta_h\) is symmetric negative semidefinite, and \(\left\|\nabla_h\boldsymbol{v}\right\|_h^2:=-(\Delta_h\boldsymbol{v},\boldsymbol{v})_h\). Skewness and linearity of \(C_h\) give \[\label{eq:Jh-properties} (J_h(\boldsymbol{u}),\boldsymbol{u})_h=-(C_h(\boldsymbol{u})\boldsymbol{u},\boldsymbol{u})_h=0, \quad J_h(\lambda\boldsymbol{u})=\lambda^2J_h(\boldsymbol{u}).\tag{21}\]
Theorem 12 (Semi-discrete rank-two structure). On \(V_h\setminus\{0\}\), the bivector \(L_2^h=Y_h\wedge J_h\), with \(Y_h(\boldsymbol{u})=\tfrac{\boldsymbol{u}}{\left\|\boldsymbol{u}\right\|_h^2}\), is Poisson and satisfies \(L_2^h\,\mathrm dF_h=J_h\).
Proof. By 21 , \(\mathrm dF_h(Y_h)=1\) and \(\mathrm dF_h(J_h)=0\), so 9 gives \(L_2^h\,\mathrm dF_h=J_h\). Quadratic homogeneity of \(J_h\) and \((\boldsymbol{u},J_h)_h=0\) give \[DJ_h[Y_h]=\tfrac{2J_h}{\left\|\boldsymbol{u}\right\|_h^2}, \quad DY_h[J_h] =\tfrac{J_h}{\left\|\boldsymbol{u}\right\|_h^2}-\tfrac{2(\boldsymbol{u},J_h)_h\,\boldsymbol{u}}{\left\|\boldsymbol{u}\right\|_h^4} =\tfrac{J_h}{\left\|\boldsymbol{u}\right\|_h^2}.\] Hence \([Y_h,J_h]=\tfrac{J_h}{\left\|\boldsymbol{u}\right\|_h^2}\), so 9 applies. ◻
With \(\boldsymbol{u}^{n+1/2}=\tfrac12(\boldsymbol{u}^n+\boldsymbol{u}^{n+1})\), the implicit midpoint discretization of 20 is \[\label{eq:NS-midpoint} \tfrac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^n}{\tau} =\nu P_h\Delta_h\boldsymbol{u}^{n+1/2}+J_h(\boldsymbol{u}^{n+1/2}) .\tag{22}\]
Theorem 13. At each midpoint, \(J_h=L_2^h\,\mathrm dF_h\); hence 22 is a GENERIC discretization in the rank-two sense. It satisfies the fully discrete energy law \[\label{eq:NS-discrete-energy} F_h(\boldsymbol{u}^{n+1})-F_h(\boldsymbol{u}^n) =-\nu\tau\left\|\nabla_h\boldsymbol{u}^{n+1/2}\right\|_h^2 .\qquad{(6)}\]
Remark 14. The midpoint rule is not generally a Poisson map for the state-dependent bracket \(L_2^h\) [43], and the MAC scheme does not preserve the full Euler bracket 18 . The GENERIC statement in 13 concerns the rank-two degeneracy and energy law, not preservation of either bracket.
On a periodic domain \(\Omega\), consider the CHNS system [40], [44]–[52] \[\label{eq:chns-system} \left\{ \begin{align} & \partial_t\boldsymbol{u}+P_\sigma (\boldsymbol{u}\cdot\nabla)\boldsymbol{u}+P_\sigma \phi\nabla\bar\mu=\nu P_\sigma\Delta\boldsymbol{u}, \\ & \partial_t\phi+\operatorname{div}(\phi\boldsymbol{u})=m\Delta\bar\mu, \\ & \mu=-\gamma\varepsilon\Delta\phi+\tfrac{\gamma}{\varepsilon}f'(\phi), \\ & \bar\mu=\mu-\tfrac{1}{\left|\Omega\right|}\int_\Omega\mu\,\mathrm dx, \\ & f(\phi) = \tfrac{1}{4}(\phi^2 - 1)^2. \end{align} \right.\tag{23}\] with mobility \(m>0\) and free energy \[F(\boldsymbol{u},\phi)=\tfrac12\left\|\boldsymbol{u}\right\|^2 +\tfrac{\gamma\varepsilon}{2}\left\|\nabla\phi\right\|^2 +\tfrac{\gamma}{\varepsilon}(f(\phi),1).\] On the divergence-free space, the reversible field \[J(\Phi)=-\bigl(P_\sigma(\boldsymbol{u}\cdot\nabla)\boldsymbol{u}+P_\sigma\phi\nabla\bar\mu,\;\operatorname{div}(\phi\boldsymbol{u})\bigr)\] is energy neutral for \(G(\Phi)=(\boldsymbol{u},\bar\mu)\) by periodic integration by parts.
Set \(D_2a^{n+1}:=\tfrac{3a^{n+1}-4a^n+a^{n-1}}{2\tau}\) and \(\hat{a}^{n+1}:=2a^n-a^{n-1}\) for \(a=\boldsymbol{u},\phi\). In continuous spatial notation, with a compatible periodic summation-by-parts discretization understood, the GSGE–BDF2 scheme is \[\label{eq:chns-scheme} \left\{ \begin{align} & D_2\boldsymbol{u}^{n+1}+\nabla p^{n+1}=\nu\Delta\boldsymbol{u}^{n+1}-\mathcal{R}_{\boldsymbol{u}}^{n+1}, \\ & \operatorname{div}\boldsymbol{u}^{n+1}=0, \\ & D_2\phi^{n+1}=m\Delta\bar\mu^{n+1}-\mathcal{R}_{\phi}^{n+1}, \\ & \mu^{n+1}=-\gamma\varepsilon\Delta\phi^{n+1} +\tfrac{\gamma}{\varepsilon} [\chi(\tfrac{3 \phi^{n+1} - \phi^n}{2}, \tfrac{3 \phi^n - \phi^{n-1}}{2}) - \hat{\phi}^{n+1}], \\ & \chi (a, b) := \tfrac{1}{4}(a^2 + b^2)(a + b), \end{align} \right.\tag{24}\] The ZEC residual uses the regularized gauge ?? , with \(\sigma=\ell_3\tau^2\) and the differential weight specified below: \[\bigl(\mathcal{R}_{\boldsymbol{u}}^{n+1},\mathcal{R}_{\phi}^{n+1}\bigr) =\lambda_A^{n+1}\,\hat{J}^{n+1} -\lambda_J^{n+1}\,A_{\ell_1,\ell_2} \hat{G}^{n+1},\] where \[\begin{align} \hat{J}^{n+1} & =\bigl((\hat{\boldsymbol{u}}^{n+1}\!\cdot\!\nabla)\hat{\boldsymbol{u}}^{n+1}+\hat{\phi}^{n+1}\nabla\hat{\mu}^{n+1},\;\operatorname{div}(\hat{\phi}^{n+1}\hat{\boldsymbol{u}}^{n+1})\bigr), \\ \hat{G}^{n+1} & =(\hat{\boldsymbol{u}}^{n+1},\bar{\hat{\mu}}^{n+1}), \\ A_{\ell_1,\ell_2} \hat{G}^{n+1} & =\bigl(\ell_1^2\hat{\boldsymbol{u}}^{n+1}-\ell_2^2\Delta\hat{\boldsymbol{u}}^{n+1},\;\ell_1^2\bar{\hat{\mu}}^{n+1}-\ell_2^2\Delta\hat{\mu}^{n+1}\bigr), \\ \hat{\mu}^{n+1} & = -\gamma \varepsilon \Delta \hat{\phi}^{n+1} + \tfrac{\gamma}{\varepsilon} \left((\hat{\phi}^{n+1})^3 - \hat{\phi}^{n+1}\right). \end{align}\] Let \(e^{n+1}=(\boldsymbol{u}^{n+1},\bar\mu^{n+1})\) and \[\label{eq:chns-denominator} \begin{align} d^{\,n+1} & =\left\langle A_{\ell_1,\ell_2}\hat{G}^{n+1},\hat{G}^{n+1}\right\rangle+\ell_3\tau^2 \\ & =\ell_1^2\left\|\hat{\boldsymbol{u}}^{n+1}\right\|^2+\ell_1^2\left\|\bar{\hat{\mu}}^{n+1}\right\|^2 +\ell_2^2\left\|\nabla\hat{\boldsymbol{u}}^{n+1}\right\|^2 \\ & \quad+\ell_2^2\left\|\nabla\hat{\mu}^{n+1}\right\|^2+\ell_3\tau^2. \end{align}\tag{25}\] Then \[\label{eq:chns-lambdas} \lambda_A^{n+1}=\tfrac{\left\langle A_{\ell_1,\ell_2}\hat{G}^{n+1},e^{n+1}\right\rangle}{d^{\,n+1}}, \quad \lambda_J^{n+1}=\tfrac{\left\langle \hat{J}^{n+1},e^{n+1}\right\rangle}{d^{\,n+1}} .\tag{26}\]
Remark 15. Let \(\ell_1,\ell_2,\ell_3\ge0\) with \(\ell_1^2+\ell_2^2>0\). If \(\ell_3>0\), then \(d^{\,n+1}>0\) even when the extrapolated force vanishes. Where the unregularized denominator is bounded away from zero, \(\ell_3\tau^2=O(\tau^2)\) retains second-order consistency by 5. If \(\ell_2=\ell_3=0\), the factors \(\ell_1^2\) cancel and the closure reduces to the SGE–SBDF2 method of [1]. If \(\ell_1=\ell_3=0\), the phase component of \(A_{\ell_1,\ell_2}\hat{G}^{n+1}\) is \(-\ell_2^2\Delta\hat{\mu}^{n+1}\) and has zero mean without projecting \(\hat{\mu}^{n+1}\); the displayed skew formula therefore gives a natural mass-preserving gauge whenever \(d^{\,n+1}>0\), although the differential weight is semidefinite in this limit. The rank-two GSGE closure weakly couples the Navier–Stokes and phase-field subproblems only through the two scalar coefficients \(\lambda_A^{n+1}\) and \(\lambda_J^{n+1}\). Hence 24 admits the decoupled implementation of [1]. Since the efficiency of the SGE framework and the BDF2 discretization has already been demonstrated in [1], [53], we do not present additional numerical results here.
Proposition 16. The scheme 24 –26 satisfies \[(\mathcal{R}_\phi^{n+1},1)=0, \quad (\mathcal{R}_{\boldsymbol{u}}^{n+1},\boldsymbol{u}^{n+1})+(\mathcal{R}_{\phi}^{n+1},\bar\mu^{n+1})=0,\] and \((\phi^{n+1},1)=(\phi^0,1)\) for all \(n\) if \((\phi^1,1)=(\phi^0,1)\).
Proof. Periodicity gives \((\operatorname{div}(\hat{\phi}^{n+1}\hat{\boldsymbol{u}}^{n+1}),1)=0\) and \((\ell_1^2\bar{\hat{\mu}}^{n+1}-\ell_2^2\Delta\hat{\mu}^{n+1},1)=0\), hence \((\mathcal{R}_\phi^{n+1},1)=0\). Testing the phase equation with \(1\) gives mass conservation. Moreover, \[\begin{align} (\mathcal{R}_{\boldsymbol{u}}^{n+1},\boldsymbol{u}^{n+1}) +(\mathcal{R}_{\phi}^{n+1},\bar\mu^{n+1}) & =\left\langle \lambda_A^{n+1}\hat{J}^{n+1} -\lambda_J^{n+1}A_{\ell_1,\ell_2} \hat{G}^{n+1},e^{n+1}\right\rangle \\ & =\lambda_A^{n+1}d^{\,n+1}\lambda_J^{n+1} -\lambda_J^{n+1}d^{\,n+1}\lambda_A^{n+1}=0 . \end{align}\] ◻
Theorem 17. Under the compatible spatial discretization, the GSGE–BDF2 scheme satisfies \[\begin{align} & \widetilde{F}^{\,n+1} - \widetilde{F}^{\,n} + \tau m \left\|\nabla\mu^{n+1}\right\|^2 + \nu \tau \left\|\nabla\boldsymbol{u}^{n+1}\right\|^2 \\ & \quad + \tfrac{1}{4} \left\|\boldsymbol{u}^{n+1} - 2 \boldsymbol{u}^n + \boldsymbol{u}^{n-1}\right\|^2 + \tfrac{\gamma \varepsilon}{4} \left\|\nabla(\phi^{n+1} - 2\phi^n + \phi^{n-1})\right\|^2 \\ & \quad + \tfrac{3 \gamma}{4 \varepsilon} \left\|\phi^{n+1} - 2\phi^n + \phi^{n-1}\right\|^2 = 0. \end{align}\] Here \[\begin{align} \widetilde{F}^{\,n+1} & =\tfrac14\left\|\boldsymbol{u}^{n+1}\right\|^2+\tfrac14\left\|2\boldsymbol{u}^{n+1}-\boldsymbol{u}^n\right\|^2 +\tfrac{\gamma\varepsilon}{4}\left\|\nabla\phi^{n+1}\right\|^2 +\tfrac{\gamma\varepsilon}{4}\left\|\nabla(2\phi^{n+1}-\phi^n)\right\|^2 \\ & \quad+\tfrac{\gamma}{\varepsilon}\bigl(f(\tfrac{3\phi^{n+1} - \phi^n}{2}),1\bigr) + \tfrac{3 \gamma}{8 \varepsilon} \left\|\phi^{n+1} - \phi^n\right\|^2. \end{align}\]
Proof. Test the momentum, phase, and potential equations with \(2\tau\boldsymbol{u}^{n+1}\), \(2\tau\bar\mu^{n+1}\), and \(2\tau D_2\phi^{n+1}\), respectively. The pressure term vanishes, and the time differences satisfy \[\begin{align} 2\tau(D_2a^{n+1},a^{n+1}) & =\tfrac12\bigl[\left\|a^{n+1}\right\|^2+ \left\|2a^{n+1}-a^n\right\|^2\bigr] \\ & \quad-\tfrac12\bigl[\left\|a^{n}\right\|^2+ \left\|2a^{n}-a^{n-1}\right\|^2\bigr] +\tfrac12\left\|a^{n+1}-2a^n+a^{n-1}\right\|^2. \end{align}\] Mass conservation gives \((D_2\phi^{n+1},\bar\mu^{n+1})=(D_2\phi^{n+1},\mu^{n+1})\), while 16 cancels the reversible terms: \[2\tau\bigl[(\mathcal{R}_{\boldsymbol{u}}^{n+1},\boldsymbol{u}^{n+1}) +(\mathcal{R}_{\phi}^{n+1},\bar\mu^{n+1})\bigr]=0.\] The remaining terms give the stated identity as in [53]. ◻
We developed a generalized skew-gradient embedding (GSGE) framework for thermodynamically consistent systems with zero-energy contributions. The compatibility condition determines an affine space of admissible skew two-forms. Weighted least squares selects a unique representative and recovers SGE in the native metric; regularization controls vanishing force profiles, and projection enforces prescribed invariants. For rank-two gauges, the Jacobi criterion identifies representatives that define low-rank Poisson structures and satisfy the corresponding GENERIC degeneracy conditions.
For the incompressible Navier–Stokes equations, a compatible MAC discretization yields a finite-dimensional rank-two Poisson–GENERIC formulation, and the implicit midpoint rule satisfies the exact discrete energy law. For the CHNS system, the GSGE–BDF2 scheme uses a three-parameter regularized differential gauge. The parameter \(\ell_3\) prevents a vanishing denominator, while the limits \(\ell_2=\ell_3=0\) and \(\ell_1=\ell_3=0\) recover the SGE gauge and a gradient-weighted mass-preserving gauge, respectively. The latter requires no mean projection of the chemical-potential profile. The rank-two closure also permits a decoupled solution of the Navier–Stokes and phase-field subproblems through two scalar coefficients. The scheme preserves mass and satisfies the unconditional discrete energy law.
Future work includes selecting the gauge parameters for efficiency and accuracy, extending regularized differential gauges to higher-order schemes, and constructing fully discrete low-rank Poisson gauges for more general systems.
No data were generated in this work.
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.