Unconditionally energy stable second-order numerical scheme for the Allen–Cahn equation with a high-order polynomial free energy

In this article, we consider a temporally second-order unconditionally energy stable computational method for the Allen–Cahn (AC) equation with a high-order polynomial free energy potential. By modifying the nonlinear parts in the governing equation, we have a linear convex splitting scheme of the energy for the high-order AC equation. In addition, by combining the linear convex splitting with a strong-stability-preserving implicit–explicit Runge–Kutta (RK) method, the proposed method is linear, temporally second-order accurate, and unconditionally energy stable. Computational tests are performed to demonstrate that the proposed method is accurate, efficient, and energy stable.


Introduction
Phase-field equations have arisen as a significant numerical framework in modeling and studying evolution of pattern formation in materials [1]. In general, phase-field models are driven by gradient flows for the governing total energy functionals [2,3], e.g., the Ginzburg-Landau (GL) free energy functional: where is a domain in R d (d = 1, 2, 3), φ is the order parameter, F(φ) = 1 4 (φ 2 -1) 2 is the GL potential, and > 0 is a constant. We assume the homogeneous Neumann boundary condition for φ: ∇φ · n = 0 on ∂ , where n is normal to ∂ . The L 2 -gradient flow for (1) is the Allen-Cahn (AC) equation [4].
Recently, the AC equation with a high-order polynomial F p (φ) = 1 4 (φ p -1) 2 , referred to as the hAC equation, was introduced to better represent the interfacial dynamics [24]. Here, p is an even integer.
For more background information about the hAC equation, let us consider the origin of the standard quartic polynomial free energy, F(φ) = 1 4 (φ 2 -1) 2 . In the original derivation of the bulk free energy in the total free energy functional E(φ), the following logarithmic free energy was used: where θ and θ c are the absolute and the critical temperatures, respectively [25]. However, the logarithmic bulk free energy function has singularities at φ = ±1. Therefore, for computational efficiency, a quartic polynomial approximation F(φ) = 1 4 (φ 2 -1) 2 has been used instead of the logarithmic potential. The AC equation has been successfully applied as a building block equation for modeling many scientifically and industrially important problems such as crystal growth, image segmentation, motion by mean curvature, tissue growth, volume repairing and smoothing, topology optimization, volume reconstruction from point cloud and slice data, and multiphase fluid flows [24]. However, the AC equation with the classical quartic polynomial function has a limitation in preserving structures. For example, as we will show that with a numerical experiment in the later section, if two separated components are close to each other, then they may merge with each other. To resolve this problem, we use the hAC equation which has a good structure preserving property. As shown in Fig. 1, the higher the order p, the larger the free energy barrier. This fact implies that once two different phases are separated, they remain separated, which is a good feature in modeling the dynamics of complex shapes.
Compared to the AC equation, the hAC equation is less studied numerically [24]. The main purpose of this article is to develop a second-order energy stable method for the hAC equation, which is based on the CS idea. For F p (φ), one can split into 1 4 φ 2p + 1 4 and Figure 1 Plot of the free energy F p (φ) = 1 4 (φ p -1) 2 with different orders: p = 2, 4, 10 -1 2 φ p , however, the resulting method is highly nonlinear thus its numerical implementation is complicated. In this study, we propose a linear CS that can be applied to all p, by placing 1 4 (φ p -1) 2 in the concave part. In addition, we combine the linear CS with a strong-stability-preserving (SSP) IMEX-RK method [26] to obtain temporally secondorder accurate and unconditionally energy stable scheme.
The contents of this article are as follows. We present the linear CS scheme with an auxiliary term in Sect. 2. In Sect. 3, we propose the second-order CS method, prove unconditional energy stability, and describe the computational implementation. In Sect. 4, we present computational tests to demonstrate the performance of the method. Conclusions are given in Sect. 5.

Linear convex splitting with an auxiliary term
It is well known that the AC equation preserves the maximum principle [27]. Thus, like a common practice to consider the AC and Cahn-Hilliard equations with the truncated Ginzburg-Landau double-well potential [5,[28][29][30][31][32], we can replace 1 4 Remark 1 With the replacement of 1 4 (φ p -1) 2 byF p (φ), we have Therefore, the hAC equation can be rewritten as follows: Next, we consider the following splitting where β ≥ 0 is a constant.
Then, we obtain For E e (φ), we get Then, we have

Lemma 2
The convexity condition of E c (φ) and E e (φ) results in the following inequality: where (·, ·) L 2 denotes the L 2 -inner product with respect to .

Unconditionally stable second-order scheme
Next, we propose a second-order CS scheme for the hAC equation by combining the linear CS (3) with the specially designed second-order (three-stage) SSP-IMEX-RK scheme: where a 10 , a 11 , a 20 , a 21 , a 22 , b 1 , b 2 satisfy the second-order conditions and the stability condition [26], and are as follows:

Theorem 1
The method (7) with β ≥ A is unconditionally energy stable, i.e., we have the following inequality: for any time step t > 0.
This completes the proof.

Numerical implementation
We can rewrite the first-step of (7) as where L = β/ 2 -. Then, we have where we have used the zero Neumann boundary condition for φ. In the same way, (1) and (2) .

Convergence test
We investigate the performance of the proposed scheme with an initial condition As we can observe in Fig. 2, the spatial convergence of the scheme under grid refinement is evident.
Next, to calculate the temporal convergence rate, computations are done by fixing x = 1/128 and changing t = t f /2 12 , t f /2 11 , . . . , t f /2 7 . We use the quadruply over-resolved computational result as the reference numerical solution. Figure 3

Energy stability test
To study the energy stability of the proposed method, let us consider the following initial condition on = [0, 1] × [0, 1]: φ(x, y, 0) = rand(x, y), where rand(x, y) is a random number between [-0.9, 0.9]. We use = 0.005, β = p 2 /2, and x = y = 1/128. Figure 4 displays the temporal evolution of the discrete energy with various t for p = 4, 6, 8, 10. And the temporal evolution of the discrete energy with t = 2 19 2 ≈ 13 for p = 4, 6, 8, 10 is displayed in Fig. 5. For all p and large time steps, all the discrete total energies are temporally decreasing. Figure 6 displays the temporal evolution of φ(x, y, t) with t = 2 4 2 for all p.

Motion by mean curvature
Suppose a radially symmetric initial condition is given as follows:

Effect of p on the interfacial dynamics in 2D
To investigate the effect of p in the hAC equation on the interfacial dynamics in 2D, we consider two circles on = [0, 2]×[0, 1] (see the first column of Fig. 8(a)). Their centers are (0.667, 0.5) and (1.333, 0.5), respectively, and radii are 0.3, and we set φ(x, y, 0) to 1 inside the circles and -1 otherwise. We choose = 0.02, β = 0, x = y = 1/128, and t = 10 -5 . Figure 8 displays the evolution of φ(x, y, t) for p = 2 and 10. The initial two circles merge into one when p = 2, whereas shrink without merging when p = 10.

Effect of p on the interfacial dynamics in 3D
To examine the effect of p in the hAC equation on the interfacial dynamics in 3D, we consider a three-dimensional spiral on = [0, 1] × [0, 1] × [0, 1] (see the first column of Fig. 9) [24]. The width and height of the spiral are 8 x and 58 x, respectively, and we set φ(x, y, z, 0) to 1 inside the spiral and -1 otherwise. We use x = y = z = 1 64 , = x, β = p 2 /2, and t = 10 2 . Figures 9 and 10 display the evolution of φ(x, y, z, t) and energy for p = 4, 6, 8, 10. As p increases, the high-order polynomial free energy term F p (φ) 2 dx influences the evolution of φ more strongly than the interfacial energy term 1 2 |∇φ| 2 dx. As a result, φ loses its spiral shape by merge as p decreases, whereas φ still shows a spiral shape maintaining its interface as p increases.

Conclusions
We developed a linear, second-order, and unconditionally energy stable method for the hAC equation. In order to handle F p (φ) dx, β 2 2 φ 2 dx is added, which yields the linear convex-concave decomposition of the total energy for β ≥ A. In addition, we combined the linear CS with the specially designed second-order SSP-IMEX-RK scheme. We demonstrated that the proposed method is efficient and unconditionally energy stable, and second-order temporally accurate. In addition, we confirmed that the hAC equation has different interfacial dynamics depending on the value of p.