Second-Order Unconditionally Stable Direct Methods for Allen–Cahn and Conservative Allen–Cahn Equations on Surfaces

: This paper describes temporally second-order unconditionally stable direct methods for Allen–Cahn and conservative Allen–Cahn equations on surfaces. The discretization is performed via a surface mesh consisting of piecewise triangles and its dual-surface polygonal tessellation. We prove that the proposed schemes, which combine a linearly stabilized splitting scheme, are unconditionally energy-stable. The resulting system of discrete equations is linear and is simple to implement. Several numerical experiments are performed to demonstrate the performance of our proposed algorithm.


Introduction
The aim of this paper is to develop a second-order accurate method in time and space, for Allen-Cahn equation and conservative Allen-Cahn equation on surfaces: and Here the phase field φ(x, t) approximates the density of a binary mixture on the surface Ω. ∆ is the Laplacian operator. F(φ) = 0.25(φ 2 − 1) 2 is the Helmholtz free energy density. is the gradient energy parameter related to the interfacial energy. To simplify the presentation, here we consider the periodic boundary for the order parameter. It should be noted that our method can also be extended to cases involving Neumann and Dirichlet boundary conditions. The Allen-Cahn equation [1] has been applied to a wide range of problems such as motion by mean curvature [2,3], image analysis [4][5][6][7], grain growth [8,9], crystal growth [10,11], surface reconstruction [12,13], triply-periodic minimal surface [14] and two-phase fluid flows [15]. Various numerical methods have been developed to improve the numerical stability and accuracy of Allen-Cahn equation on the flat surface [16][17][18][19][20][21]. λ(t) in a conservative Allen-Cahn equation is the Lagrange multiplier, which can be defined as λ(t) = Ω F (φ)dx/( 2 Ω dx). It is well-known that Equation (2) satisfies the mass conservation: Equation (1) is derived from a constrained gradient flow in a L 2 space of the Helmholtz free energy functional: The total energy E (φ) decreases with time. That is, Here, ∂Ω is the boundary of Ω and n is the outward normal vector at the boundary. Here we have used the divergence theorem, since the solution φ(x, t) is assumed to be smooth. Equation (2) is also from a constrained gradient flow and makes E (φ) non-increasing with time: Here Equation (3) has been used. An important feature of the Allen-Cahn equation [1] and conservative Allen-Cahn equation is that they both satisfy the energy dissipation law. The conservative Allen-Cahn equation has also been studied and applied in the various fields on the flat surface [22][23][24][25][26][27][28][29][30]. Solving these two equations on non-flat surfaces are important issues, both in the geometrical and numerical sense. Recently, J.Kim et al. proposed a finite difference method for the conservative Allen-Cahn equation on various surfaces embedded in a narrow band domain in the three-dimensional space [31]. Based on the radial basis functions, V. Mohammadl et al. [32] presented a study of the conservative Allen-Cahn equation on surfaces. To the best of the authors' knowledge, there is no research work that has attempted to prove the second-order unconditionally stabilities of Allen-Cahn and conservative Allen-Cahn equations on surfaces. The goal of this paper is to develop second-order accurate schemes for the Allen-Cahn equation and conservative Allen-Cahn equation on surfaces. We will prove that the proposed schemes, which combine a linearly stabilized splitting scheme, are unconditionally energy-stable. The resulting system of discrete equations is linear and is simple to implement. Several numerical experiments will be performed to demonstrate the performance of our proposed algorithm. This paper is organized as follows. In Section 2, we derive the second-order schemes for Allen-Cahn equation and conservative Allen-Cahn equation on surfaces. We prove their properties. Numerical experiments are performed in Section 3. In Section 4, conclusions are given.

Numerical Method
Let V = {v i |1 ≤ i ≤ N v } be a triangular discretization of surface and F = {T k |1 ≤ k ≤ N F } be the set of triangles. For a vertex v ∈ Ω, let v j be the neighboring vertices of v for j = 0, 1, ..., p, where v 0 = v p . The vertices v j are labeled counterclockwise. Assume T j to be the triangle with vertices v, v j , and let v j+1 , and G j = (v j + v j+1 + v)/3. Let φ be defined on the set of vertices V. For a small regular surfaceΩ by using Green's formula, one has Ω ∆φdv = ∂Ω ∇φ, n d∂v.
Here, n is the outer normal vector of the regular surfaceΩ. < · > and · are the inner product and L 2 normal product operators, respectively. It is easy to obtain where D(v) = ∑ p−1 j=0 |T j | and |T j | is the area ofT j , which is a triangle with v, G j , and G j+1 . ∆ d is a discrete Laplace-Beltrami operator. On the other hand, we can get The normal vectors n v (G j ) and n v (G j+1 ) are defined as Here N j and N j+1 are the normal vector of triangles T j and T j+1 , respectively.
Equations (6)- (8) suggest that the approximating Laplace-Beltrami operator [33] can be defined as: Now Equation (1) can be discretized with second-order time accuracy as follows: whereφ n+ 1 2 = (3φ n − φ n−1 )/2. Equation (2) can be discretized with second-order time accuracy as follows: where Note that we can easily design a full Crank-Nicolson scheme that is a second-order two-level scheme. However, it requires solving a nonlinear equation at each time step. Compared to the fully Crank-Nicolson scheme, only an elliptic equation with constant coefficients is required to be solved at each time step, which makes it easy to implement. We define the discrete mass M(φ n ) as The discrete energy E (φ n ) and pseudo the discrete energyĒ d (φ n+1 , φ n ) are defined as: Here, ξ is the constant satisfying The discrete L 2 inner product at the vertex and centroid of the triangle face by Here, D j = D(v j ) for j = 1, 2, . . . , N V . In this work, we restrict our attention to the order parameter φ which is bounded. That is, there exists a constant such that |φ| ≤ L, because the Allen-Cahn equation satisfies the maximum principle. This condition can be satisfied by many physically relevant potentials by restricting the growth of F(φ) to a quadratic for |φ| ≥ L. Furthermore, if the maximum norm of initial condition φ 0 is bounded by 1, then the solution for the Allen-Cahn equation is also bounded based on our numerical tests. Now we will prove several lemmas.

Lemma 2.
Let φ and ψ be defined on a surface S, which is assumed to be closed.
Here we state Lemma 2 without proof, which can be obtained simply by summation-by-parts formulas and discrete first and second Green's identities. (15), we can show that

Proof. By Equation
Here The proof is complete.
Let G k G k+1 denote the edge of G k G k+1 . It is obvious that G k G k+1 belongs to the dual neighbor triangular mesh of the vertices s and v k+1 . Since the summation in Equation (21) is for the vertex v, we can rewrite it as Since all vertices are labeled counterclockwise, we get Substituting Equations (23)-(24) into Equation (22)

Proof. For Equation (11), we have
i.e., Here Equation (13) and Theorem 1 have been used. The proof is completed. Now, we prove the unconditional energy stability of these schemes under the condition α ≥ F (L).

Theorem 3.
If α ≥ F (L), three solutions (φ n+1 , φ n , φ n−1 ) of the scheme (11) satisfy: Proof. By multiplying Equation (11) by φ n+1 − φ n and summing by parts, we have By combining the above relations and using the definition of energy (15), we obtain Here, we have used the condition that α ≥ F (L). The proof is complete.
Theorem 5. Suppose that {φ n } are a sequence of solutions of the schemes (11) or (12) with the starting values Proof. By the last Theorem 3 and φ 0 = φ −1 ∈ [−1, 1], we have a chain of inequalities, Thus, the discrete version of the original energy is bounded and non-increasing in time, which follows from the above by the energy of the initial condition.
Therefore, our proposed method is unconditionally energy-stable with the suitable stabilizing parameter α ≥ F (L) = 3L 2 − 1, which is chosen to be α = 2 in practical simulation.

Spinodal Decomposition on Surfaces
Spinodal decomposition is a mechanism by which a solution of two or more components separates into different phases. The system separates into spatial regions rich in one species and poor in the other species and evolves into an equilibrium state with a lower overall free energy [34,35]. In this section, we will consider this problem for the Allen-Cahn equation and conservative Allen-Cahn equation with the initial condition φ(x, y, z, 0) = rand() on a sphere surface. Here rand() is a random number between −1 and 1. The calculation is run until t = 4 with time step ∆t = 0.01. We chose h = 0.025 and = 0.05. Figure 1a,b show the evolutions of spinodal decomposition of the Allen-Cahn equation and conservative Allen-Cahn equation, respectively. Observing these results, we can see that the positive phases are separated and gathered together. At the early stage, the evolutions with two methods are similar because the effect of the phase separation is dominant. Once the positive phases are fully separated, the dynamics of two equations become significantly different. On the later stage, the pattern of the phase deforms to a bar shape and shrinks with increasing speed, because the Allen-Cahn equation has the motion of mean curvature flow [1]. On the other hand, the bar shape remains due to the mass conservation of the conservative Allen-Cahn equation [31].

Convergence Test
To test the convergence rates of the scheme, we used a same smooth initial condition φ(x, y, z, 0) = sin(2πx) sin(2πy) sin(2πz) on the sphere surface. Here the time step ∆t = 0.01h is fixed. We set a series of grids as h = 0.2, 0.1, 0.05, and 0.025, which implies the spatial and temporal grids are refined both by a factor of 2. We expect the ratio of successive errors to increase by a factor of 2 because of O(∆t 2 + h 2 ) in the proposed scheme. Here we consider the Cauchy error, which is widely used to compute the error when we do not know the exact solution. The results are presented in Table 1, which confirm the second-order convergence rates in time and space.

Unconditional Stability
As mentioned above, our proposed schemes (11) and (12) are unconditionally stable and allow the use of large time steps. To demonstrate that, we performed a numerical experiment using a large time step, ∆t = 100. The initial condition and parameters were set as in Section 3.1. The calculations were run until time t = 1000. Observing the numerical solutions in Figure 3, we can see that the total energy evolutions are non-increasing, which suggests that our proposed schemes are indeed unconditionally stable.
(a) Allen-Cahn equation (b) Conservative Allen-Cahn equation Figure 3. The non-increasing discrete energy for our proposed method with a larger time step ∆t = 100. Note that we have normalized the total energy and mass at the initial time, respectively.

Mean Curvature Flow
Numerical simulations of surfaces evolving according to their mean curvature are presented. The conservative Allen-Cahn equation has the motion of mean curvature [31]. Following the reference [31], we consider two spherical caps on a unite sphere surface. Let R 1 (t) and R 2 (t) be the radii of two spherical caps at time t, respectively. We perform the initial condition as The calculation is run until t = 0.25 with a time step ∆t = 0.01. Here h = 0.025 and = 0.05 are chosen. Figure 4 shows the temporal evolution of two spherical caps on the unit sphere surface. Observing the results, we can see that the larger spherical cap grows and the smaller one shrinks.  The comparison R 1 and R 2 between exact solution and numerical solution are drawn in Figure 5. Note that the analytic solutions of R 1 and R 2 are obtained by the fourth-order Runge-Kutta method [31]. The result shows that the numerical solution is close to the exact solution.

Conclusions
In this work, we proposed second-order unconditionally stable direct methods for the Allen-Cahn and conservative Allen-Cahn equations on surfaces. The discretization is performed via a surface mesh consisting of piecewise triangles and its dual-surface polygonal tessellation. We have proved that the proposed schemes are unconditionally energy-stable. The proposed schemes are easy to implement because only linear elliptic equations should be solved at each time step. We performed several numerical experiments such as evolutions on simple and complex surfaces, convergence test, and stability test. The numerical experiments confirmed the efficiency of the proposed algorithm.
Author Contributions: All authors contributed equally to this work; B.X., Y.L. and Z.L. critically reviewed the manuscript. All authors have read and agreed to the published version of the manuscript.