{'fraction_non_alphanumeric': 0.05242768595041322, 'fraction_numerical': 0.02365702479338843, 'mean_word_length': 4.278353326063249, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Energy Stable Second Order Linear Schemes for the Allen-Cahn Phase-Field Equation

Phase-field model is a powerful mathematical tool to study the dynamics of interface and morphology changes in fluid mechanics and material sciences. However, numerically solving a phase field model for a real problem is a challenge task due to the non-convexity of the bulk energy and the small interface thickness parameter in the equation. In this paper, we propose two stabilized second order semi-implicit linear schemes for the Allen-Cahn phase-field equation based on backward differentiation formula and Crank-Nicolson method, respectively. In both schemes, the nonlinear bulk force is treated explicitly with two second-order stabilization terms, which make the schemes unconditional energy stable and numerically efficient. By using a known result of the spectrum estimate of the linearized Allen-Cahn operator and some regularity estimate of the exact solution, we obtain an optimal second order convergence in time with a prefactor depending on the inverse of the characteristic interface thickness only in some lower polynomial order. Both 2-dimensional and 3-dimensional numerical results are presented to verify the accuracy and efficiency of proposed schemes.


Introduction
In this paper, we consider numerical approximation for the Allen-Cahn equation with Neumann boundary condition Here Ω ∈ R d ,d = 2,3 is a bounded domain with a locally Lipschitz boundary, n is the outward normal, T is a given time, φ(x,t) is the phase-field variable. f (φ), the bulk force, is the derivative of a given energy function F (φ), which is usually non-convex with two or more than two local minima. One commonly used energy function for two-phase problem is the double-well potential F (φ) = 1 4 (φ 2 − 1) 2 . ε is the thickness of the interface between two phases. γ, called mobility, is related to the characteristic relaxation time of the system. The homogeneous Neumann boundary condition implies that no mass loss occurs across the boundary walls. The equation (1.1) is introduced by Allen and Cahn [1] to describe the process of phase separation in multi-component alloy systems. It can be regarded as the L 2 gradient flow with respect to the Ginzburg-Landau energy functional The corresponding energy dissipation is given as Another popular phase field model is the Cahn-Hilliard equation, which is the the H −1 gradient flow with respect to the Ginzburg-Landau energy functional. It was originally introduced by Cahn and Hilliard [4] to describe the phase separation and coarsening phenomena in non-uniform systems such as alloys, glasses and polymer mixtures. The Allen-Cahn equation and the Cahn-Hilliard equation are widely used in modeling many interface problems due to their good mathematical properties (cf. e.g. [7, 14-16, 44, 53] ). However, the small parameter ε and the non-convexity of energy function F make the numerical approximation of a phase field equation a challenging task, especially the design of time marching schemes. It is well-known that if a fully explicit or implicit time marching scheme is used, a tiny time step-size is required for the semi-discretized scheme to be stable or uniquely solvable since the nonlinear function F is neither convex nor concave. A very popular approach to obtain unconditional stable time marching schemes is the so called convex splitting method which appears to be introduced by Elliott and Stuart [17], and popularized by Eyre [18], in which, the convex part of F (φ) is treated implicitly and the concave part of F (φ) is treated explicitly. This method has been applied to various gradient flows (see e.g. [2,19,21,22,40]). Traditional convex splitting schemes are first order accurate. Recently, several extensions to second order schemes were proposed based on either the Crank-Nicolson scheme (see e.g. [2,6,8,12,29,36]), or second order backward differentiation formula (BDF2) [35,47]. In all convex splitting schemes, no matter first order or second order, one usually obtains an uniquely solvable nonlinear convex problem at each time step.
There are another types of second order unconditional stable schemes for the phase field equations. In [13], Du and Nicolaides proposed a secant-line method which is energy stable and second order accurate. It is used and extended in several other works, e.g. [2,3,9,20,25,26,54]. Similar to the convex splitting method, the secant-line method leads to nonlinear semi-discretized system, which need special efforts to solve. Recently, an augmented Lagrange multiplier(ALM) method was proposed in [27,28] to get second order linear energy stable schemes. The idea is generalized as invariant energy quadratization (IEQ) by Yang et al. and successfully applied to handle several very complicated nonlinear phase-field models (see e.g. [30,[49][50][51]). Based on similar methodology, a new variant called scalar auxiliary variable (SAV) method is developed by Shen et al. [38,39]. In the ALM and IEQ approach, nonlinear semi-discretized systems are avoided, but one has to solve variable-coefficient systems, while in the SAV scheme, one only need to solve some linear systems with constant coefficients. Different to other methods, the energy in ALM, IEQ and SAV approach is a modified one which also depends on the auxiliary variable.
Our methods are based on semi-implicit discretization and stabilization skill. To improve the numerical stability of solving phase-field equations, semi-implicit schemes were proposed by Chen and Shen [5] and Zhu et al. [55]. Although not unconditionally stable, semi-implicit schemes allow much larger time step-sizes than explicit schemes. To further improve the stability, Xu and Tang proposed stabilized semi-implicit methods for epitaxial growth model in [45]. The proposed schemes have extraordinary numerical stability even though the mathematical proof of the stability is not complete. Similar schemes was developed for phase field equation by He et al. [31] and Shen and Yang [40], where the latter one adopted a mixed form for the Cahn-Hilliard equation, by using a truncated double well potential such that the assumption f ′ (φ) ∞ ≤ L is satisfied, the unconditional energy stability was proved for the first order stabilized scheme. It is worth to mention that with no truncation made to f (φ), Li et al [33,34] proved that the energy stable property can be obtained as well, but a much larger stability constant need be used.
In this paper, we develop two second-order unconditionally energy stable linear schemes for the Allen-Cahn equation based on the schemes proposed in [45] and [40]. The energy dissipation is guaranteed by including two second order stabilization terms, the first one is directly from [45], the other one is inspired by the work [44]. We also carry out an optimal error estimate for the time semi-discretized schemes. For the phase field equations, the error bounds will depend on the factor of 1/ε exponentially if one uses a standard procedure. By using a spectrum estimate result of de Mottoni and Schatzman [10,11] and Chen [7] for the linearized Allen-Cahn operator, we are able to get an optimal error estimate with a prefactor depend on 1/ε only in some lower polynomial order for small ε. This spectrum estimate argument was first used by Feng and Prohl [23,24] for an implicit first order scheme for phase field equations. It was also applied by Kessler et al. [32] to derive a posteriori error estimate for adaptive time marching. Similar analysis for a first-order stabilized semi-implicit scheme of the Allen-Cahn equation in given by Yang [48]. Recently, Feng and Li [21] , Feng et al. [22] extended this spectrum estimate argument to first order convex splitting scheme coupled with interior penalty discontinuous Galerkin spatial discretization for Allen-Cahn and Cahn-Hilliard equation, respectively. To our best knowledge, our analysis is the first such result for second order linear schemes. In summary, the proposed methods have several merits: 1) They are second order accurate; 2) They lead to linear systems with constant coefficients after time discretization; 3) The stability and error analysis bases on weak formulations, so both finite element method and spectral method can be used for spatial discretization to satisfy discretized energy dissipation law. 4) The methods can be easily used in more complicated systems. Note that, similar approach can be extended to the Cahn-Hilliard equation [42,43], where Lipschitz condition of f is assumed based on physical intuition and the analyses are more tedious.
The remain parts of the paper is organized as follows. In Section 2, we present the two second-order stabilized schemes for the Allen-Cahn equation and prove they are energy stable. The error estimate to derive a convergence rate that does not depend on 1/ε exponentially is then constructed in Section 3. Detailed implementation and numerical experiments for problems in both 2-dimensional and 3-dimensional tensorproduct domain are presented in Section 4 to verify our theoretical results. We end the paper with some conclusions in Section 5.

The two second order stabilized linear schemes
We first introduce some notations which will be used throughout the paper. We use · m,p to denote the standard norm of the Sobolev space W m,p (Ω). In particular, we use · L p to denote the norm of W 0,p (Ω) = L p (Ω); · m to denote the norm of W m,2 (Ω) = H m (Ω); and · to denote the norm of W 0,2 (Ω) = L 2 (Ω). Let (·,·) represent the L 2 inner product. For p ≥ 0, we define H p 0 (Ω) := {u ∈ H p (Ω)|(u,1) = 0}, and denote L 2 0 (Ω) := H 0 0 (Ω). For any given function φ(t) of t, we use φ n to denote an approximation of φ(nτ ), where τ is the step-size.
where L is a non-negative constant. Remark 2.1. Note that the commonly used double well potential does not satisfy the above assumption. But, thanks to the maximum principle that the Allen-Cahn equation has (cf. e.g. [1,7,23,48]), the solution to equation (1.1) will be bounded by value −1 and 1 if the initial condition is bounded by −1 and 1. So it is safe to modify the doublewell energy F (φ) for |φ| larger than 1 to be quadratic growth without affecting the exact solution if the initial condition is bounded by −1 and 1, such that assumption (2.3) is satisfied. This argument also applies to the assumption (3.1) in next section.

15)
then the SL-BDF2 scheme is stable with any any non-negative A and B, including the case A = B = 0.
On the other hand side, if we take then the SL-BDF2 scheme is unconditional stable for any τ .
for the scheme (2.17), where we define Proof. Pairing the equation (2.17) with δ t φ n+1 , we get We use Taylor expansion atφ (2.25) (2.26) then the SL-CN scheme is unconditional stable for any τ .
On the other hand, by using the inequality Remark 2.4. To make SL-BDF2 and SL-CN scheme be unconditionally stable, i.e. stable for any time step size τ > 0, we need take A ∼ O(γ/ε 2 ). This seems that A need to be very large in a real simulation. But actually, it is not necessary. It was showed by Magaletti et al. [37] and Xu et al. [46] that the phase-field Cahn-Hilliard-Navier-Stokes model for binary fluids has a fast convergence with respect to ε when the phenomenological mobility γ ∼ O(ε 2 ). We expect a similar scaling of γ in the Allen-Cahn model, which means A is not necessary as large as O(1/ε 2 ). Furthermore, it is proved that the numerical interface for the Allen-Cahn equatoin converges with the rate O(ε 2 |lnε| 2 ) before the singularities appear [21,23], which suggests that we don't need to take ε as small as a physical interface. Remark 2.5.
Recently, Li, Qiao and Tang [34], Li and Qiao [33] studied several first order and second order stabilized semi-implicit Fourier schemes, respectively, for the Cahn-Hilliard equation with double-well potential Without a Lipschitz condition on F ′ (φ), they proved that those schemes are unconditionally stable when very large stability constant A used. For example, according to Theorem 1.3 in [33], for a classical second order semi-implicit stabilized scheme proposed by Xu and Tang [45] applied to the Cahn-Hilliard equation, the stabilization constant A need to be as large as O(|lnε| 2 /ε 8 ) to make the scheme unconditionally stable (Note that the A in [33] corresponds to εB in this paper). However, the constants A,B in our schemes are only of order O(γ/ε 2 ),O(1/ε), respectively. The reasons are in two aspects. Firstly, for the Allen-Cahn equation, since its solution satisfies a maximum principle, it is reasonable to modify the F defined in (2.29) for |φ| > 1, such that the Lipschitz condition (2.3) is satisfied. Secondly, we use two stabilization terms instead of only one stabilization term. Our approach can be extended to the Cahn-Hilliard equation with quadratic growth energy as well [42,43].

Convergence analysis
In this section, we shall establish the error estimate of the two proposed schemes for the Allen-Cahn equation in the norm of l ∞ (0,T ;L 2 ) ∩ l 2 (0,T ;H 1 ). We will shown that, if the interface is well developed in the initial condition, the error bounds depend on 1/ε only in some lower polynomial order for small ε. Let φ(t n ) be the exact solution at time t = t n to the Allen-Cahn equation (1.1) and φ n be the solution at time t = t n to the time discrete numerical scheme (2.4) (or (2.17)), we define error function e n := φ n − φ(t n ). Obviously e 0 = 0.
Before presenting the detailed error analysis, we first make some assumptions. For simplicity, we take γ = 1 in this section, and assume 0 < ε < 1. We use notation in the way that f g means that f ≤ Cg with positive constant C independent of τ,ε. Assumption 3.1. We make following assumptions on f : where L 2 is a non-negative constant.
Since the solution of Allen-Cahn equation satisfies maximum principle (see Remark 2.1), one can always modify f (φ) for large |φ| such that Assumption 3.1 hold without affecting the exact solution. Assumption 3.2.
(i) We assume that there exist non-negative constants σ 1 such that (ii) Assume that an appropriate scheme is used to calculate the numerical solution at first step, such that Then it is easy to get Given Assumption 3.1 3.2 (i), we have following estimates for the exact solution to the Allen-Cahn equation. Lemma 3.1. Let φ be the exact solution of (1.1), under the condition of Assumption 3.1 and 3.2 (i), the following regularities holds: (i) Pairing (3.11) with φ t and taking integration by parts on the second term, we get After integration over [0,∞] and using the inequality (3.2), we obtain (i). (ii) We differentiate (3.11) in time to obtain (3.14) Integrate (3.14) over [0,∞), yields ess sup The assertion then follows from (i) and the inequality (3.3). (iii) Testing (3.13) with φ tt , we get (3.16) and by using (i) and the inequality (3.4) of Assumption 3.2, we obtain (iii). (iv) We differentiate (3.13) in time to derive Testing (3.18) with −∆φ tt and using H 1 (Ω) ֒→ L 4 (Ω) for d ≤ 4, we have (3.20) The assertion then follows from (i) (ii) (iii) and the inequality (3.5). (v) Testing (3.18) with φ ttt , we have (3.21) Integrating in time yields (3.22) The assertion then follows from (i) (ii) (iii) and the inequality (3.5).
By taking τ ε 1 4 max{4σ1+7,σ0}+ 9 2(6−d) , we have So, by taking step-sizes as defined in (3.39), the finer error estimate for N + 1 step can be obtained, and the the proof is completed by mathematical induction.
We prove this by induction. Assuming that the finer estimate (3.77) holds for all first N time steps, the coarse estimate (3.55) leads to On the other hand, in Theorem 3.2 and 3.1, the step size need be smaller than ε 2 to guarantee the convergence, which is much stronger than the requirement for the schemes to be unconditionally stable. Remark 3.2. The proofs of 3.1 and Theorem 3.2 are inspired by the works [24], [32], [21] and [22] for first order convex splitting schemes. The main difference is that we use a mathematical induction to handle high order terms come from the L 3 term, while a generalized Gronwall lemma is used in [21], [22], and a continuation argument is used in [32].

Implementation and numerical results
In this section, we numerically verify our schemes are second order accurate in time and energy stable.
Given φ 0 , to start the second order schemes, we use following first order stabilized scheme with smaller time steps to generate φ 1 , We take ε = 0.075 and M = 63 and use random initial values φ 0 to test the stability and accuracy of the proposed schemes. For the 3-dimensional case, the initial value is given as {φ 0 (x i ,y j ,z k )} ∈ R 2M×2M×2M with x i ,y j ,z k are tensor product Legendre-Gauss quadrature points and φ 0 (x i ,y j ,z k ) is a uniformly distributed random number between −1 and 1 (shown in the first picture of Fig. 4.1);

Stability results
We present the required minimum values of A (resp. B) with different B (resp. A) and τ values for stably solving the Allen-Cahn equation (1.1) in 3-dimensional case in Table 4.1 (resp. 4.2). Here by "stably solving", we mean the energy keep dissipating in first 1024 time steps. The corresponding 2-dimensional results are given in Table 4.3 and 4.4. Those results are obtained by using initial value φ 0 , the results for the cases taking initial value φ 1 are similar. From those tables, we see that the maximum required A values are of order O(ε 2 ) and the maximum required B values are of order O(ε). For τ small enough, both schemes are stable with A = 0 and B = 0. This is consistent to our analysis result. On the other hand side, by using a nonzero B, e.g. B = 10, the requirement for a large A will be dramatically reduced.
To check the energy dissipation property, we present in Figure 4.2 the log-log plot of the energy versus time for two schemes using different time step-sizes. We see the energy decaying property is maintained.

Accuracy results
We take initial value φ 1 (see the second plot in Fig. 4.1) for (1.1) to test the accuracy of the two schemes in a 2-dimensional domain Ω = [−1,1] × [−1,1]. The Allen-Cahn equation with the time relaxation parameter γ = 0.5 are solved from t = 0 to T = 1.28. To calculate the numerical error, we use the numerical result generated using τ = 10 −4 as a reference of exact solution. The results are given in Table  4.5 and Table 4. 6. We see that the schemes are second order accuracy in both L 2 and H 1 norm.  We proposed two second order stabilized linear schemes, namely the SL-BDF2 and the SL-CN scheme, for the phase-field Allen-Cahn equation. In both schemes, the nonlinear bulk forces are treated explicitly with two additional linear stabilization terms to guarantee unconditionally energy stable. The schemes lead to linear systems with constant coefficients thus can be efficiently solved. An optimal error estimate is given by using a spectrum argument to remove the exponential dependence on 1/ε. The error analysis also holds for the special cases when one of the stabilization constants or both of them take zero values. Numerical results verified the stability and accuracy of the proposed schemes.     Table 4.5. The convergence of the SL-BDF2 scheme with B = 5, A = 10 for the Allen-Cahn equation with initial value φ 1 , parameter γ = 0.5,ε = 0.05. The errors are calculated at T = 1.28.