A High-Order Convex Splitting Method for a Non-Additive Cahn–Hilliard Energy Functional

Various Cahn–Hilliard (CH) energy functionals have been introduced to model phase separation in multi-component system. Mathematically consistent models have highly nonlinear terms linked together, thus it is not well-known how to split this type of energy. In this paper, we propose a new convex splitting and a constrained Convex Splitting (cCS) scheme based on the splitting. We show analytically that the cCS scheme is mass conserving and satisfies the partition of unity constraint at the next time level. It is uniquely solvable and energy stable. Furthermore, we combine the convex splitting with the specially designed implicit–explicit Runge–Kutta method to develop a high-order (up to third-order) cCS scheme for the multi-component CH system. We also show analytically that the high-order cCS scheme is unconditionally energy stable. Numerical experiments with ternary and quaternary systems are presented, demonstrating the accuracy, energy stability, and capability of the proposed high-order cCS scheme.


Introduction
The Cahn-Hilliard (CH) equation was originally introduced as a phenomenological model of phase separation in a binary alloy [1] and has been applied to a wide range of problems [2]. The CH equation is derived from the Ginzburg-Landau energy functional: where c is the concentration field defined in Ω (⊂ R d , d=1, 2, 3) and w, > 0 are the free energy and gradient energy coefficients. The CH equation is a gradient flow for E GL (c) in the H −1 -inner product, thus E GL (c) is nonincreasing in time.
Generalizations of E GL (c) for more than two components can be applied to wide range of problems, thus have been studied intensively [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. A great deal of research has been focused on the ternary system [10] such as: where c i is the concentration field of the phase i. There are many other forms of energy functional for the ternary system and some of them are equivalent. For example, under the constraint c 1 + c 2 + c 3 = 1, the first term in E FP (c 1 , c 2 , c 3 ) can be rewritten [11] as follows: One of the most important criteria for the multi-component model is to avoid the generation of spurious phases. To be more precise, we think a physically reasonable model must satisfy the following two fundamental criteria: (A) (Consistency of null-phase) If a phase is absent at the initial time, it should not appear at any time.

(B) (No additional phase on interface)
The interface including multiple junction should be free of additional phases.
The model with E FP (c 1 , c 2 , c 3 ) and many other models with three components obey both criteria; however, it is not well-known how to construct an energy functional for more than three components satisfying mathematical and physical criteria including these two. For example, the following generalization introduced by Lee and Kim [12] for the vector-valued concentration field c = (c 1 , . . . , c N ) T : does not satisfy Criteria A and B when N > 3.
In this paper, we consider the following energy functional satisfying Criteria A and B: where The energy functional E (c) introduced by Tóth et al. [14] is a non-trivial extension of E GL (c) for more than three component system in the sense that it is equivalent to E FP (c 1 , c 2 , c 3 ) when N = 3 and E GL (c) when N < 3. We develop a high-order energy stable numerical method for this energy with quadratically mixed terms ∑ i<j c 2 i c 2 j to study phase separation in multi-component systems. The H −1 -gradient flow for E (c) is given by under the partition of unity constraint, where µ = (µ 1 , . . . , µ N ) T is the vector-valued chemical potential, µ i is the chemical potential of the phase i, δ δc denotes the variational derivative with respect to c, ∂F ∂c = ∂F is a Lagrange multiplier to ensure the constraint [8,[10][11][12][13]16,17], and 1 = (1, . . . , 1) T ∈ R N . We consider the boundary conditions for c i and µ i as the zero Neumann boundary conditions: where n is a unit normal vector to ∂Ω. We refer to Equation (2) as the vector-valued CH (vCH) equation. Because the vCH equation is of gradient type, E (c) is nonincreasing in time as the constraint holds: The vCH equation is a fourth-order nonlinear partial differential equation and the N unknowns c 1 , . . . , c N are linked through the constraint. Therefore, accurate and efficient numerical methods are desirable to study the dynamics of the vCH equation. In this paper, we propose a constrained Convex Splitting (cCS) scheme for the vCH equation, which is based on a convex splitting of E (c) under the constraint. For N = 2 and 3, E (c) has a straightforward convex-concave splitting. However, there is a difficulty with N > 3 since Ω wF(c) dx in E (c) is neither convex nor concave. To apply the convex splitting idea [18][19][20] for all N, we add and subtract an auxiliary term in E (c). Then, a convex-concave decomposition is available. We show analytically that the cCS scheme is mass conserving and satisfies the constraint at the next time level. It is uniquely solvable and energy stable. Furthermore, we combine the convex splitting with the implicit-explicit Runge-Kutta (RK) method [21,22] to develop a high-order (up to third-order) cCS scheme. We employ the specially designed implicit-explicit RK tables [23] to have both high-order time accuracy and unconditional energy stability. We also show analytically that the high-order cCS scheme is unconditionally energy stable.
This paper is organized as follows. In Section 2, we describe the convex splitting with an auxiliary term. We propose the (first-order) cCS scheme for the vCH equation and prove its unconditional unique solvability and energy stability. In Section 3, we construct the high-order cCS scheme with a proof of unconditional energy stability. In Section 4, we present numerical examples showing the accuracy, energy stability, and capability of the proposed scheme. Finally, conclusions are drawn in Section 5.

Constrained Convex Splitting Scheme
The convex splitting of E (c) when N < 3 is trivial; however, it is not well-known how to split E (c) for N ≥ 3 into convex and concave parts. Therefore, we here propose to split E (c) for N ≥ 3 where .
From the positive semi-definiteness of H(F c (c)), we obtain , Thus, the convexity of E c (φ) and E e (φ) is assured.
We now present the cCS scheme for the vCH Equation (2) by treating E c (c) implicitly and E e (c) explicitly under the partition of unity constraint in Equation (3): (5) is mass conserving.

Lemma 2. The cCS scheme in Equation
Proof. Let c n+1 be a solution of the cCS scheme. From Equation (5), we have where we used the zero Neumann boundary condition for µ i . It follows that Ω c n+1 i dx = Ω c n i dx.
Lemma 3. The cCS scheme satisfies the constraint at any time t n , i.e., where I denotes the identity operator. Since I + ∆t 2 ∆ 2 with the zero Neumann boundary condition for c i is an invertible operator, Equation (6) has a unique solution. Thus, Equation (6) Theorem 1. The cCS scheme with an initial condition satisfying ∑ N i=1 c 0 i = 1 is uniquely solvable for any time step ∆t > 0.

Proof. We consider the following functional on
where (·, ·) H −1 and (·, ·) L 2 denote the H −1 -and L 2 -inner products, respectively, and (−∆c, d) [24] for the definition of the H −1 -inner product). It may be shown that c n+1 ∈ H is the unique minimizer of G(c) if and only if it solves, for any d ∈ H 0 , because G(c) is strictly convex by the convexity of E c (c). Equation (7) is true for any d ∈ H 0 if and only if Equation (5) holds. Hence, minimizing the strictly convex functional G(c) is equivalent to solving Equation (5).

Lemma 4.
The convexity of E c (c) and E e (c) yields the following inequality: Proof. Since both E c (c) and E e (c) are convex, we obtain Theorem 2. The cCS scheme with an initial condition satisfying ∑ N i=1 c 0 i = 1 is unconditionally energy stable, meaning that, for any ∆t > 0, Proof. Setting c = c n+1 and d = c n in Equation (8), we have

Extension of the Constrained Convex Splitting Scheme To High-Order Time Accuracy
The cCS scheme in Equation (5), is first-order accurate in time and its order of time accuracy can be improved by various approaches. One of them is to combine with an s-stage implicit-explicit RK method [21]: let c (0) = c n , where a kl andâ kl are RK coefficients for k = 1, . . . , s, and then c n+1 = c (s) by using the stiffly accurate condition.

Proof.
The proofs are similar to Lemmas 2 and 3 and Theorem 1, thus we omit the details here.
Before proving the energy stability of the s-stage high-order cCS scheme, we define an s × s matrix R as R kl = r kl for l ≤ k and R kl = 0 for l > k, and an s × s matrix R as R kl =r kl = r kl − r k−1,l with r 0l = 0. Theorem 3. Suppose that R is positive definite. The s-stage high-order cCS scheme with r kk ≥ 0 for k = 1, . . . , s and an initial condition satisfying ∑ N i=1 c 0 i = 1 is unconditionally energy stable, meaning that for any ∆t > 0, Proof. The analogous proof can be found in [23]. Using Lemma 4, we have where the last equality follows from the fact that Let ∇µ = (∇µ (1) , . . . , ∇µ (s) ) T . Since R is positive definite, (1) , . . . , φ (s) ) T and ψ = (ψ (1) , . . . , ψ (s) ) T . It follows that E (c n+1 ) ≤ E (c n ). Remark 1. The first-order cCS scheme can be viewed as the one-stage cCS scheme with R = ( 1 ) and R = ( 1 ).

Numerical Experiments
The s-stage high-order cCS scheme in Equation (10) can be rewritten as follows: for k = 1, . . . , s, where S (k) = r kk ∆ −w ∂F e (c (k−1) ) ∂c + α e (c (k−1) )1 + ∑ k−1 l=1 r kl ∆µ (l) . The nonlinearity of the scheme comes from ∂F c (c (k) ) ∂c i and α c (c (k) ) and these can be handled using a Newton-type linearization [24][25][26][27]: for m = 0, 1, . . ., . We then develop a Newton-type fixed point iteration method for the scheme as where c (k−1,0) = c (k−1) , for i = 1, . . . , N, and we set if a relative l 2 -norm of the consecutive error is less than a tolerance tol. In this paper, the biconjugate gradient (BICG) method is used to solve the system in Equation (11) 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. For first-, second-, and third-order accuracy, we use the following matrices R, respectively [23]: and The positive definiteness of R is easily seen by showing eigenvalues of 1 2 ( R + R T ) are all positive. The eigenvalues of 1 2 2 3 , and 2 3 + √ 26 8 for Equation (13), and approximately 0.0063, 0.1105, 0.3582, 0.5722, 0.9225, and 1.0303 for Equation (14).
We used the Fourier spectral method for the spatial discretization and the discrete cosine transform in MATLAB was applied for the whole numerical simulations to solve the vCH equation with the zero Neumann boundary condition.

Convergence Test
We demonstrate the convergence of the proposed schemes with the initial conditions on Ω = [0, 2π]. We set = 0.25 and compute c(x, t) for 0 < t ≤ 280. The grid size is fixed to h = 2π/64, which provides enough spatial accuracy. To estimate the convergence rate with respect to ∆t, simulations are performed by varying ∆t = 2 −10 , 2 −9 , . . . , 2 2 . We take the quadruply over-resolved numerical solution using the third-order scheme as the reference solution. Figure 1a,b shows the evolution of E (t) for the reference solution and the relative l 2 -errors of c(x, 120) (this time is indicated by a dashed line in Figure 1a) for various time steps, respectively. It is observed that the schemes give desired order of accuracy in time.

Energy Stability of the Proposed Schemes
To investigate the energy stability of the proposed schemes, we consider the phase separation of a ternary system with the initial conditions Here, rand(x, y) is a random number between −0.1 and 0.1, and we use = 0.1 and h = 2π/64. Figure 2 shows the evolution of E (t) using the first-, second-, and third-order schemes with different time steps. All energy curves are nonincreasing in time even for sufficiently large time steps. Figure 3 shows the evolution of c(x, y, t) using the third-order scheme with ∆t = 2 −2 . In each snapshots, the red, green, and blue regions indicate c 1 , c 2 , and c 3 , respectively, and contour lines represent c i = 0.5.

Consistency of Null-Phase
To confirm whether consistency of null-phase guarantees, we consider that only three phases are present but the simulation is performed using a quaternary system, i.e., we take the initial conditions as . For E LK (c), we employ the convex splitting in [17] and also apply the third-order scheme in Equation (10) with R in (14). We use = 0.25, h = 2π/64, and the third-order scheme and compute c(x, t) for 0 < t ≤ 120. Figure 4a,b shows c 4 (x, 120) obtained with ∆t = 2 −3 and max x |c 4 (x, 120)| for various time steps, respectively, for two models E LK (c) and E (c). As shown in Figure 4, E LK (c) generates c 4 , even though the initial condition for c 4 is zero. We believe that this generation is not a result of numerical computation but a consequence of the model error not satisfying Criterion A. On the other hand, for E (c), the maximum of c 4 is only controlled by the accuracy of the numerical scheme.

No Additional Phase Generation on Interface
To test whether spurious phase generation takes place on interfaces, we consider the phase separation of a quaternary system with the initial conditions Here, rand(x, y) is a random number between −0.1 and 0.1. For E LK (c), we employ the convex splitting in [17] and also apply the third-order scheme in Equation (10) with R in Equation (14). We use = 0.05, h = 2π/64, and ∆t = 2 2 . Figures 5 and 6 show c 4 (x, y, T=3072) and local maximum of c 4 (x, y, t) for 0 ≤ t ≤ T using the third-order scheme for two models E LK (c) and E (c), respectively. As shown in Figure 5, E LK (c) generates spurious phase at two triple junction points where the first, second, and third components meet. On the other hand, E (c) in Figure 6 suppresses the formation of spurious phases on interfaces almost completely. To quantify the spurious phase generation, we define LocalMax(c i ) at time t as a set of local maxima of c i (x, y, t) in Ω. We observe many local maxima near 0.25 at the beginning of evolution but only two maxima 0 and 1 are expected for the fully separated phases. The model with E (c) shown in Figure 6 gives well separated phases over time, whereas the model with E LK (c) has another local maximum at about 0.043 due to the spurious phases at the junction points. We believe that this spurious phases generation is not a result of numerical computation but a consequence of the model error not satisfying Criterion B.   Figure 6. c 4 (x, y, T) and local maxima of c 4 (x, y, t), 0 ≤ t ≤ T for the model E (c).

Conclusions
In this paper, we consider the multi-component CH system where all phase variables are nonlinearly coupled. To study the dynamics of this system, we propose the high-order energy stable scheme based on the convex splitting idea. To handle the nonconvex, nonconcave term in the energy, we add an auxiliary term, which yields a convex-concave decomposition of the energy. We combine the convex splitting with the specially designed implicit-explicit RK method thus developed the high-order cCS scheme. We confirmed that the schemes give desired order of accuracy in time and are unconditionally energy stable. By using the scheme, we also demonstrated that the use of E (c) is crucial and gives very significant qualitative improvement of the results compared to the additive model E LK (c).