A Second-Order Crank–Nicolson Leap-Frog Scheme for the Modified Phase Field Crystal Model with Long-Range Interaction

In this paper, we construct a fully discrete and decoupled Crank–Nicolson Leap-Frog (CNLF) scheme for solving the modified phase field crystal model (MPFC) with long-range interaction. The idea of CNLF is to treat stiff terms implicity with Crank–Nicolson and to treat non-stiff terms explicitly with Leap-Frog. In addition, the scalar auxiliary variable (SAV) method is used to allow explicit treatment of the nonlinear potential, then, these technique combines with CNLF can lead to the highly efficient, fully decoupled and linear numerical scheme with constant coefficients at each time step. Furthermore, the Fourier spectral method is used for the spatial discretization. Finally, we show that the CNLF scheme is fully discrete, second-order decoupled and unconditionally stable. Ample numerical experiments in 2D and 3D are provided to demonstrate the accuracy, efficiency, and stability of the proposed method.


Introduction
It is well known that crystal phenomena (such as edge dislocation [1], deformation and plasticity in nanocrystalline materials [2], epitaxial growth, and zone refinement [3]) are usually described by phase-field crystal (PFC) models that preserve the original properties. The initial PFC model was presented by Elder and Grant [3,4], which mainly studied the growth dynamics of crystals at both spatial and time scales. In order to introduce a mechanism for simulating elastic interaction, the MPFC model was presented by Stefanovic et al. [5]. In addition, the MPFC model mainly describes the separate situation of time scales, which allows for the elastic relaxation of the crystal lattice on fast time scales, but other evolutionary processes are usually on slower time scales. The MPFC equation is given as follows: where Ω is a domain in R d (d = 1, 2, 3), f (φ) = F (φ) = φ 3 − φ, 0 < < 1, φ is the density field, µ is the variation of F known as the chemical potential, M > 0 is a mobility coefficient and β denotes a non-negative constant. The PFC model is obtained from the nonlocal free energy functional with respect to the gradient flows and the equilibrium state of system commonly is the minimum of the energy functional. Hence, the total free energy functional E tol (φ) of the above governing Equation (1) is given as follows: The energy functional E tol (φ) decreases over time, i.e., the MPFC model satisfies the dissipation of energy law. In this study, we mainly focus on the MPFC model with longrange interaction, which is used to model the micro-phase separation process in diblock copolymers. It is well-known that the Ohta-Kawasaki (OK) model [6] mainly simulates the macro-phase separation process in diblock copolymer melts, which are chain molecules that makes up of two different and chemically incompatible segments connected by a covalent chemical bond. Therefore, on the basis of the Ohta-Kawasaki model, one can add a nonlocal operator that means the long-range interaction between the chain molecules. That is, we add a wave operator for the MPFC equation with long-range interaction. Nevertheless, it is challenging to design an accurate, easy-to-implement, efficient, and unconditionally stable scheme for this meaningful problem.
In this paper, we assume the periodic boundary condition and focus on time and space discretization. Therefore, we use the Fourier spectral method [7,8] for the spatial discretization for the system. In recent decades, many scholars focus on the numerical solution for the MPFC model and have conducted extensive related works and proposed various time-marching methods. For example, the semi-implicit method [9][10][11][12][13][14], the backward differentiation formula (BDF), and the Crank-Nicolson (CN) [15][16][17] method. In this paper, we mainly consider the CNLF method to deal with the MPFC model with long-range interaction. N. Hurl proved the stability of the unstable mode of the CNLF method, that is to say that CNLF does indeed control the unstable mode [18]. In order to eliminate the CFL time step retrictions of the CNLF scheme, N. Jiang proposed a new stabilized CNLF approach to the geophysical flow in [19,20], and the convergence and unconditional energy stability are proved in these works. Layton et al. [21] proposed a time-marching CNLF method for the uncoupling systems and demonstrated the conditional stability. The CNLF scheme was used to successfully solve geophysical flows, fast slow waves, and uncoupling groundwater-surface water flows in [22,23]. In addition, the CNLF scheme has some advantages, such as being a large, easy to decouple system, permitting the use of a large time step in actual numerical computations and parallel computation in decoupled sub-domain equations and avoiding some difficulties caused by nonlinear terms. Furthermore, we add a wave operator to the MPFC equation with long-range interaction, which is more beneficial to solving the CNLF scheme.
In recent years, it has become extremely important to develop efficient and easy-toimplement algorithms for gradient flow models. There are some popular strategies, for example, the convex splitting (CS) [24,25] methods, the invariant energy quadratization (IEQ) [26,27] and scalar auxiliary variable(SAV) [28][29][30] methods, etc. One can design first-and second-order unconditionally energy stable schemes by using the CS method; however, the second-order scheme of the CS method is nonlinear and implicit, which can potentially cause difficulties in numerical computation and some normally iterative steps. The IEQ method requires solving some coupled linear equations with complicated variable coefficients at each time step, which leads to high computational costs. Compared with the above two methods, the SAV method is more efficient and easy to implement, as it just needs to solve a purely linear system with constant coefficients at each time step. In this paper, we derive an equivalent PDE system by replacing the original variable with some new auxiliary variables. After developing a second-order CNLF scheme for the reformulated PDE system, we only need three steps to deal with two fully decoupled, sixth-order, purely linear equations for the MPFC equation with long-range interaction.
The main goal of this work is to construct a CNLF approach for solving the MPFC model with long-range interaction and prove the properties of the CNLF scheme. Our main contributions are: • We construct a second-order decoupled CNLF scheme by combining it with the SAV method in time and the Fourier spectral method in space and carry out the unique solvability, mass conservation, and stability for the MPFC equation with long-range interaction.
• The presented scheme is very easy to implement and highly efficient. We only need to solve two six-order, fully decoupled and linear equations with some constant coefficients at each time step. • We compare the second-order accuracy and errors of the presented scheme with the CN and BDF2 schemes numerically. In addition, we display the dynamic evolution of the energy of phase variable in experiments with long-range interaction.
The rest of this paper can be split into four sections. In Section 2, we introduce the governing systems and state the corresponding energy dissipation law. In Section 3, we design a second-order fully discrete and decoupled CNLF scheme by applying the SAV approach in time and the Fourier spectral method in space for the MPFC model with long-range interaction. The mass conservation and unconditional stability of the proposed scheme are proven rigorously, and we demonstrate the uniform-in-time H 2 bound for the numerical solution. In Section 4, the effectiveness and unconditional stability of the proposed method are depicted through various numerical experiments. Finally, conclusions can be described in Section 5.

Governing Systems
In this section, we focus on the following MPFC model with long-range interaction [17]: which satisfies the initial condition φ| t=0 = φ 0 (x) and φ t | t=0 = 0. In addition, Let us defineφ as the mean value of φ over Ω, i.e.,φ = 1 |Ω| Ω φdx.β and β are two positive parameters, and σ denotes the strength of the long-range interaction. If σ = 0 andβ = 0, the above model is the classical PFC model. If σ = 0, σ(φ −φ) is a nonlocal potential term. The nonlocal term is derived from the following second term of the free energy functional E(φ): where the first term is the so-called Swift-Hohenberg-type [31] functional and the second term is the nonlocal potential with long-range interaction. We first give some notations here. Let (·, ·) and · denote the standard L 2 inner product and the corresponding induced L 2 norm. For each s ≥ 0, we use (·, ·) s and · s to denote the inner product and norm in H s (Ω). Note that H 0 (Ω) = L 2 (Ω). We define Hilbert spaces L 2 0 (Ω) = {u ∈ L 2 (Ω)|(u, 1) = 0}, L 2 per (Ω) = {u ∈ L 2 (Ω)|u is periodic on ∂Ω}, H s per = {u ∈ H s (Ω)|u is periodic on ∂Ω}, and H −s per (Ω) is the dual space of H s per (Ω). Then, let us introduce the inverse Laplace operator (−∆) −1 and the H −1 per inner product. Suppose f ∈ L 2 0 (Ω) and let u ∈ H 2 per (Ω) ∩ L 2 0 (Ω) be the solution of the periodic boundary value problem: − ∆u = f in Ω.
Then, we define u = (−∆) −1 f . Given f , g ∈ L 2 0 (Ω), the H −1 per inner product can be defined by: where u 1 , u 2 ∈ H 2 per (Ω) ∩ L 2 0 (Ω) are the solutions of the periodic boundary value problems −∆u 1 = f and −∆u 2 = g, respectively. The H −1 norm is defined by We can easily obtain: In Formula (4), G denotes the Green's function of −∆ with periodic boundary conditions, which has the property of ∆G(x − y) = −δ(x − y) [32]. Thus, we note that the second term of E(φ) can be rewritten as follows: and then (4) can be rewritten as: The MPFC model with long-range interaction (3) satisfies the following two properties.

Lemma 1.
Mass conservation: assuming that the initial condition satisfies φ t (x, 0) = 0, we have: Proof. By integrating the first equation of (3) over Ω with the periodic boundary condition for µ, we derive:β Then, the solution of the Equation (11) is: Therefore, as long as the initial condition satisfies φ t (x, 0) = 0, then Ω φ t (x, t)dx = 0, and we also have: Therefore, we complete the proof.

Numerical Scheme
In this section, we firstly convert system (3) to an equivalent form by employing the SAV approach, which aims to transform the nonlinear term of total energy into a simple quadratic form in terms of some new auxiliary variables, and we can also demonstrate the energy stability of the new equivalent system. To this end, two auxiliary variables are introduced as follows: where Therefore, the free energy functional is rewritten as: and the pseudo energy is obtained as: Now, a new but equivalent PDE system is derived as follows: where The variables are equipped with the compatible initial boundary conditions: It is worth noting that U depends only on time; thus, U does not require any boundary conditions, and other variables adopt periodic boundary conditions. Furthermore, the system (21)-(24) still keeps the energy dissipation law: Remark 1. It should be noted that system (21)-(24) is equivalent to the original system (3) with the same initial conditions. The energy (20) and its dissipation law (26) for the system (21)- (24) are exactly the same as the original energy (3) and its dissipation law (14), respectively. Then, in the following subsection, we will construct an efficient and fully discrete numerical scheme for the transformed system (21)-(24).

The Semi-Discrete CNLF Scheme
In this subsection, we design a second-order, time-marching scheme to deal with the MPFC problem with long-range interaction by employing the SAV method. In the following, we use δt > 0 to denote the time step and set t n = nδt for 0 ≤ n ≤ N with the final time T = Nδt. Assuming that (φ n−1 , ψ n−1 , U n−1 ) and (φ n , ψ n , U n ) are already calculated with n ≥ 1, we update (φ n+1 , ψ n+1 , U n+1 ): Theorem 1. The linear system resulting from the scheme (27)-(30) is uniquely solvable.
In addition, it is easy to derive that for any φ with the periodic boundary condition, we have Therefore, we have the unique solution: where γ n ≥ 0, since all the eigenvalues of the symmetric operator ℘ −1 (∆H(φ n )) are non-positive. Then, substituting (42) into (36) leads to: which implies the existence and uniqueness of the scheme (27)- (30) due to the fact that ℘ is a symmetric positive definite operator. Furthermore, it is very easy to solve the scheme (27)- (30) in the following procedure: (42) and U n+1 from (30) by solving another six-order decoupled linear equation ν 2 = ℘ −1 (q n−1 ); • Finally, based on the above two steps, we can obtain φ n+1 from (37) as: Therefore, we just need to compute ℘ −1 (∆H(φ n )) and ℘ −1 (q n−1 ); in other words, the total price of computing the CNLF scheme (27)-(30) at each time step is merely need to solve two completely decoupled, six-order, purely linear equations with the periodic boundary conditions. Remark 2. The above solution processes are similar to reference [17]. Since the linear operator ℘ is a constant coefficient matrix, we show that ℘ is symmetric and positive definite, and it also follows that the operator ℘ −1 (∆(·)) is non-positive. Thus, the existence and uniqueness of formula (36) are obtained. In addition, a detailed explanation of the operator ℘ −1 is provided in reference [17], and we leave it to the interested readers.

Remark 3.
In this paper, we rigorously demonstrate the unique solvability of the scheme (27)-(30) at the semi-discrete level. For the fully discrete scheme, we could be proven with similar analysis in Theorem 1. We leave the detailed proof of the unique solvability to the interested reader.

The Fully Discrete CNLF Scheme with Fourier Spectral Method in Space
We describe the Fourier spectral framework [7]. We denote the domain Ω = (0, L x ) × (0, L y ) uniformly with h x = L x /N x , h y = L y /N y , where N x and N y are positive integers. The Fourier approximation space is: where i = √ −1, ζ k = 2πk/L x , η l = 2πl/L y . Then, any function u(x, y) ∈ L 2 (Ω) can be approximated by: k,l e iζ k x e iη l y , where the Fourier coefficients are denoted: u k,l = u, e iζ k x e iη l y = 1 Ω |Ω| ue −i(ζ k x+η l y) dx.
In what follows, we take N = N x = N y for simplicity. The full discrete scheme with the Fourier spectral method in space can be constructed as follows: Assuming φ n h , ψ n h , U n and φ n−1 h , ψ n−1 h , U n−1 are known, and for ∀Θ, Λ and Υ ∈ S h , we update φ n+1 h , ψ n+1 h , U n+1 by solving with the initial conditions Remark 4. The MPFC model with long-range interaction is a sixth-order, nonlinear, damped wave equation, and the stiffness of the system exists in a higher-order linear operator. The stiff terms of the proposed scheme are treated implicitly, and the non-stiff terms are treated explicitly. Then, we can use the fast Fourier transform to directly invert the equation, so that we can obtain the numerical solution of the equation easily and quickly.

Remark 5.
In order to obtain the second-order accuracy of the fully discrete CNLF scheme, we need to compute all the initial values at t = t 1 , which can be derived by applying the first-order backward Euler and corrector methods to obtain (φ 1 h , ψ 1 h , U 1 ). For ∀Θ, Λ and Υ ∈ S h , we obtain: Then, the corrector method is used as follows: where φ h dx. In addition, (φ 0 h , ψ 0 h , U 0 ) is given by the initial conditions for the MPFC model with long-range interaction.

Mass Conservation and Unconditional Stability
In this subsection, we demonstrate the unconditional stability and mass conservation for the second-order fully discrete scheme (48)-(51).
We first show that the property of mass conservation still holds. Letting Θ = 1 in (48) and (56), we obtain: Letting Υ = 1 in (50) and (58), we obtain: Then, we establish the unconditional stability in the following theorem.
As a direct consequence of the unconditional stability, a uniform-in-time H 2 bound for the numerical solution is derived as follows.

Lemma 3. Let φ n
h be the solution of (48)-(51). C 0 and C 1 are independent of δt and h. Suppose that the initial data are sufficiently regular such that Then, we have φ n h H 2 ≤ C 1 for all n and h.
Proof. As a result of (62), we have: for any n ≥ 1. In addition, Using integration by parts and Cauchy inequality yields: Therefore, we obtain: This completes the proof of Lemma 3.

Numerical Experiments
In this section, we perform various numerical simulations to test the accuracy, stability, and efficiency of the CNLF scheme (48)-(51) for solving the MPFC problem with long-range interaction and show the energy stability. We mainly use the SAV approach in time and the Fourier spectral method in space; in addition, the fast Fourier transform is applied for all numerical simulations.
In order to verify the accuracy and efficiency of the proposed CNLF scheme, we show the L 2 errors at T = 1 of the phase variable φ between the numerical solution and the exact solution in Figures 1 and 2. Based on the same parameters, we compare the temporal convergence rates of the Crank-Nicolson scheme described by SAV-CN. The BDF2 scheme is described by SAV-BDF2 [17], and the second-order Crank-Nicolson and BDF2 schemes proposed based on the IEQ approach are denoted by IEQ-CN and IEQ-BDF2 [33], respectively. In Figure 1, it is easy to observe that the SAV-CNLF, SAV-CN, and SAV-BDF2 schemes are stable for all tested time steps and result in good approximations and corresponding orders of accuracy. This further proves that our scheme is second-order stable.
When using the IEQ method, we need to solve the linear system with variable coefficients at each time step. On the contrary, when using the SAV method, we only need to solve two linear systems with constant coefficients at each time step. Therefore, although IEQ-CN, IEQ-BDF2, and SAV-CNLF all have good approximation effects in Figure 2, the total CPU time including all various time step δt of the IEQ-CN and IEQ-BDF2 schemes is 2.88 s and 3.25 s, respectively, and the SAV-CNLF scheme only needs 0.94 s. That is, the IEQ-CN and IEQ-BDF2 schemes require more CPU time than SAV-CNLF. Hence, we can find that the SAV approach proposed in this paper is more effective than the IEQ approach. Combining the two figures, we can verify that our proposed scheme is accurate, effective, and easy to implement.
We set the appropriate parameters, which are = 0.025, M = 1, T = 100, and the computational domain is Ω = [0, 32] 2 . For spatial discretization, the 128 2 Fourier modes are used. In Figure 3, we display the evolutionary process of the energy E CNLF with different time steps δt = 0.01, 0.1, 1, 5, 10, 20, where the other parameters areβ = 1, β = 2, σ = 0.05. It means that the energy E CNLF is still decreasing over time at large steps. Meanwhile, we show that the mass is invariant over time. In Figure 4, we show how the evolutionary process of original energy E, pseudo energy ε, and modified discrete energy E CNLF changes over time with σ = 0, σ = 0.1, and σ = 0.5, respectively, where β = 0.01, β = 1, δt = 0.01. We discover that ε and E CNLF are decreasing in time, while the original energy E is becoming more perturbed as σ becomes bigger. That indicates the strength of the long-range interaction between the chain molecules is stronger and the vibrations are more intense. In addition, we show how the evolution of original energy E, pseudo energy ε, and modified discrete energy E CNLF changes over time with β = 0.05, β = 2, and β = 20, respectively, where σ = 0.1, δt = 0.01 in Figure 5. We observe that the behaviors of the MPFC model are consistent with the PFC model and ε is the same as E when β is large (β = 20, in the high damping case). In particular, when β is small (β = 0.05, lower damping case), E shows an oscillatory behavior, and ε differs from E. In Figure 6, we depict the changes in the three energies with differentβ using the same parameters σ = 0.1, β = 1. It can be observed that the energy has a great performance. Due to the second-order time derivative of φ giving rise to short wavelength oscillations, the original energy E has a little oscillation asβ increases. As shown in the above figures, the discrete energy E CNLF of the CNLF scheme is non-increasing in time.   Figure 4. The curves of energy E, pseudo energy ε, and modified discrete energy E CNLF with different parameters σ = 0, σ = 0.1, σ = 0.5 in sequence.  Figure 6. The curves of energy E, pseudo energy ε, and modified discrete energy E CNLF with different parametersβ = 0.1,β = 1,β = 10 in sequence.

Phase Transition Behaviors in 2D and 3D
In this subsection, we display the phase transition behaviors of the density field by using the second-order CNLF scheme in 2D and 3D, respectively.
Firstly, we mainly simulate the phase transition behaviors with time evolution. In the 2D case, with the numerical simulation domain of [0, 128] 2 , the mesh size is h = 1. For the proposed scheme (48)-(51), we select the initial data φ 0 = 0.1 + rand(), where rand() is a uniformly distributed random number between −0.01 and 0.01 at the grid points in 2D. In addition, the fixed parameters are set as = 0.25, T = 200, δt = 0.0005, M = 1, β = 1. In Figures 7 and 8, we display the evolutionary process of the phase transition with σ = 0.1, 0.2 andβ = 0.1, 0.2 at t = 20, 80, 150, 200. Moreover, one can see that under the appropriate parameters, the pattern will become hexagonal over time. The numerical results are consistent with the transition behavior in [26,34]. With the decrease in σ andβ, the phase variable reaches the steady-state faster and makes the steady-state pattern of the phase transition more uniform. Finally, we show the corresponding energy and mass associated with the phase variable φ in Figures 9 and 10, respectively. We demonstrate that the energystable and mass-conserving properties of the proposed scheme are consistent with the theoretical result.
For the 3D simulation, we firstly set the numerical simulation domain [0, 50] 3 and use 64 3 Fourier modes. For the CNLF scheme, we set the appropriate initial data φ 0 = −0.35 + rand(), where rand() is a uniformly distributed random number between −0.01 and 0.01 at the grid points in 3D. In addition, the parameters are chosen as = 0.56, T = 80, δt = 0.0005, M = 1, σ = 0.1, β = 1. In Figures 11 and 12, we show the steady-state microstructure of the phase transition behaviors withβ = 0.1, 0.2, and we can obtain similar results as the energy E CNLF and mass in Figure 13. It can be seen that the pattern becomes a more uniform lattice structure from stripes and the steady-state microstructure of the phase transition behaviors is achieved faster as the value ofβ increases. All of the numerical results are consistent with the transition behavior in [26,34]. In addition, the curves depict that the free energy is monotonically non-increasing over time.

Conclusions
In this paper, we have presented an efficient, fully discrete and unconditionally stable CNLF scheme for the MPFC equation with long-range interaction. Our proposed scheme conquer the inconvenience from nonlinearities by linearizing the nonlinear cubic term, and the advantage of the presented scheme is that, in each time step, only two decoupled six-order linear equations with constant coefficients need to be solved. Thus, it is easy to implement and efficient. In addition, we established the unique solvability at the semidiscrete level, and we rigorously proved the mass conservation and unconditional stability at the fully discrete level. Finally, numerical results for some benchmark simulations were presented to verify the efficiency and stability of the developed scheme.
Overall, the proposed CNLF scheme is accurate, efficient, and easy to implement, and the presented idea can be readily extended to study a broader class of multiphase hydrodynamic models for developing fully discrete, linear, and unconditionally stable schemes.
Author Contributions: Conceptualization, C.W., X.F. and L.Q.; methodology, C.W., X.F. and L.Q.; writing-original draft preparation, C.W.; writing-review and editing, C.W., X.F. and L.Q.; supervision X.F. and L.Q. All authors have read and agreed to the published version of the manuscript.