Novel, linear, decoupled and unconditionally energy stable numerical methods for the coupled Cahn–Hilliard equations

This paper uses a novel numerical approach to approximate the coupled Cahn–Hilliard equations, which are a highly nonlinear system depicting the phase separation of the homopolymer and copolymer mixtures. The new method is named 3S-IEQ, and its construction and calculation are more straightforward than the invariant energy quadratization and scalar auxiliary variable methods. Notably, we only need to solve two linear decoupled constant-coefficient equations at each time step. Numerical simulations are shown


Introduction
The phase field model is a mathematical method used to solve interface problem. It is mainly applied to solidification dynamics, but it has also been applied to other situations, such as viscous fingering, crack dynamics, vesicle dynamics, etc. In this method, the boundary conditions at the interface are replaced with partial differential equations to obtain the evolution of the auxiliary field (phase field) as an ordered parameter. The typical phase field models include Allen-Cahn equation, Cahn-Hilliard equation, phase field crystal equation, and Ohta-Kawasaki equation [1][2][3][4]. There are many works devoted to approximating phase field models, such as the nonlinear convex splitting method, the linear stabilized semi-implicit method, the invariant energy quadratization (IEQ) method, and the scalar auxiliary variable (SAV) method [5][6][7][8].
The coupled Cahn-Hilliard equations are a highly nonlinear system depicting the phase separation of the homopolymer and copolymer mixtures [9]. Avalos [10] et al. have introduced two phase variables to describe the macro phase separation between homopolymer and copolymer as well as micro phase separation between two components of diblock copolymers, respectively. The authors [11] have designed two efficient, decoupled, and second-order unconditionally energy stable numerical schemes for this system by utilizing the SAV method.
In this paper, we apply the 3S-IEQ method, which stands for a step-by-step solving approach based on IEQ and constructed by Liu and Chen in [12], to design three linear, decoupled, and unconditionally energy stable numerical schemes for the coupled Cahn-Hilliard system. The IEQ method sometimes requires specific matrix solvers, such as iteratively solving the elliptic equations with complex variable coefficients at each time step. Although the SAV method does not suffer such a requirement, it is still necessary to solve four fourth-order linear equations at each time step, while we only solve two fourth-order linear equations by 3S-IEQ in terms of the coupled Cahn-Hilliard system.

Governing equations
We first give some notations, which will be used later. For every k ≥ 0, denote (·, ·) k and · k to be the H k ( ) inner product and norm, separately. In particular, we abbreviate the L 2 inner product (·, ·) 0 and norm · 0 as (·, ·) and · , separately. We further define some Sobolev spaces is the unique solution of the following elliptic problem: Consequently, we define ψ f := (-) -1 f . Assuming f , g ∈ L 2 0 ( ), the H -1 per ( ) inner product and norm can be written as Integrating by parts with periodic boundary condition, one has Next, we provide a brief description of the coupled Cahn-Hilliard equations proposed in [9,10], which depict the phase transition of the mixture of a homopolymer and a copolymer. Consider the following Swift-Hohenberg type free energy function: where is a smooth, open, bounded, connected domain in R d (d = 1, 2, 3), ε u and ε v are two real coefficients, proportional to the thickness of the interface between macro phase u and micro phase v. A(u, v) = 1 4 (u 2 -1) 2 + 1 4 (v 2 -1) 2 + αuv + βuv 2 , α and β are two coupled parameters. As a consequence, the coupled Cahn-Hilliard equations can be derived by taking the variational derivative of (2.2) in H -1 ( ) with respect to u and v, separately.
where M u and M v are the motion parameters controlling the moving speed, f (u, v) = u 3u + αv + βv 2 and g(u, v) = v 3v + αu + 2βuv. Throughout this paper, we assume that all variables satisfy the periodic boundary condition on ∂ . Suppose that there is no external forcing other than gravity, and the coupled Cahn-Hilliard equations satisfy the following energy dissipation law and mass conservation law: Before providing the discrete formulations, we let N > 0 be a positive integer. For n ≤ N , set δt = T/N , t n = nδt, u n denotes the numerical approximation of u(t n ) and u n+ 1 2 = u n+1 +u n 2 .

Numerical schemes
In this section, we apply the 3S-IEQ method, which was first introduced by Liu and Chen in [12], to design linear, decoupled, and unconditionally energy stable numerical schemes for the coupled Cahn-Hilliard system. Both the first-order backward Euler (BDF1), the second-order Crank-Nicolson (CN), and the second-order backward differentiation formula (BDF2) methods can be successfully applied to this system. Without losing generality, we mainly focus on the second-order CN scheme. We introduce an auxiliary variable η as follows: Then, we define the following new functions H(u, v, η) and G(u, v, η) which are equivalent to the nonlinear functional f (u, v) and g (u, v).
Thus, the original coupled Cahn-Hilliard system (2.3) can be rewritten as By taking the L 2 inner product of (3.3a) and (3.3b) with (-) -1 u t and (-) -1 v t , we can obtain the following modified energy dissipation law:

Theorem 3.2 Scheme (3.15a)-(3.15c) is unconditionally energy stable satisfying the following discrete energy dissipation law:
Proof The proof of this theorem is similar to Theorem 3.1. Indeed, we rely on the identity So we omit the proof.

Theorem 3.3
The second-order Scheme 1 and Scheme 2 are uniquely solvable. (u n+1 , v n+1 , η n+1 ) can be solved step by step, the fast Fourier transform is highly efficient for solving the fully discrete schemes.
Proof Taking the second-order CN Scheme 1 as an example, the proof of BDF2 Scheme 2 is similar. We rewrite Scheme 1 as the following formulations: (3.18) The two linear systems are uniquely solvable because the coefficient matrices are not symmetric positive definite but constant coefficient matrices. Meanwhile, u n+1 , v n+1 , and η n+1 can be solved totally decoupled. We can first calculate u n+1 , v n+1 , then η n+1 can be updated by Remark 3.1 In [11], the authors proposed two time second-order schemes based on the CN and the BDF2 methods. Although the second-order schemes derived by SAV are de-coupled and unconditionally energy stable, the authors need to solve four fourth-order linear equations with constant coefficients at each time step. Here we only solve two fourthorder linear equations that significantly improve the efficiency.
Remark 3.2 Both the second-order Scheme 1 and Scheme 2 include three time levels, and (u n+1 , v n+1 , η n+1 ) can be updated recursively by the recurrence relation as long as the initial values (u 0 , v 0 , η 0 ) and (u 1 , v 1 , η 1 ) have been calculated. The former is directly derived by initial conditions, and the latter is computed by BDF1.

Numerical experiments
We now present various two-dimensional numerical simulations for the coupled Cahn-Hilliard equations. Usually, the proposed time discretization schemes can be combined with any spatial discretization, provided the spatial discretization reaches the desired tolerance error. In this paper, we assume that all variables satisfy the periodic boundary condition and merely focus on time discretization. Therefore, we can apply the Fourier spectral method for spatial discretization, and the fast Fourier transform is easy to solve the fully discrete schemes.

Temporal convergence test
We test the convergence rates of the three proposed schemes. The parameters are The computational domain is = [-3 * π, π] 2 , and we use 128 × 128 Fourier modes. In Table 1 and Table 2, we list the L 2 errors for u and v at final time T = 0.05 with different time step sizes δt = 0.005/2 k , k = 0, 1, . . . , 5. These exact solutions and parameters are completely consistent with those in [11] for convenient comparison. From Table 1 and Table 2, we can observe that all schemes reach the desired order of accuracy in time. BDF1 is first-order accurate in time, while CN and BDF2 are secondorder accurate in time. Compared with the results of [11], the error is a little bit bigger because 3S-IEQ is a total explicit method.

Energy stability test and phase transition
We select the following smooth initial conditions to demonstrate the energy stability of the proposed schemes: u(x, y, 0) = sin 2x(x -1)y(y -1) , v(x, y, 0) = cos 10x(x -1) y(y -1), (4.2) where = [0, 1] 2 , and we use 128 × 128 Fourier modes for spatial discretization. From Fig. 1, the original energy E, the modified energy E, the discrete energies E bdf 1 , E cn , and E bdf 2 are nonincreasing in time. In Fig. 2, we plot the evolution curves of modified discrete energies E bdf 1 , E cn , and E bdf 2 , which are calculated by BDF1, CN, and BDF2 for various time steps. Overall, the schemes developed in this paper are unconditionally energy stable. For the same time step, the discrete energy of CN scheme can obtain a better approximation, especially δt = 0.001, 0.0005, 0.0001, which is consistent with the results in [11]. Moreover, Fig. 2 provides a basis for us to select an appropriate time step when simulating the longterm phase transition process, that is, δt ≤ 0.001. Figure 3 displays the snapshots of the phase variables u and v at different times. Moreover, we can see that u describes a macro phase separation between homopolymer and copolymer, while v describes a micro phase separation that occurs within the separate domain. Significantly, we choose the same parameters and initial data as [11] to demonstrate the consistent snapshots of u and v.
In Table 3, we give the CPU time consumed by 3S-IEQ compared with the SAV method in [11]. The three schemes of the former method take significantly less time than the latter, especially for the minimal time step simulation δt = 0.0001, because the authors in [11]  needed to solve four fourth-order linear equations with constant coefficient at each time step. Here we only solve two fourth-order linear equations, which greatly improves the efficiency.

Conclusions
This paper applies the 3S-IEQ method to design three linear, decoupled, and unconditionally energy stable numerical schemes for the coupled Cahn-Hilliard system. Numer-ical simulations are shown to demonstrate the time accuracy and unconditional energy stability. The most significant advantage of this method is that we only solve two fourthorder linear equations, which dramatically improves the efficiency, and it can be applied to various gradient flows, including L 2 and H -1 gradient flows. However, the error and convergence analysis of the proposed schemes is still an open challenge due to the coupling and nonlocality. We will carefully consider the bottleneck in the future.