Energy-Stable PFEM in Geometric Evolution

• Energy-stable PFEM is a numerical method that uses parameterized finite elements to directly simulate curvature-driven interface evolution while rigorously preserving energy dissipation and conservation laws.
• It employs robust variational formulations and discretization schemes, such as mass-lumping and backward Euler, to ensure unconditional energy stability across diverse geometric flows.
• The method effectively handles anisotropy, solid-state dewetting, and high-order time integration, making it a versatile tool for advanced simulations in materials science and geometric modeling.

An energy-stable parametric finite element method (PFEM) is a class of numerical algorithms developed to simulate geometric evolution problems—such as curvature-driven interface motion, surface diffusion, Willmore flow, and solid-state dewetting—by directly tracking the evolution of interfaces using parameterized finite element representations. The defining property of an energy-stable PFEM is the rigorous preservation at the discrete level of two fundamental features of the underlying continuum model: dissipation of the total interfacial energy (or another relevant energy functional) and the conservation of physical invariants (such as enclosed area or volume), both unconditionally with respect to time step size. Recent advances have extended the methodology to a broad array of geometric flows with isotropic and anisotropic surface energies, for both planar curves and surfaces in higher dimensions, including flows with contact point/migration, solid-state dewetting, and flows involving kinetic or elastic effects.

1. Mathematical Foundations and Model Formulations

Energy-stable parametric FEMs are designed for sharp-interface models governed by geometric evolution equations of the form

$$V_n = \text{GeometricOperator}[\kappa, \mathbf{n}, \gamma(\theta), \ldots],$$

where $V_n$ is the normal velocity of the evolving curve or surface, $\kappa$ the (mean) curvature, $\mathbf{n}$ the outward normal, and possibly an anisotropic surface energy $\gamma(\theta)$. For surface diffusion, the prototypical equation is

$$\partial_t \mathbf{X} = -\partial_s^2\mu\,\mathbf{n}, \qquad \mu = (\gamma(\theta) + \gamma^{\prime\prime}(\theta)) \kappa,$$

and for Willmore flow,

$$\partial_t \mathbf{X}\cdot \mathbf{n} = -\Delta_s \kappa - \frac{1}{2}\kappa^3.$$

Boundary and contact line conditions are often encoded via Robin-type or variationally natural terms to allow for mass (area/volume) conservation and to couple contact line mechanics (e.g., dynamic contact angle) to the bulk dynamics (Zhao et al., 2020, Zhao et al., 2020, Bao et al., 2020).

The energy dissipation law at the continuum level takes the form

$\frac{dW}{dt} = -\mathcal{D}(t) \le 0,$

where $W$ is the total energy (combining interfacial, line, and sometimes bulk energies), and $\mathcal{D}(t)$ denotes a sum of positive-definite dissipation terms (e.g., $L^2$ norm of $\partial_s \mu$, kinetic energy dissipation, contact line friction, etc).

2. Variational and Weak Formulation

The construction of energy-stable PFEMs is anchored on variational formulations of the evolution equations. This involves casting the evolution law and any auxiliary equations (such as for curvature or chemical potential) into a system of weak formulations suitable for finite element approximation:

• For surface diffusion (isotropic/anisotropic):

$$(n \cdot \partial_t X, \phi)_\Gamma + (\partial_s\mu, \partial_s\phi)_\Gamma = 0$$

$$(\mu n, \omega)_\Gamma - (G_k(n)\partial_s X, \partial_s \omega)_\Gamma = 0,$$

where $G_k(n)$ is a (possibly symmetrized) surface energy matrix encoding anisotropy and stabilization (Li et al., 2020, Bao et al., 2021, Bao et al., 2022, Bao et al., 2023, Zhang et al., 2 Apr 2024, Li et al., 28 Jan 2025).

• For Willmore flow:

$$V{\bf n} = \partial_s(\partial_s \kappa {\bf n} - \frac{1}{2}\kappa^2 \partial_s {\bf X})$$

$$\partial_t\kappa = -\partial_s({\bf n}\cdot\partial_s(\partial_t {\bf X})) - (\partial_s{\bf X}\cdot \partial_s(\partial_t {\bf X}))\,\kappa,$$

together with a normal velocity update ${\bf n}\cdot\partial_t {\bf X} = V$ (Bao et al., 24 Jan 2024, Bao et al., 26 Jun 2025, Garcke et al., 30 Jun 2025).

Weak treatment of boundary and contact point conditions—for instance, a Robin-type formulation of the relaxed contact angle for solid-state dewetting—is crucial for ensuring energy dissipation and efficient numerical handling (Zhao et al., 2020, Bao et al., 1 Oct 2024).

3. Discretization Schemes and Energy Stability

The full discretization in space is carried out via parametric finite elements, usually with piecewise linear ($\mathbb{P}_1$) elements on a reference domain (curve parameter or triangulated surface). Integration is performed using mass-lumped inner products for stability and efficiency: $$(u, v)_{\Gamma^h}^h = \frac{1}{2}\sum_j |h_j| [u(\rho_j^-)v(\rho_j^-) + u(\rho_{j-1}^+)v(\rho_{j-1}^+)].$$ Time discretization is typically achieved via:

• Backward Euler (first order, unconditionally energy-stable) (Zhao et al., 2020, Li et al., 2020, Bao et al., 2021).
• High-order schemes using BDF2 or higher, or predictor–corrector structure, for increased accuracy while preserving energy stability via appropriate SAV or extrapolation corrections (Li et al., 12 Jul 2024, Gan et al., 18 Oct 2025).

For anisotropic and/or strongly anisotropic energies, scheme stabilization is enforced through a minimal stabilizing function $k_0(\theta)$ in the surface energy matrix, guaranteeing positivity and coercivity in the discrete energy law (Bao et al., 2021, Bao et al., 2022, Bao et al., 2023, Zhang et al., 2 Apr 2024, Li et al., 28 Jan 2025). The necessary and sufficient condition for unconditional energy stability is often $3\gamma(\theta) - \gamma(\theta-\pi) \ge 0$ (Li et al., 28 Jan 2025).

The resulting discrete energy dissipation law generally reads

$$W^{m+1} - W^m \le -\tau \|\partial_s\mu^{m+1}\|^2 - \text{(boundary dissipation terms)}.$$

Area or volume conservation is shown to be satisfied either exactly or up to machine precision (for suitably constructed fully implicit schemes), and mesh equidistribution properties arise naturally from the variational structure (Zhao et al., 2020, Bao et al., 2021, Bao et al., 2020).

4. Treatment of Anisotropy and Strong Anisotropy

Handling arbitrary, non-symmetric, and strongly anisotropic surface energies is accommodated by introducing the Cahn–Hoffman vector $\boldsymbol{\xi}$ and forming surface energy matrices (e.g., $G_k(n)$, $Z_k(n)$) that incorporate both $\gamma(\theta)$ and its derivatives. Key forms include: $G_k(n) = \gamma(n)I - n\xi^T + \xi n^T + k(n)nn^T,$ where $\xi = \nabla\gamma(p)|_{p=n}$, and $k(n)\geq k_0(n)$ is the minimal stabilizing function ensuring local energy estimates (Bao et al., 2022, Bao et al., 2023, Zhang et al., 2 Apr 2024, Li et al., 28 Jan 2025).

For strong anisotropy—where the interfacial energy may become nonconvex or the problem ill-posed—regularization via a small Willmore energy term is employed, leading to an equivalent conservative model allowing for energy-stable discretizations (Li et al., 5 Jul 2024). The geometric identities used for recasting the high-order terms into divergence-form weak equations are essential for practical implementation and unconditional stability.

5. Application Areas and Benchmark Results

Energy-stable PFEMs have been successfully applied to simulate:

• Solid-state dewetting of thin films, including the migration of contact points and evolution toward Wulff shapes under various anisotropies (Zhao et al., 2020, Li et al., 2020, Bao et al., 2020, Li et al., 5 Jul 2024, Gan et al., 18 Oct 2025).
• Surface diffusion and anisotropic crystalline flows for both closed and open interfaces, including 3D evolution (Bao et al., 2021, Bao et al., 2022, Bao et al., 2023).
• Willmore flow and curvature-driven geometric flows in planar and three-dimensional settings (spheres, tori, Clifford torus), including variants with spontaneous curvature, open surfaces, and boundary conditions (Bao et al., 24 Jan 2024, Bao et al., 26 Jun 2025, Garcke et al., 30 Jun 2025).
• Complex substrate-coupled problems, such as dewetting on curved substrates with careful mass conservation via substrate arclength parameterization (Bao et al., 1 Oct 2024).

Key metrics, including error in manifold distance, convergence rates, mesh ratio indicators, and exact preservation of energy/area, are routinely documented. For all variants, energy dissipation and, where designed, structure preservation (such as area/volume invariance) are observed to hold up to roundoff error, even during long-time or strongly anisotropic evolution.

6. Practical Considerations: Solver Strategies, Mesh Quality, and Extensions

The discrete systems arising from energy-stable PFEMs are frequently linear (semi-implicit schemes) or weakly nonlinear (fully implicit or midpoint-variant schemes); nonlinear systems are efficiently handled with Newton or hybrid Picard–Newton iterations, requiring only a few iterations per time step (Bao et al., 1 Oct 2024, Bao et al., 2022). Explicit control of tangential velocity or the use of “structure-preserving” discretizations—such as mass-lumped inner products and midpoint normals—ensures sustained mesh equidistribution, obviating remeshing even for long computations (Bao et al., 2020, Bao et al., 2021, Bao et al., 26 Jun 2025).

The approach also extends to high-order time integration (BDF2, BDF3, BDF4) and scalar auxiliary variable (SAV) reformulations, combining accuracy with rigorously provable unconditional energy stability (Li et al., 12 Jul 2024, Gan et al., 18 Oct 2025). The use of minimal stabilizing functions and local energy estimates is now recognized as a sharp and broadly applicable criterion in the field (Zhang et al., 2 Apr 2024, Li et al., 28 Jan 2025, Bao et al., 2023).

7. Impact and Future Directions

The development of energy-stable parametric FEMs systematically resolves the challenge of designing geometric interface evolution algorithms that are both accurate and robust over long times, especially for anisotropic, high-order, or strongly nonlinear flows. This framework provides a general methodological template for other geometric PDEs—including general anisotropic gradient flows, multi-phase and multiphysics problems, and geometric PDEs in higher codimension—by reducing the problem to stabilization via local energy estimates and suitable variational formulations.

Ongoing and future developments focus on robust and efficient extension to:

• Further general anisotropies and singular energies (e.g., crystalline cases).
• Multi-interface and multi-phase coupling with junctions and topology changes.
• Three-dimensional Willmore and biharmonic flows with codimension two features.
• Integration with kinetic or elastic effects, bulk-interface coupling, and direct application to nanostructure evolution and advanced materials design.

The energy-stable PFEM, with its structure-preserving, variationally consistent, and computationally efficient approach, now serves as a reference standard for simulation of geometric evolution in materials science, fluid interfaces, and related areas (Zhao et al., 2020, Li et al., 2020, Bao et al., 2021, Bao et al., 2022, Bao et al., 2023, Zhang et al., 2 Apr 2024, Li et al., 28 Jan 2025, Bao et al., 1 Oct 2024, Bao et al., 26 Jun 2025, Gan et al., 18 Oct 2025).