{'fraction_non_alphanumeric': 0.05572926974074807, 'fraction_numerical': 0.024421135958836333, 'mean_word_length': 4.315800547022933, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical approximations of the Cahn-Hilliard and Allen-Cahn Equations with general nonlinear potential using the Invariant Energy Quadratization approach

In this paper, we carry out stability and error analyses for two first-order, semi-discrete time stepping schemes, which are based on the newly developed Invariant Energy Quadratization approach, for solving the well-known Cahn-Hilliard and Allen-Cahn equations with general nonlinear bulk potentials. Some reasonable sufficient conditions about boundedness and continuity of the nonlinear functional are given in order to obtain optimal error estimates. These conditions are naturally satisfied by two commonly used nonlinear potentials including the double-well potential and regularized logarithmic Flory-Huggins potential. The well-posedness, unconditional energy stabilities and optimal error estimates of the numerical schemes are proved rigorously.

1. Introduction. In this paper, we carry out stability analyses and error estimates for two first-order, semi-discrete time-stepping numerical schemes, that are based on the newly developed Invariant Energy Quadratization approach, for solving the well-known Cahn-Hilliard and Allen-Cahn equations with general nonlinear bulk potentials. These two equations, produced by the gradient flow (a diffusive process) in L 2 and H −1 spaces respectively, are typical representatives of the phase field (diffusive interface) method, which is a well-known efficient modeling and a numerical tool to resolve the motion of free interfaces between multiple material components. A main advantage of this approach is that it permits to solve the free interface problem by integrating a set of partial differential equations, thus avoiding the explicit treatment of the boundary conditions at the interface, see [3,4,10,19,20,23,25,27,32] and the references therein about the theoretical/numerical work, as well as their wide applications in various science and engineering fields.
For algorithms design of any phase field models related to the Allen-Cahn and Cahn-Hilliard equations, a significant goal is to verify the energy stable property at the discrete level irrespectively of the coarseness of the discretization in time and space. In what follows, those algorithms will be called unconditionally energy stable. Schemes with this property are especially preferred for solving phase field models since it is not only critical for the numerical scheme to capture the correct long-time dynamics of the system, but also provides sufficient flexibility for dealing with the stiffness issue, induced by the thin interface. Meanwhile, since the dynamics of coarse-graining (macroscopic) process may undergo rapid changes near the interface, the noncompliance of energy dissipation laws may lead to spurious numerical solutions if the mesh and time step size are not carefully controlled.
However, it is a challenging task to develop unconditionally energy stable schemes to resolve the stiffness issue from the thin interface. It is shown that the traditional fully-implicit or explicit discretization for the nonlinear term can cause very severe time step constraint (called conditionally energy stable) on the interfacial width [1,13,31]. Many efforts had been done (cf. [5, 12-14, 18, 24, 26, 28, 29, 31, 33, 34, 34, 38, 45, 47, 50] and the references therein) in order to remove this constraint. Among those, the convex splitting approach [12,29,33,35,46] and linear stabilization approach [5,24,31,34,38,45,50] are the two most commonly used numerical techniques. In the convex splitting approach, the convex part of the potential is treated implicitly and the concave part is treated explicitly. The scheme is unconditionally energy stable, however, it usually produces nonlinear type schemes, thus the implementation is complicated and the computational cost might be high. The linear stabilization approach treats the nonlinear term explicitly thus it is efficient and very easy to implement. In order to remove the time step constraint dependence on the interfacial width, a linear stabilizing term is added and its magnitude usually depends on the upper bound of the second order derivative of the nonlinear potential. Unfortunately, the upper bound is usually infinite, therefore the remedy is to modify the nonlinear potential, that in turn results in additional accuracy issues.
Moreover, it appears difficult to design second-order unconditionally energy stable schemes with a stabilized approach, although some progress has been made recently in [21,22].
Recently, a novel numerical method, called Invariant Energy Quadratization (IEQ) approach, has been developed and successfully applied to solve a variety of gradient flow models (cf. [6,8,15,17,36,37,[39][40][41][42][43][44]44,48,49]). Its essential idea is to transform the bulk potential into a quadratic form (since the nonlinear potential is usually bounded from below) using a set of new variables via the change of variables. Then, for the reformulated model in terms of the new variables that still retains the identical energy dissipation law, all nonlinear terms can be treated semi-explicitly, which then yields a linear system. This method bypasses those typical challenges such as the justification/adjustment of convexity or implicit/explicit terms, to arrive at first, second order or even higher order in time unconditionally energy stable schemes readily.
Although one might think that it would be natural to derive the corresponding error analysis for the IEQ type schemes by analogy with the proof of stability, the reality is quite the opposite. An exceptional case is for the double well potential where many terms can be simplified (cf. Remark 3.2). For general nonlinear potentials, an essential difficulty comes from the way of quadratization to introduce the new variable, that actually leads the new variable to act as a encapsulation, making it difficult to estimate the quantitative relation between the new and original variables. To the best of the authors' knowledge, we are not aware of any results about the error analysis of IEQ type schemes with general nonlinear potentials and almost all works had been focused on its remarkable unconditional energy stability. In view of the scarce of error analysis, the main objective of this paper is to derive optimal error estimates for the IEQ schemes for solving the Cahn-Hilliard and Allen-Cahn equations. We give some reasonable sufficient conditions about boundedness and continuity for the nonlinear functionals in order to obtain optimal error estimates. These conditions are naturally satisfied by the commonly used double well potential and regularized logarithmic Flory-Huggins potential. Moreover, the analytical approach developed in this paper is general enough and thus it can work as a standard framework to derive error estimates of IEQ type schemes for various gradient flow models with diverse nonlinear potentials.
The rest of paper is organized as follows. In Section 2, we give a brief introduction to the Cahn-Hilliard and Allen-Cahn systems. In Section 3, for the IEQ scheme of solving the fourth order Cahn-Hilliard equation, we study the unconditional energy stability, practical implementation, wellposedness, and derive the optimal error estimates. In Section 4, similar analytical work is performed for the second order Allen-Cahn equation. In Section 5, some concluding remarks are given.
2. PDE Models and their energy laws. We consider the following Lyapunov energy functional, is the nonlinear bulk potential, is an interface/penalty parameter causing stiffness issue into the PDE system when 1. There are two commonly used nonlinear bulk potentials for F (φ): (i) Ginzburg-Landau double-well type potential, cf. [3,23,28]: (ii) Logarithmic Flory-Huggins potential, cf. [2,3,9,11]: For either of these two nonlinear potentials, we note there always exists a positive constant A such that where we can simply choose A = 1 for both cases.
By applying the varational approach for the total free energy (2.1) in H −1 (Ω), we obtain the Cahn-Hilliard type system that reads as: where M is the mobility constant, w is the chemical potential, and f (φ) = F (φ). The initial condition is φ| (t=0) = φ 0 . For simplicity, we choose the suitable boundary conditions so that all complexities from the boundary integrals can be removed, i.e., By applying the varational approach for the total free energy (2.1) in L 2 (Ω), we obtain the Allen-Cahn type system that reads as: An important feature of the Cahn-Hilliard and Allen-Cahn equations is that they both satisfy energy dissipation law. For the Cahn-Hilliard system (2.5)-(2.6), by taking the L 2 inner product of (2.5) with −w, of (2.6) with φ t , performing integration by parts and summing up two equalities, we obtain For the Allen-Cahn system (2.8), by taking the L 2 inner product of with φ t , performing integration by parts, we obtain 3. Cahn-Hilliard equation. We first introduce some notations that will be used throughout the paper. We let L p (Ω) denote the usual Lebesgue space on Ω with the norm · L p . The inner product and norm in L 2 (Ω) are denoted by (·, ·) and · , respectively. W k,p (Ω) stands for the standard Sobolev spaces equipped with the standard Sobolev norms φ k,p . For p = 2, we write H k (Ω) for W k,2 (Ω), and the corresponding norm is φ k . We define two Sobolev spaces: 3.1. Unconditional energy stable linear scheme using the IEQ approach. We recall that the main challenge to develop efficient, unconditionally energy stable schemes for solving the system (2.5)-(2.6) lies in how to discretize the nonlinear term f (φ). We now consider the newly developed IEQ approach, see [6, 8, 15, 17, 36, 37, 39-44, 44, 48, 49]. Note F (φ) is bounded from below as (2.4), we choose a positive constant B such that B > A, and introduce a new variable U (φ) through the following quadratization formula, that is then the Cahn-Hilliard equation (2.5)-(2.6) can be rewritten as: with the initial conditions The new transformed system (3.3)-(3.5) also follows an energy dissipative law in terms of φ and the new variable U . By taking the L 2 inner product of (3.3) with −w, of (3.4) with φ t , of (3.5) with −2U , performing integration by parts and summing up all obtained equalities, we can obtain the energy dissipation law of the new system (3.3)-(3.5), that reads as Note the transformed system (3.3)-(3.5) is exactly equivalent to the original system (2.5)-(2.6) since (3.1) can be easily obtained by integrating (3.5) with respect to the time. Therefore, the energy law (3.7) for the transformed system is exactly the same as the energy law (2.10) for the original system for the time-continuous case. We emphasize that we will develop energy stable numerical schemes for time stepping of the new transformed system (3.3)-(3.5). Consequently, the proposed schemes should follow the new energy dissipation law (3.7) formally instead of the energy law for the original system (2.10).
Let δt > 0 denote the time step size and set t n = nδt for 0 ≤ n ≤ N with T = N δt, then a first-order, semi-discrete time discretization IEQ scheme for solving the new transformed system (3.3)-(3.5) reads as, (3.11) where H n = H(φ n ). The boundary conditions are as follows, We show the unconditionally energy stablilty of the scheme (3.9)-(3.11) as follows.
Theorem 3.1. The scheme (3.9)-(3.11) is unconditionally energy stable in the sense that Proof. By taking the L 2 inner product of (3.9) with −δtw n+1 and performing integration by parts, we derive By taking the L 2 inner product of (3.10) with φ n+1 − φ n , using the identity of and performing integration by parts, we get By taking the L 2 inner product of (3.11) with −2U n+1 and using (3.15), we get By combining the above equations together, we have which concludes the energy stability (3.13) by dropping some unnecessary positive terms.
Remark 3.1. The essential idea of the IEQ method is to transform the complicated nonlinear potentials into a simple quadratic form in terms of some new variables via a change of variables. When the nonlinear potential is the double well potential, one can choose B = 0 and thus this method is exactly the same as the so-called Lagrange multiplier method developed in [15]. We remark that the Lagrange multiplier method in [15] only works for the fourth order polynomial potential φ 4 since its derivative φ 3 can be decomposed into λ(φ)φ with λ(φ) = φ 2 which can be viewed as a Lagrange multiplier term. However, for other type potentials, the Lagrange multiplier method is not applicable. For example, one cannot separate a factor of φ from the logarithmic term.
On the contrary, the IEQ approach can handle various complex nonlinear terms as long as the corresponding nonlinear potentials are bounded from below, see the author's recent works in [6, 8, 15, 17, 36, 37, 39-44, 48, 49]). This simple way of quadratization provides some great advantages including (i) the complicated nonlinear potential is transformed into a quadratic polynomial that is much easier to handle; (ii) the derivative of the quadratic polynomial is linear, which makes it possible to develop linear scheme; and (iii) the quadratic formulation in terms of new variables can automatically maintain this property of positivity (or bounded from below) of nonlinear potentials.

Implementations and well-posedness.
The introduction of the new variable U may lead to an increase of the computational cost if one attempts to solve the coupled system (3.9)-(3.11). In fact, note the nonlinear coefficient H of the new variable U is treated explicitly in (3.10), thus we can rewrite it as follows, Using this equality, (3.10)-(3.11) can be rewritten as Therefore, in practice one can solve φ n+1 and w n+1 directly from (3.17)- (3.18). Once φ n+1 is obtained, U n+1 is automatically given by (3.16).
The scheme (3.18) includes a variable-coefficient 1 2 H n H n for the term φ n+1 that leads to timedependent dense matrices. Explicitly building those time-dependent dense matrices are extremely expensive (Note that, if one use finite element methods, the corresponding matrices will be sparse but time-dependent). So in practice, an efficient way is to use a conjugate gradient type solver with preconditioning (PCG), that only needs a subroutine to calculate the matrix-vector product instead of building the mass matrix explicitly. An efficient preconditioner is to replace the variable coefficient term 1 2 H n H n φ n+1 by a term with a constant coefficient being the maximum, i.e., 1 2 H n H n L ∞ φ. Furthermore, for any φ, ψ satisfy the boundary conditions (3.12), we have that means the linear operator P (φ) is symmetric (self-adjoint). Meanwhile, for φ that satisfies the boundary conditions (3.12) and Ω φdx = 0, we have where " = " is valid if and only if φ ≡ 0. These facts imply that the linear operator P is symmetric positive definite in H per (Ω) or H(Ω).
We now show the well-posedness of the scheme (3.17)-(3.18). In the following arguments, we will only consider the periodic boundary condition for convenience. For the case of homogenous Neumann boundary conditions, as long as the space of H per (Ω) is replaced by H(Ω), all theoretical results are still valid.
By taking the L 2 inner product of (3.17) with 1, we derive Ω Ω µ n+1 , and we define then the weak form for (3.17)-(3.18) is the system with unknowns φ, w ∈ (H per , H per )(Ω): We denote the above linear system as

Proof. (i) For any
where c 1 is a postive constant dependent on δt, M, and H n L ∞ . Therefore, the bilinear form (LX, Y ) is bounded.

Error estimates.
We now focus on the error estimates. To simplify the notations, without loss of generality, in the below, we let M = = 1. We use x y to denote there exists a constant C that is independent on δt and n such that x ≤ Cy.
To derive the error estimates of (3.9)-(3.11), we first give the following two Lemmas to establish the quantitative relation between the L 2 and H 1 norms of H(φ(t n )) − H(φ n ) and φ(t n ) − φ n under some reasonable assumptions.
for n ≤ M , where C 0 is a positive constant that is only dependent of C 0 , A and B.
Proof. First, for any θ ∈ [0, 1], (3.28) can ensure ψ n = θφ(t n ) + (1 − θ)φ n is uniformly bounded, i.e., ψ n ∈ [−2C 0 , 2C 0 ] for n ≤ M . Thus from assumption (ii), we can always find a positive constant Second, for any x, y ∈ (−A, ∞), by applying the intermediate value Theorem, there exists some value ξ ∈ (−A, ∞) that is between x and y, such that Thus, using (3.30), (3.31), and applying the intermediate value Theorem again, we derive where we have used Remark 3.2. Lemma 3.1 establishes a quantitative relation between the L 2 norm of H(φ(t n ))− H(φ n ) and φ(t n ) − φ n under some reasonable assumptions, where the Lipschitz property (3.31) of the quadratization function √ x + B (x > −A) plays a critical role. We note assumptions (i) and (ii) are automatically valid for the fourth order polynomial type double-well potential. Indeed, for the double well potential, one can choose B = 0 and H n = φ n , thus Lemma 3.1 is trivial and the error analysis is straightforward for this case, see [7,15,16]. However, for the logarithmic Flory-Huggins potential, B = 0 and assumption (ii) is not true since the domain is the open interval (0, 1) instead of (−∞, ∞). This issue can be overcome by extending the logarithmic functional near the domain boundary with a continuous, convex, piecewise function, see [9,11,36]. Such a regularized method is also a common practice to remove the difficulty about that any small fluctuation near the domain boundary (0, 1) can cause the overflow, numerically. In this way, the domain is regularized to (−∞, ∞) and thus the assumptions (ii) will become valid.
We now establish the error estimates for scheme (3.9)-(3.11). To this end, we formulate the Cahn-Hilliard system (3.3)-(3.5) as a truncation form: (3.47) We assume the exact solution φ, w, U of the system (3.3)-(3.5) possesses the following regularity conditions, Proof. Since the proof is rather straight forward, we leave this to the interested readers.
To derive the error estimates, we denote the error functions as e n φ = φ(t n ) − φ n , e n w = w(t n ) − w n , e n u = U (t n ) − U n , e n H = H(φ(t n )) − H(φ n ). We first prove the L ∞ stability of φ n , which plays the key point in the error estimates. From (3.48), we define a positive constant κ such that The preliminary result is given in the following lemma. Proof. We use the mathematical induction to prove this Lemma.
For n = 0, φ 0 L ∞ ≤ κ is true naturally. Assuming that φ n L ∞ ≤ κ is valid for all n ≤ M , we show φ M +1 L ∞ ≤ κ is also valid through the following two steps.
(i). By taking the L 2 inner product of (3.51) with δte n+1 w , we obtain   (3.62) where, H n L ∞ is bounded since φ n L ∞ is bounded, f is continuous, and F (x) > −A; (3.74) where we have used the Sobolev inequality φ L ∞ ≤ C Ω φ where C Ω is a constant that only depends on Ω. Thus Then we obtain the conclusion (3.55) by induction provided that δt ≤ s 0 = min(s 1 , Proof. Since φ n L ∞ ≤ κ for any 0 ≤ n ≤ T δt when δt ≤ s 0 , by following the first step in the proof of Lemma 3.4, we obtain the conclusion (3.76).

Allen-Cahn equation.
4.1. Unconditional energy stable linear scheme using the IEQ approach. For the Allen-Cahn equation, by using the same quadratization formula, we obtain a transformed PDE system as: where the new variable U is defined as (3.1). The initial conditions φ| t=0 = φ 0 , U | t=0 = F (φ 0 ) + B and boundary conditions are (2.9). By taking the L 2 inner product of (4.1) with φ t , and of (4.2) with −2U , performing integration by parts and summing up two equalities, we can obtain the energy dissipation law of the new system (4.1)-(4.2), that reads as The first-order, semi-discrete in time, IEQ scheme for solving the Allen-Cahn system (4.1)-(4.2) reads as follows, where H n = H(φ n ). The boundary conditions are: The unconditional energy stability of the scheme (4.4)-(4.5) is shown as follows.
Theorem 4.1. The scheme (4.4)-(4.5) is unconditionally energy stable in the sense that Proof. By taking the L 2 inner product of (4.4) with 1 M (φ n+1 − φ n ) and using (3.15), we get (4.8) By taking the L 2 inner product of (4.5) with 2U n+1 and using (3.15), we get By combining the above equations together, we have which concludes the energy stability (4.7) by dropping some unnecessary positive terms.

4.2.
Implementations and well-posedness. Similar to the Cahn-Hilliard equation, we can rewrite (4.5) as follows, Therefore, in practice one can solve φ n+1 directly from (4.11) and then update U n+1 by (4.10).
The weak form for (4.11) can be written as the following system with unknowns φ ∈ H 1 (Ω), M δt φ n − H n U n + 1 2 H n H n φ n . We denote the above linear system as The well-posedness of the linear system (4.13) is shown as follows. Proof. (i) For any φ, ψ ∈ H 1 (Ω), we have where c 1 is a positive constant dependent on δt, M, and H n L ∞ . Therefore, the bilinear form (Lφ, ψ) is bounded.

Error estimates.
For simplicity, we still assume = M = 1, and then formulate the Allen-Cahn system (4.1)-(4.2) as a truncation form: δt . (4.18) We assume the exact solution φ, U of the system (4.1)-(4.2) possesses the following regularity One can easily establish the following estimates for the truncation errors, provided that the exact solutions of (4.1)-(4.2) satisfy the regularity conditions (4.19).
Lemma 4.1. If the exact solutions of (4.1)-(4.2) satisfy the regularity conditions (4.19), then the truncation errors satisfy Proof. Since the proof is rather standard, due to the page limit, we leave it to the interested readers.
To derive the error estimates, we denote the error functions as By subtracting (4.4)-(4.5) from (4.16)-(4.17), we derive the error equations: Let κ = max 0≤t≤T φ(t) L ∞ + 1, we first prove the L ∞ stability of solution φ n . Proof. For n = 0, φ 0 L ∞ ≤ κ is true naturally. Assuming that φ n L ∞ ≤ κ is valid for all n ≤ M , we show φ M +1 L ∞ ≤ κ is also valid through the following two steps.
(i). By taking the L 2 inner product of (4.22) with e n+1 φ − e n φ , we get  (4.29) (4.32)   (ii). By using the H 2 regularity of elliptic problem (4.4) and the estimate (4.35), there exists a positive constant s 3 , such that the following inequality holds, that implies, there exists a constant s 4 such that, Therefore, by (4.37) and (4.35), we obtain   Proof. When δt ≤ s 0 , we have φ n L ∞ ≤ κ for any 0 ≤ n ≤ T δt . Thus, by following the first step in the proof of Lemma 4.2, we obtain the conclusion (4.39).

Concluding
Remarks. We carry out the stability and error analysis of two first order, semi-discrete time stepping schemes for solving the Cahn-Hilliard and Allen-Cahn equations. The idea to introduce some new variables through quadratization of the IEQ method is quite different from those traditional methods like implicit, explicit, nonlinear splitting, or other various tricky Taylor expansions. In the reformulated model in terms of the new variables, all nonlinear terms can be treated semi-explicitly, which then yields a well-posed linear system that still retains the energy dissipation law. For time continuous case, the energy law in terms of the new variables is identical to that in terms of old variables, while for time discrete case, it is the first order approximation. This method is widely applicable as long as the nonlinear potential is bounded from below. The induced scheme is linear, well-posed and unconditionally energy stable. For the Allen-Cahn equation, the linear system is further symmetric positive definite.
We recall that almost all numerical schemes for various gradient flow models based on the IEQ approach or its modified version SAV approach developed by Shen et. al. in [30], had been focused on their great unconditionally energy stable properties. The quadratic formula acts as an encapsulation and thus it is difficult to obtain the two quantitative relations derived in Lemma 3.1 and 3.2, that leads to the scarce of their corresponding error estimates. In this paper, some general, sufficient conditions about the boundedness/continuity of the nonlinear functional are given in order to obtain the optimal error estimates. These conditions are naturally satisfied by the commonly use polynomial type double well potential. For the logarithmic Flory-Huggins potentials of the regularized version, these conditions are appropriate as well. By utilizing the Lipschitz property of the quadratic formula together with the mathematical inductions, we rigorously derive the optimal error estimates for the first order IEQ schemes. Moreover, the analytical approach developed in this paper is general enough and thus it can work as a standard framework to derive error estimates of IEQ type schemes for various gradient flow models with diverse nonlinear potentials. Although we consider only time discrete schemes in this paper, the results here can be carried over to any consistent finite-dimensional Galerkin type approximations since the analyses are based on the variational formulation with all test functions in the same space as the space of the trial functions.