Stability Condition of the Second-Order SSP-IMEX-RK Method for the Cahn–Hilliard Equation

: Strong-stability-preserving (SSP) implicit–explicit (IMEX) Runge–Kutta (RK) methods for the Cahn–Hilliard (CH) equation with a polynomial double-well free energy density were presented in a previous work, speciﬁcally H. Song’s “Energy SSP-IMEX Runge–Kutta Methods for the Cahn–Hilliard Equation” (2016). A linear convex splitting of the energy for the CH equation with an extra stabilizing term was used and the IMEX technique was combined with the SSP methods. And unconditional strong energy stability was proved only for the ﬁrst-order methods. Here, we use a nonlinear convex splitting of the energy to remove the condition for the convexity of split energies and give a stability condition for the coefﬁcients of the second-order method to preserve the discrete energy dissipation law. Along with a rigorous proof, numerical experiments are presented to demonstrate the accuracy and unconditional strong energy stability of the second-order method.


Introduction
In this paper, we consider the Cahn-Hilliard (CH) equation [1]: where φ is the order parameter, M > 0 is a mobility, F(φ) = 1 4 (φ 2 − 1) 2 is a polynomial double-well free energy density, and > 0 is the gradient energy coefficient. We assume that φ and ∇φ are periodic along the normal to the boundary of a domain Ω ⊂ R d (d = 1, 2, 3). The CH equation has been applied to a wide range of problems [2] and is an H −1 -gradient flow of the following free energy functional: i.e., ∂φ ∂t = M∆ δE δφ , where δ δφ denotes the variational derivative. Thus, E (φ) is nonincreasing in time. The interesting coarsening process of large systems usually occurs on a very long time scale. Therefore, energy stable schemes with high-order time accuracy are highly desirable to perform long time simulations for the coarsening process and there are various related works. Gomez and Hughes introduced a second-order semi-implicit method based on the Crank-Nicolson method [3]. In [4], Guan et al. presented a second-order convex splitting scheme by combining the convex splitting idea [5,6] and the secant method [7]. Yang developed first-and second-order schemes based on the invariant energy quadratization idea [8]. In [9], Shin et al. proposed high-order (up to third-order) convex splitting schemes by combining the convex splitting idea and the specially designed implicit-explicit (IMEX) Runge-Kutta (RK) method. Shen et al. [10] presented second-order backward differentiation and Crank-Nicolson formulas based on the scalar auxiliary variable approach. In [11][12][13], Gong et al. proposed high-order (up to sixth-order) schemes by combining the energy quadratization technique and a specific class of RK methods. Recently, strong-stability-preserving (SSP) IMEX-RK methods for the CH equation were presented [14]. The main idea of the methods was to use a linear convex splitting of E (φ) with an extra stabilizing term (convex and concave parts of E (φ) are treated implicitly and explicitly, respectively) and combine the IMEX technique [15,16] with the SSP methods [17][18][19]. In [14], unconditional strong energy stability (the energy is bounded by its value at the previous time step) was proved only for the first-order time-accurate methods.
In our work, we concentrate mainly on energy stability of the second-order time-accurate (three-stage) SSP-IMEX-RK method for the CH equation. The main issue, which is different from that in [14], consists of two aspects. First, we use a nonlinear convex splitting of E (φ) to remove the condition for the convexity of split energies. Second, we give a stability condition for the coefficients of the second-order method to preserve the discrete energy dissipation law.
This paper is organized as follows. In Section 2, we present the second-order method with the nonlinear convex splitting and prove the method with the stability condition is unconditionally strongly energy stable. In Section 3, we present numerical examples showing the accuracy and energy stability of the method. Finally, conclusions are drawn in Section 4.

Second-Order SSP-IMEX-RK Method and Its Stability Condition
In order to present the second-order SSP-IMEX-RK method for the CH equation, we split E (φ) into convex and concave parts [5,6]: Then, we have the following lemma.

Lemma 1.
The convexity of E c (φ) and E e (φ) yields the following inequality: Proof. We refer to [20].

Definition 1. (Stability Condition). Define a matrix M given by
The stability condition is defined as Next, we show that the scheme (5) with their coefficients satisfying the stability condition is unconditionally strongly energy stable. Theorem 1. The scheme (5) with their coefficients satisfying the stability condition (6) is unconditionally strongly energy stable, i.e., it satisfies for any time step ∆t > 0.
Proof. Using Lemma 1, we have δφ . The second-step of (5) can be rewritten as follows: And the third-step of (5) is (1) ). Then, This completes the proof.
Examples of the coefficients of the scheme (5), satisfying the stability condition (6), are

Numerical Implementation
In order to make order of accuracy in space compatible with second-order in time, we employ the Fourier spectral method [21][22][23] in space for the scheme (5) to arrive at fully discrete second-order SSP-IMEX-RK method. Then, the fully discrete second-order SSP-IMEX-RK method with (6) can be proved similarly to preserve the energy dissipation law in the fully discrete level.
We consider a two-dimensional space Ω = [0, L x ] × [0, L y ] for simplicity and clarity of exposition. One-and three-dimensional cases are defined analogously. Let N x and N y be positive integers and ∆x = L x /N x and ∆y = L y /N y be the space step sizes. In order to solve with the periodic boundary condition, we employ the discrete Fourier transform: for k x = 0, 1, . . . , N x − 1 and k y = 0, 1, . . . , N y − 1, where x l x = l x ∆x, ξ k x = 2πk x /L x , y l y = l y ∆y, and ξ k y = 2πk y /L y . Then, the first-step of (5) can be rewritten in the form where F denotes the discrete Fourier transform and F −1 its inverse transform. The nonlinearity in Equation (11) comes from (φ (1) ) 3 and this can be handled using the truncated Taylor expansion [21][22][23] (φ n,m+1 ) 3 ≈ (φ n,m ) 3 + 3(φ n,m ) 2 (φ n,m+1 − φ n,m ) for m = 0, 1, . . .. We then develop a fixed point iteration method as where φ n,0 = φ n , and we set if a relative l 2 -norm of the consecutive error φ n,m+1 −φ n,m 2 φ n,m 2 is less than a tolerance tol. In this paper, the biconjugate gradient (BICG) method is used to solve the system (12), and we use the following preconditioner P to accelerate the convergence speed of the BICG algorithm: The stopping criterion for the BICG iteration is that the relative residual norm is less than tol.
The second-and third-steps of (5) are implemented analogously.
The computational domain is Ω = [0, 1] × [0, 1]. We set M = 10 −3 , = 10 −2 , and tol = 10 −7 ∆t, and compute φ(x, y, t) for 0 < t ≤ T = 0.4. The coefficients in (7) are used and the grid size is fixed to ∆x = ∆y = 1/256 which provides enough spatial accuracy. In order to estimate the convergence rate with respect to ∆t, simulations are performed by varying ∆t = T/2, T/2 2 , . . . , T/2 6 . Figure 1 shows the relative l 2 -errors of φ(x, y, T) for various time steps. Here, the error is computed by comparison with a reference numerical solution using ∆t = T/2 8 , i.e., is defined as where φ ∆t is a solution with a time step ∆t and φ ref is the reference solution. It is observed that the method is second-order accurate in time.

Energy Stability of the Proposed Method
In order to investigate the energy stability of the proposed method, we take an initial condition as φ(x, y, 0) = rand(x, y).

Conclusions
In this work, we investigated energy stability of the second-order (three-stage) SSP-IMEX-RK method for the CH equation with the polynomial double-well free energy density F(φ) = 1 4 (φ 2 − 1) 2 . Under the stability condition, unconditional strong energy stability of the second-order method was proved theoretically. And, since the nonlinear convex splitting of E (φ) was used compared to the linear convex splitting presented in [14], the restriction for the convexity of split energies was removed. We carried out numerical experiments to verify the accuracy and energy stability of the method.
We note that a choice of coefficients of the method may affect a convergence constant. An optimal choice of coefficients in terms of accuracy can be considered as the scope of future research. And we also note that there is a difficulty associated with the singularity as φ approaches −1 or 1 for a logarithmic free energy density F(φ) = θ 2 (1 + φ) ln 1+φ 2 + (1 − φ) ln 1−φ 2 + θ c 2 (1 − φ 2 ). In order to avoid this, one can consider to apply a regularization to the logarithmic function [24,25] or develop a positivity-preserving scheme [26]. An extension of our method to the case of logarithmic free energy density using such approaches requires further investigation.