A second order accurate scalar auxiliary variable (SAV) numerical method for the square phase field crystal equation

Motivation

• The SPFC model replaces the standard PFC model’s φ^4 energy term with |∇φ|^4, and its analysis and numerical approximation are more challenging.The introduction notes that relatively few numerical results for SPFC were available.Source: §1 · p. 2

Energy decomposition

• A direct SAV application fails because the concave diffusion energy −∥∇φ∥² lacks a global lower bound, alone or combined with the quadratic convex terms.The proposed split places this concave term in the nonlinear energy, where its combination with the 4-Laplacian energy has lower bound −|Ω|.Source: §1 · p. 3

Numerical method

• The energy split yields a second-order BDF2 SAV scheme whose linear diffusion terms have positive eigenvalues and whose steps require only a few constant-coefficient Poisson-like FFT solvers.The nonlinear chemical potential and nonlinear energy value use explicit extrapolation.Source: §1 · p. 3

Analysis

• The scheme has an optimal O(∆t^2 + h^m) convergence rate in the ℓ∞(0,T; H²_N) energy norm.The error analysis uses discrete H² and H³ bounds and aliasing-error control for the 4-Laplacian term.Source: §1 · pp. 3–4

2.1 Review of Fourier pseudo-spectral approximations

• Fourier pseudo-spectral approximation represents functions on a quadrature grid, enabling efficient operator evaluation with FFT-based algorithms.The section also introduces discrete Fourier operators and estimates used to control collocation interpolation and nonlinear terms.Source: §2.1 · pp. 4, 6, 7

2.2 The fully discrete numerical scheme

• The fully discrete scheme uses an SAV reformulation and second-order BDF time stepping, with explicit extrapolation to approximate nonlinear terms.A ghost-point approximation initializes the two-step method with O(∆t^2 + h^m) accuracy.Source: §2.2 · pp. 7–8

2.3 Unique solvability and efficient numerical solver for the proposed scheme

• For given prior phase fields and auxiliary-variable values, the scheme has a unique solution; with constant initial data, it conserves that mass at every step.The solvability argument reduces to a scalar equation with a positive coefficient.Source: §2.3 · p. 9

3.1 Modified energy stability for the proposed numerical scheme

• The numerical solution satisfies dissipation properties for the proposed scheme’s modified energy.The energy law is derived by combining the BDF temporal terms for the phase variable and auxiliary scalar.Source: §3.1 · p. 10
• Under sufficiently regular initial data, energy stability yields a uniform-in-time discrete H2_N bound with a constant independent of h, Δt, and time.The stability estimate is for the modified energy, which is a second-order approximation to the original discrete energy under regularity assumptions.Source: §3.1 · p. 11
• A uniform-in-time bound for the original energy is also established, but it is less sharp than the primitive-variable estimate.The comparison attributes the weaker bound to preservation of dissipation for the reformulated energy involving the auxiliary variable.Source: §3.1 · pp. 12–13

3.2 The ℓ∞(0, T; H3) bound estimate for the numerical solution

• The scheme admits a uniform-in-time H3 estimate whose constant depends on the domain and initial H3 data, but not Δt, h, or T.The estimate uses the uniform H2 bound and discrete Sobolev and aliasing-error controls.Source: §3.2 · pp. 13, 14, 15, 17

4.1 The consistency analysis

• The consistency analysis gives a truncation error bounded by C(∆t^2 + h^m) using the extrapolated phase variable and spectral approximation estimates.The derivation of the truncation estimate is omitted from the paper excerpt.Source: §4.1 · p. 18

4.2 A few preliminary estimates

• Preliminary estimates control nonlinear energy and gradient-error terms using exact-solution regularity, discrete W^1,6 bounds, and Sobolev embeddings.These bounds are used in the convergence proof to estimate the nonlinear terms.Source: §4.2 · pp. 19, 21

4.3 Proof of the convergence theorem

• The convergence proof cancels the first nonlinear error inner-product terms and applies a discrete Gronwall inequality to obtain the error estimate.The cancellation is identified as crucial to the analysis.Source: §4.3 · pp. 25–26
• The scheme has an optimal-rate convergence estimate in ℓ∞(0, T; H2_N) with no scaling-law constraint between ∆t and h, provided ∆t and h are sufficiently small and the exact solution has regularity class R.The estimate constant is independent of ∆t and h; the proof uses a higher-order H3 bound for the numerical solution.Source: §4 · p. 18, §4.3 · pp. 26–27
• For sufficiently small ∆t and h, the refined estimate is sharper than the uniform-in-time bound, but is local in time because its convergence constant depends on the final time.The rough bound (3.20), by contrast, is global.Source: §4.3 · p. 27
• The SAV scheme requires a linear solve at each step, while the modified BDF primitive-formulation scheme in requires a nonlinear solve costing approximately three times as much at the same resolutions.The comparison attributes the nonlinear solve to implicit treatment of the 4-Laplacian term for unconditional energy stability.Source: §4.3 · p. 27
• The analysis adds technical tools for higher-order stability and convergence; its uniform-in-time H3_N bound (3.26) is not needed in.The passages contrast the article’s ℓ∞(0,T;H2_N) error estimate with the estimate in, but the excerpt does not fully preserve the norm label for the latter.Source: §4.3 · p. 27

5.1 Convergence test for the numerical scheme

• The discrete ℓ2 and ℓ∞ errors show an approximate slope of -2 as temporal resolution increases, confirming second-order temporal convergence.The test fixes N = 128 and uses NT from 100 to 1000.Source: Figure 1 · p. 28, §5.1 · p. 28

5.2 Numerical simulation of square symmetry patterns

• The simulations show emergence of rotationally invariant square-symmetry patterns from random initial data on a 100-by-100 domain.The time step increases from 0.01 to 0.02 after t = 1000.Source: §5.2 · p. 28
• The original and SAV-introduced energy curves overlap, and energy dissipation is observed through t = 1000.The comparison is shown for the one-nucleation-site example.Source: §5.2 · pp. 28–29
• For a magnitude-10 singular nucleation at (50, 50), direct SAV does not produce a reasonable solution, while stabilized SAV can handle the rough initial data.The reported limitation contrasts with the scheme's performance for smooth flows with mild random perturbations.Source: §5.2 · p. 29

• The proposed SAV scheme uses an energy decomposition whose nonlinear part has a global lower bound and whose remaining terms have positive-eigenvalue constant-coefficient diffusion.The decomposition addresses a key difficulty in applying SAV to the SPFC model.Source: §6 · p. 29
• The scheme can be implemented efficiently with constant-coefficient Poisson-like solvers via FFT, and modified-energy stability yields a uniform-in-time H2 bound.Source: §6 · p. 29
• A uniform-in-time H3 bound supports the optimal-rate error estimate in ℓ∞(0, T; H2) ∩ ℓ2(0, T; H5), and experiments demonstrate accuracy and efficiency, including square-pattern simulations.Source: §6 · p. 29

• A grid function f is extended to a continuous trigonometric polynomial, with the space specified by a maximum degree.Source: A Proof of Proposition 2.4 · pp. 29–30
• Parseval's identity and elliptic regularity for the continuous extension are used to obtain estimates for the discrete setting.Source: A Proof of Proposition 2.4 · p. 30
• The discrete elliptic regularity inequality is obtained by combining the two preceding estimates.Source: A Proof of Proposition 2.4 · p. 30