How can the surface Stokes-Cahn-Hilliard system be analyzed and discretized?

A membrane model that drops inertia but keeps geometry

The system couples an inertia-free surface Stokes momentum equation to a Cahn-Hilliard phase equation through advection and capillary forces. The unknowns are a tangential velocity u, a reduced pressure p, a phase field c, and a chemical potential μ, posed on a closed, stationary hypersurface Γ of class C2 embedded in R^(d+1) with d=2,3. The parameters are constant viscosity η>0, line (surface for d=3) tension σ>0, mobility M>0, and interface width ε>0.

The inertia-free form is not an arbitrary simplification. The paper states that it follows from the surface Navier-Stokes-Cahn-Hilliard model by neglecting inertia at small surface Reynolds numbers, an assumption the authors describe as consistent with the material parameters of multicomponent lipid bilayers. The same inertia-free surface flow equations also arise in one-phase membrane models with Helfrich bending forces. The intended physical setting is lateral phase separation in viscous fluid membranes, motivated by the lipid raft hypothesis, coexisting fluid domains in multicomponent vesicles, and engineered liposomes and lipid bilayers.

Killing fields leave the velocity determined only up to a rigid motion

The surface viscous operator has the Killing fields as its kernel. These are the tangential fields with vanishing surface rate of strain, defined as E = {v in H1(Γ)^(d+1) : v·n = 0, D_Γ v = 0 on Γ}. Because they lie in the kernel, the Stokes momentum equation determines the velocity only up to the addition of a Killing field. The paper removes this indeterminacy by seeking the tangential velocity that is L2(Γ)-orthogonal to all Killing fields, and it shows that the capillary force is orthogonal to the Killing fields, which makes the momentum equation compatible.

This choice is not cosmetic for the coupled problem. Adding a Killing field leaves the viscous stress unchanged but generally changes the transport of the phase field, so the selection rule changes the solution. The authors also note that an alternative used elsewhere is to add a zeroth-order friction term γu with constant γ>0 to the momentum equation. The paper's own route is the orthogonality constraint u ⟂ E, and the discrete method must control spurious velocity components associated with these fields.

Global weak solutions exist and satisfy an energy equality

For constant viscosity and mobility and a class of regular double-well potentials, the paper proves global existence and uniqueness of weak solutions for initial phase fields in H1(Γ). The double-well potential f0 is assumed to satisfy f0 ≥ 0 and the growth bound |f0'''(s)| ≤ L(1+|s|) for a.e. s, which the authors note is satisfied by the Ginzburg-Landau potential f_quart(s) = (s^2-1)^2/4. Existence is proved by a Faedo-Galerkin scheme with a priori bounds, compactness via the Aubin-Lions-Simon lemma, and passage to the limit.

The weak solution satisfies an energy equality and conserves mass: (c(t),1) = (c0,1) for all t in [0,T]. The energy inequality stated in the existence theorem bounds E[c(t)] plus the time integral of 2η‖D_Γ u‖^2 + σM‖∇_Γ μ‖^2 by E[c0] for a.e. t. The paper also establishes Lipschitz continuous dependence in energy norms on every finite time interval, including for initial data with different means, and derives higher regularity results for d=2 under additional smoothness assumptions on the surface, potential, and initial data.

A centered convection form removes the need for pointwise divergence-free velocity

The weak formulation uses a centered trilinear form b(u,v,w) = (u·∇_Γ v, w#), where w# = w - ⟨w⟩_Γ is the mean-zero part. The paper shows that b(u,v,1) = 0 for all u,v, and that when u·n = 0 and div_Γ u = 0, b agrees with the usual non-centered convection form and is skew-symmetric in its second and third arguments.

Using the centered form b(w,c,μ) in the momentum equation instead of the term corresponding to the strong form amounts only to a pressure shift p ↦ p + σ(cμ - ⟨μ⟩_Γ c)#. The paper explains that this consistent modification creates a direct cancellation between the capillary force term b(w,c,μ) and the convective term b(u,c,v) in the phase field equation when testing with w = u and v = μ, thereby circumventing the need for div_Γ u = 0. The authors describe this as critical for the structure-preserving numerical method, because the divergence constraint is hard to enforce pointwise in the discrete setting.

The unfitted finite element method conserves phase and energy exactly

The discretization uses the trace finite element method (TraceFEM) with an implicitly represented surface and a background mesh that need not fit the surface. Velocity and pressure use bulk Taylor-Hood elements, while the phase field and chemical potential use continuous finite elements. The method combines a penalty on approximate Killing fields with projections in the transport and capillary terms to control these modes while preserving the energy balance.

The fully discrete scheme uses a midpoint time stepping with a secant approximation of the potential derivative, motivated by the rapid spinodal decomposition stage followed by slow coarsening in Cahn-Hilliard dynamics. Testing with w_h = u_h^(n+1/2), q_h = p_h^(n+1/2), v_h = σμ_h^(n+1/2), and g_h = σδ_τ c_h^(n+1) gives the exact discrete energy equality δ_τ E_h[c_h^(n+1)] + 2η a_h(u_h^(n+1/2), u_h^(n+1/2)) + s_0^ext(p_h^(n+1/2), p_h^(n+1/2)) + σM a_h(μ_h^(n+1/2), μ_h^(n+1/2)) = 0. The paper states that this identity holds for every τ_(n+1) > 0 for which a solution of the discrete system exists, so time discretization introduces no additional dissipation beyond that present in the semidiscrete energy law. Testing the phase equation with v_h = 1 gives total phase conservation: (c_h^(n+1),1)_Γh = (c_h^n,1)_Γh = (c_h^0,1)_Γh.

The paper contrasts this with an alternative regularization and filtering approach that adds an artificial friction term δu with δ = h^α, 1 ≤ α ≤ 2, and identifies approximate Killing fields through a discrete Stokes eigenvalue problem. The authors note that artificial friction introduces additional dissipation and thus modifies the energy balance of the friction-free model, and that removing the Killing component only from the transport velocity generally leaves an extra term in the energy balance. Their method instead applies the same projection in the transport and capillary terms and penalizes the approximate Killing component of the velocity.

Approximate Killing fields come from small eigenvalues or ambient rigid motions

The paper recovers approximate Killing fields from eigenmodes associated with small eigenvalues of a discrete surface viscous operator. For d=2, it also introduces an alternative approach based on restrictions of ambient rigid motions to Γ. That alternative reduces the identification to a six-dimensional eigenvalue problem and recovers the full Killing space when all Killing fields are extrinsic.

The numerical experiments examine the identification of Killing fields and illustrate phase separation on a torus and a liposome with experimentally informed material parameters. The liposome simulation uses a PC-DPPC-cholesterol liposome of composition 25:50:25 mol%, with snapshots ordered from left to right and top to bottom, colors showing the liquid-ordered phase fraction c and arrows indicating the velocity. The paper also notes that eigenvalue-based identification and removal of Killing fields for the surface Stokes problem was previously studied, and that for the coupled problem the new element is the penalty on the approximate Killing component.

What the guarantees cover and where they stop

The well-posedness and continuous dependence results are proved for constant viscosity and mobility and for a class of regular double-well potentials satisfying the stated growth bound. The analysis assumes a closed, stationary C2 hypersurface, and the higher regularity results are derived only for d=2 and only under additional smoothness assumptions on the surface, potential, and initial data. The paper explicitly notes that the weak surface NSCH theory in prior work assumes C3 geometry and the strong theory uses C6 geometry, while this paper assumes C2 surface smoothness for the well-posedness analysis and accounts explicitly for the geometric assumptions needed for higher regularity.

The exact discrete energy identity holds for the stabilized discretization and for every time step for which a solution of the discrete system exists. The paper also notes that midpoint time stepping provides weak damping of stiff dissipative modes, which can reduce accuracy for rough initial data, and that a Rannacher startup replacing the first time interval by two backward Euler half-steps can damp these modes. The startup steps conserve the total phase but do not satisfy the exact energy identity. The numerical experiments are limited to a torus and a liposome, so the reported phase separation behavior is demonstrated in those two geometries rather than across a broad class of surfaces.

The decisive choice is the Killing-field treatment, not the mesh

For readers comparing discretizations of the surface Stokes-Cahn-Hilliard system, the paper's central distinction is how the Killing-field component is handled. The paper's method penalizes the approximate Killing component of the velocity and applies the same projection in the transport and capillary terms, which yields the cancellation between the coupling terms and preserves the energy balance of the stabilized semidiscrete problem. The alternative friction-based approach adds dissipation and modifies the energy balance, and removing the Killing component only from the transport velocity generally leaves an extra term in the energy balance.

The practical priority that follows from the paper's evidence is therefore to check whether a candidate scheme controls the Killing component in both the transport and capillary terms while keeping an exact discrete energy identity and total phase conservation. The paper's own guarantees are tied to constant viscosity and mobility, regular double-well potentials, a closed stationary C2 surface, and, for higher regularity, d=2 with extra smoothness. Within those conditions, the combination of global weak well-posedness, Lipschitz continuous dependence in energy norms, and an unfitted finite element method with exact discrete energy identity is the paper's specific contribution.

Sources used in this answer