Convergence Analysis of a Linear, Unconditionally Energy-Stable SAV Finite Element Method for the Cahn-Hilliard Equation

1. Problem Statement

The paper addresses the numerical simulation of the Cahn-Hilliard (CH) equation, a fourth-order nonlinear partial differential equation widely used to model phase separation and coarsening dynamics in binary mixtures.

• Challenges: The CH equation presents two primary numerical difficulties:
  1. Nonlinear Stiffness: The double-well potential term ($F'(\phi)$) imposes strict restrictions on time-step sizes for explicit schemes.
  2. High-Order Derivatives: The fourth-order nature requires high regularity in spatial discretization.
• Existing Limitations: Traditional fully implicit schemes are stable but computationally expensive due to the need to solve nonlinear systems at every time step. Explicit or semi-implicit schemes are cheaper but often lack unconditional energy stability or require very small time steps.
• Gap: While the Scalar Auxiliary Variable (SAV) method has been successful in spectral frameworks, a rigorous optimal-order convergence analysis for a fully discrete finite element method (FEM) combining backward Euler time stepping and linear spatial elements was lacking.

2. Methodology

The authors propose a fully discrete finite element scheme based on the Scalar Auxiliary Variable (SAV) approach.

A. Mathematical Reformulation

The original CH system is reformulated by introducing a scalar auxiliary variable $r(t) = \sqrt{E(\phi) + C_0}$, where $E(\phi)$ is the free energy. This transforms the nonlinear potential term into a quadratic form involving $r$, allowing the system to be written as:
$$\phi_t = \Delta \mu$$
$$\mu = -\Delta \phi + \frac{r}{\sqrt{E(\phi)}} F'(\phi)$$
$$r_t = \frac{1}{2\sqrt{E(\phi)}} \int_\Omega F'(\phi) \phi_t \, dx$$
This reformulation decouples the nonlinearity from the linear operators, enabling the construction of a linear scheme.

B. Discretization Scheme

• Temporal Discretization: Backward Euler method (first-order in time).
• Spatial Discretization: Linear finite elements ($P_1$) on a quasi-uniform triangular mesh.
• Algorithm Design: To avoid solving a large coupled linear system, the authors employ a decoupling strategy. The variables $(\phi_h, \mu_h)$ are expressed as a linear combination of two auxiliary solutions $(\phi_{1h}, \mu_{1h})$ and $(\phi_{2h}, \mu_{2h})$. This reduces the problem to solving two linear systems with constant coefficient matrices and a single scalar equation for $r_h$ at each time step.

3. Key Contributions

The paper makes three primary theoretical and practical contributions:

1. Unconditional Energy Stability:
   The authors prove that the proposed scheme is unconditionally energy stable. They derive a discrete energy inequality showing that the modified total free energy decreases monotonically at each time step, regardless of the time step size ($\Delta t$) or mesh size ($h$).
   $$E^{n+1}_{total} - E^n_{total} + \text{positive terms} = -\|\nabla \mu^{n+1}_h\|^2 \Delta t$$

2. Rigorous Convergence Analysis:
   The paper provides a complete a priori error analysis for the fully discrete scheme. By introducing a specific elliptic projection operator, the authors derive optimal-order error estimates in the $H^1$-norm for:

   • The phase variable ($\phi$).
   • The chemical potential ($\mu$).
   • The auxiliary variable ($r$).

3. Optimal Convergence Rates:
   Under suitable regularity assumptions, the scheme is proven to achieve:

   • First-order convergence in time ($O(\Delta t)$).
   • First-order convergence in space for $\phi$ and $\mu$ ($O(h)$).
   • Second-order convergence in space for the auxiliary variable $r$ ($O(h^2)$).

4. Results

The theoretical findings are validated through numerical experiments on the domain $\Omega = [0, 1] \times [0, 1]$.

• Energy Stability Verification: Simulations demonstrate that both the original free energy and the SAV-modified energy decay monotonically over time, confirming the unconditional stability of the scheme. Mass conservation is also verified.
• Phase Separation Dynamics: The numerical solution successfully captures the physical process of phase separation, showing the evolution from a homogeneous state to the formation of distinct phases and subsequent coarsening.
• Convergence Rates:
  • Time Convergence: With a fixed fine mesh, the errors for $\phi$, $\mu$, and $r$ decrease linearly as $\Delta t$ is reduced, confirming $O(\Delta t)$ accuracy.
  • Space Convergence: With a fixed small time step, the errors for $\phi$ and $\mu$ decrease linearly with mesh refinement ($O(h)$), while the error for $r$ decreases quadratically ($O(h^2)$).
  • These results align perfectly with the theoretical bounds established in Theorem 3.3.

5. Significance

This work is significant for the computational mathematics community for several reasons:

• Bridging Theory and Practice: It fills a critical gap by providing a rigorous convergence proof for SAV methods within the finite element framework, moving beyond spectral methods.
• Efficiency and Stability: The proposed scheme offers a "best of both worlds" solution: it retains the unconditional energy stability required for long-time simulations while maintaining computational efficiency by solving only linear systems with constant coefficients.
• Robustness: The decoupling algorithm makes the method highly scalable and easy to implement for complex phase-field models coupled with other physics (e.g., fluid flow), as the linear systems can be solved using standard, efficient solvers.
• Foundation for Future Work: The established error analysis framework provides a solid theoretical basis for applying SAV methods to more complex, coupled phase-field systems where rigorous convergence proofs are currently missing.