High-order energy stable schemes of incommensurate phase-field crystal model

This article focuses on the development of high-order energy stable schemes for the multi-length-scale incommensurate phase-field crystal model which is able to study the phase behavior of aperiodic structures. These high-order schemes based on the scalar auxiliary variable (SAV) and spectral deferred correction (SDC) approaches are suitable for the L 2 gradient flow equation, i.e., the Allen-Cahn dynamic equation. Concretely, we propose a second-order Crank-Nicolson (CN) scheme of the SAV system, prove the energy dissipation law, and give the error estimate in the almost periodic function sense. Moreover, we use the SDC method to improve the computational accuracy of the SAV/CN scheme. Numerical results demonstrate the advantages of high-order numerical methods in numerical computations and show the influence of length-scales on the formation of ordered structures.


1.
Introduction. Aperiodic crystals, such as quasicrystals, are an important class of materials whose Fourier spectra cannot be all expressed by a set of basis vectors over the rational number field. The irrational coefficients give rise to the denseness of Fourier spectra which results in the difficulties in the theoretical study. Theoretically, a multiple characteristic length-scale model which possesses, at least, an irrational scale, has been widely applied to study the formation and thermodynamic stability of the aperiodic structures [1,6,13,10,14,7]. The early model could trace back to Bak's work on three-dimensional icosahedral quasicrystals. Since then, many related models have been proposed to study aperiodic structures, including for multicomponent systems [7]. Among these models, Lifshitz and Petrich (LP) modified the Swift-Hohenberg model and explicitly added an incommensurate twolength-scale potential into a Lyapunov functional to explore quasiperiodic patterns that emerged in Faraday experiments [13]. Recently, Savitz et al. extended the LP model from two-length-scale potential to multiple (≥ 3) length-scale potential and studied more kinds of quasicrystals [14]. The m-length-scale model actually is an incommensurate multi-length-scale phase-field-crystal (iPFC) model who owns a 4m (m ∈ N) order differential operator and nonlinear term in the energy functional. A high precision computation is helpful to study the phase behaviors of aperiodic structures. In this article, we will pay attention to the development of high-order numerical methods for the iPFC model.
Recently, various numerical methods have been proposed to solve phase-field equations including the convex splitting methods [17], the linear stabilized schemes [16], the invariant energy quadratization (IEQ) [18] and the scalar auxiliary variable (SAV) approaches [15,11]. The convex splitting method splits the energy functional into the convex and concave parts. The method treats the convex part implicitly and the concave one explicitly to keep the unconditional energy stability. While the application of this method is restricted by the form of the energy functional, such as double-well bulk energy. The linear stabilized scheme adds a penalty term to improve its stability and deals with the nonlinear terms explicitly for implementing it easily. However, such a stabilized approach makes it difficult to design secondorder unconditionally energy stable schemes. Assume that the nonlinear part has a lower bound, via introducing an auxiliary variable, the IEQ method transforms the energy functional into a quadratic form to keep the unconditional energy dissipation property. Similarly, the SAV approach introduces a scalar auxiliary variable by supposing the bounded bulk energy and obtains an unconditionally energy stable system. Besides these methods, the spectral deferred correction (SDC) [4,5] algorithm is an efficient strategy to improve the accuracy of the above schemes. In the paper, we will apply the SAV approach to solve the time-dependent equations and further use the SDC strategy to improve the numerical accuracy.
For aperiodic structures, two kinds of numerical methods, including the crystalline approximant method (CAM) and the projection method (PM) are usually used to discretize the quasiperiodic functions [9]. The CAM uses a big periodic structure to approximate an aperiodic structure and corresponds to the Diophantine approximation problem which studies how to approximate irrational numbers by rational numbers [3]. To evaluate aperiodic structures accurately, the CAM needs an extremely big computational region with an unacceptable computational burden to reduce the error of Diophantine approximation. To avoid the Diophantine approximation problem, the PM accurately describes aperiodic structures based on the fact that the aperiodic structure can be regarded as a periodic crystal in an appropriate higher-dimensional space. The PM uses one higher-dimensional periodic region to capture the essential characteristics and greatly reduces computational complexity.
The rest of the paper is organized as follows. In Section 2, we first outline some useful preliminaries of the almost periodic functions and then present the iPFC model. The energy dissipation principle of the L 2 gradient flow for the iPFC model in the almost periodic sense is also given. In Section 3, we propose a second-order energy stable scheme for the iPFC model and give the corresponding error estimate. And we use the SDC approach to improve its computational accuracy efficiently. Section 4 presents the convergence rates of these numerical schemes and discusses the advantages of high-order numerical approaches in simulating dynamic evolution. Moreover, we also show the influence of multiple length-scales on the thermodynamic stability of aperiodic structures. There are some conclusions in Section 5.
2.1. Preliminary. Aperiodic structures are space-filling phases without decay. A useful mathematical theory to describe aperiodic structures is the almost periodic function theory which is a generalization of continuous periodic functions. We define the notation of a d-dimensional almost periodic function. Definition 2.1. Let f (r) be a real-valued or complex-valued function defined on R d and let > 0. We say that ζ ∈ R d is an -almost period of f if A function f is almost periodic on R d if it is continuous and if for every there exists a number L = L( , f ) such that any cube with the side length of L on R d contains an -almost period of f .
The almost periodic L 2 inner product is defined by where Q(R) = [−R, R] d and |Q(R)| is the measure of Q(R). It is known that this inner product is well-defined [2]. We denote the corresponding norm by · 2 AP = ·, · AP . For simplicity, we denote the average spacial integral over the whole space as .
Some useful properties of d-dimensional almost periodic functions are presented as follows.
(1). An almost periodic function is uniformly continuous and bounded. (2). If f (r) and g(r) are almost periodic functions, then f (r) + g(r) and f (r) · g(r) are almost periodic functions.
These properties can be easily proven from the one-dimensional results [2].
2. If f (r) and g(r) are almost periodic functions and differentiable, then the Green's identity holds in the almost periodic sense, i.e., f (r), ∇g(r) AP = − ∇f (r), g(r) AP .
Proof. Since f (r) and g(r) are almost periodic functions, then there exists a constant M such that sup where n is the outward normal of ∂Q(R). In the d-dimensional space, we have Therefore, we obtain the desired conclusion by f (r), ∇g(r) AP = − f (r)∇g(r) dr = −− ∇f (r)g(r) dr = − ∇f (r), g(r) AP .

2.2.
Incommensurate phase-field crystal (iPFC) model. The simplest iPFC model may be the LP model which was originally proposed to study the bi-frequency excited Faraday wave [13]. Concretely, the free energy functional of the LP model can be written as where q is an irrational number depending on the property of quasiperiodic structures. The essential feature of the energy functional is the existence of two characteristic length-scales, 1 and q, which is a critical factor to stabilize quasi-periodic structures. The LP model has been also used to study the soft-matter quasicrystals [12]. However, this model could be able to globally stabilize the dodecagonal and decagonal rotational quasicrystals [8]. To generalize the iPFC model, Savitz et al. extended the interaction potential from two characteristic length-scales to multiple characteristic length-scales. In particular, the Lyapunov functional of an m-length-scale incommensurate system can be written as , is the order parameter corresponding to the density profile of the system. q j is the j-th characteristic length-scale which depends on the property of aperiodic structures. c is a positive model parameter to ensure that the principle wavelengths are near the critical wavelengths. ε and α are both model parameters related to physical conditions, such as temperature, pressure. To reduce the number of model parameters, the iPFC energy functional can be rescaled by To solve the iPFC free energy functional, we consider the following Allen-Cahn dynamic equation The initial value is φ| t=0 = φ 0 . From Theorem 2.2, we can prove that the system (1) satisfies the following energy dissipation law Then we impose the following mean zero constraint of order parameter on the iPFC model to ensure the mass conservation 3. Numerical methods. In this section, we propose the second-order unconditionally energy stable Crank-Nicolson scheme based on the SAV technique, and also give the error analysis. Further, we use the SDC strategy to improve the temporal accuracy of the second-order scheme.
3.1. The scalar auxiliary variable (SAV) approach. Let's suppose that We introduce a scalar auxiliary variable R = F 1 (φ) to transform (1) into an equivalent system as By taking the almost periodic inner products of (3a) with W, and (3b) with −φ t , multiplying (3c) with 2R and adding them together, we have the following energy dissipation property d dt The SAV approach can construct high-order unconditionally energy stable schemes. In this section, we discuss a second-order semi-discrete scheme based on the Crank-Nicolson method. Suppose the time interval [0, T ] is divided into N T nonoverlapping subintervals by the partition 0 = t 0 < t 1 < · · · < t n < · · · < t N T = T . The time step size is τ n = t n+1 − t n . We denote t n + τ n /2 as t n+1/2 . Let φ n = φ(t n ) and R n = R(t n ). Scheme 3.1 (SAV/CN). For n ≥ 1, given φ n , R n and φ n−1 , R n−1 , we update φ n+1 and R n+1 by where Theorem 3.1. The SAV/CN scheme satisfies the following energy dissipation mechanism Proof. We take the almost periodic inner products of (4a) with τ n W n+1/2 , and (4b) Multiplying (4c) with 2R n+1/2 yields Adding (5), (6) and (7) together, we obtain (8) can be recast as The desired conclusion is obtained from the above equation.
Remark 3.1. It is noted that the modified energy F n SAV /CN is different from the original energy F(φ n ) since R n is obtained from the iteration process.
Substituting (4b) and (4c) into (4a), we obtain Eqn. (9) can be rewritten as Taking the almost periodic inner product with (I + 1 2 τ n G 2 ) −1 u n+1/2 leads to where From (11), we get Then we can directly calculate φ n+1 via (10) and (12). 3.2. Error estimate. In this section, we will derive the error estimate of the SAV/CN scheme 3.1. Denote e n = φ n − φ(t n ), w n+1 = W n+1 − W(t n+1 ), and r n = R n − R(t n ), we have For the Allen-Cahn dynamic equation, assume that φ 0 is almost periodic and φ t 2 AP is bounded. Considering a uniform time partition, i.e., τ = τ k , k ≤ N T , we have where the constant C is independent on τ .
Proof. Let's subtract (3) from (4) where the truncation errors are given by With the Taylor expansion, the truncation errors can be rewritten as Firstly, making the almost periodic inner product of (13) with w n+1/2 yields e n+1 − e n , w n+1/2 Then its right-hand term can be bounded by Secondly, by taking the almost periodic inner products of (14) with −(e n+1 −e n ), we obtain Without the minus sign, the second term on the right-hand side in the above equation can be transformed into Note that |R(t)| ≤ C and ≤ C e n AP + e n−1 AP . The last term on the right-hand side of (18) can be estimated by Thirdly, multiplying the both sides of (15) by 2r n+1/2 , we obtain The first two terms on the right-hand side of (19) can be rewritten as where the last two terms on the right-hand side satisfy ≤ Cτ (r n+1 ) 2 + (r n ) 2 + e n 2 AP + e n−1 2 AP . Therefore the last term on the right-hand side of (19) can be bounded by With the above estimates, adding (16), (17) and (19) together leads to Then we obtain the desired conclusion by summing over n, n = 0, 1, · · · , k − 1, and using the Gronwall inequality.

Spectral deferred correction (SDC) approach.
The SDC approach [4,5] is an efficient method to improve the numerical precision for an existing scheme using spectral collocation points. Firstly, we introduce the basic idea of the SDC method. Integrating both sides of the Eqn. (1) with respect to t, we have Suppose the approximation solution φ [0] (t) has been calculated by some numerical schemes. Then we define the residual R [0] (t) and the error [0] (t) as

Replacing (22) into (20) yields
Then we insert (23) into (22) and subtract (21) By taking the derivative of both sides of the above equation, we obtain Then, we apply the SDC approach into the second-order SAV/CN scheme to improve the numerical accuracy. We denote this strategy as SAV/CN+SDC. To calculate the integral precisely, we adopt the following Chebyshev nodes, The time step size is τ n = t n+1 − t n . To obtain the approximation solution φ [0] (t), we rewrite the SAV/CN scheme as We adopt a similar strategy to discretize (24) as follows where , andφ n+1/2, . Substituting (25) and (26b) into (26a) and eliminating the residual by (21), we obtain the following linear solvable equation ) .
A more accurate solution can be updated by

3.4.
The projection method (PM) discretization. The PM is an accurate approach in computing aperiodic structures that can avoid the Diophantine approximation error. The PM is based on the fact that the d-dimensional aperiodic structure can be embedded into an n-dimensional periodic structure (n ≥ d). The dimensionality n is determined by the spectrum structures of the aperiodic structures. In particular, n is the number of linearly independent vectors over the rational number field which span the spectrum. In the PM, the d-dimensional order parameter φ(r) can be expanded as where B ∈ R n×n is an invertible matrix related to the n-dimensional primitive reciprocal lattice. The projection matrix P ∈ R d×n depends on the property of aperiodic structures. If we consider d-dimensional periodic structures, the projection matrix degenerates to a d-order identity matrix. Therefore, the PM provides a unified framework in calculating periodic and aperiodic crystals. More details about the PM can refer to [9]. The Fourier coefficientφ(h) satisfies In practice, let N = (N 1 , N 2 , . . . , N n ) ∈ N n , and The number of elements in the set is N = (N 1 + 1)(N 2 + 1) · · · (N n + 1). Using the PM, the SAV/CN scheme of full discretization readŝ In the above equations, the nonlinear terms are n-dimensional convolutions in the Fourier space. Directly computing them is extremely expensive. To avoid this, we use the pseudospectral method through the n-dimensional fast Fourier transformation to compute them efficiently in the n-dimensional time domain by simple multiplication. The mass conservation constraint (2) can be satisfied through where e 1 = (1, 0, · · · , 0) T ∈ R N . 4. Numerical results. In this section, we present several numerical examples to verify the accuracy of the SAV/CN and SAV/CN+SDC schemes and to illustrate the advantages of the higher-order scheme in dynamic evolution. We also show the influence of multiple length-scales on the thermodynamic stability of aperiodic structures.
4.1. Accuracy. In this subsection, we take the two characteristic length scale iPFC model in one-dimensional space to test the numerical accuracy of the SAV/CN and SAV/CN+SDC schemes. The model parameters are set as q 1 = √ 2, q 2 = √ 3,ε = 10 andα = 4. The computational domain is [0, 2π]. Correspondingly, the projection matrix P and B in the PM both are 1. The initial data is chosen as φ(x, 0) = sin(x). We check the temporal accuracy by taking the space discretization N = 128. The numerical solution of N T = 2048 is set as the reference solution. Tab. 1 shows the errors and convergence rates of the SAV/CN and SAV/CN+SDC schemes at T = 0.2. One can observe that the numerical accuracy of the SAV/CN scheme is second-order and can be improved to fourth-order by the SDC approach. In this subsection, we simulate the dynamic process of 2-dimensional dodecagonal quasicrystal (DDQC) using the iPFC model with two length-scales. The model parameters are q 1 = 1, q 2 = 2 cos(π/12),ε = −2 andα = 2. The initial value is the DDQC whose spectral distribution and real morphology are shown in Fig. 1. The big blue dot represents the origin and the others are the 24 basic Fourier modes located on the circles of radii q 1 and q 2 , respectively. When calculating the DDQC, we adopt the PM to discretize the  spatial functions in 4-dimensional space with 24 4 trigonometric functions. The 24 4 basis functions bring a negligible spatial error comparing the temporal error. The projection matrix in the PM is P = 1 cos(π/6) cos(π/3) 0 0 sin(π/6) sin(π/3) 1 , and the B is a 4-order identity matrix. We use the second-order SAV/CN scheme with N T = 256 to simulate the dynamic evolution for the DDQC. Fig. 2 shows the change tendency of the energy value. The corresponding morphologies at t = 0, 50, 100, 150, 200 are presented in Fig. 3. As one can see, the proposed scheme satisfies the energy dissipation law. It should be noted that the value C 1 in the SAV approach plays an important role in computational simulations. In our computation, the C 1 chosen as 10 16 which maintains the consistency of original and modified energy values. To show the role of the high-order methods in dynamic evolution, we give the reference energy E s which is calculated by the SAV/CN+SDC scheme in the case of N T = 2048. Using the reference value as the baseline, Fig. 4 shows the energy difference of the SAV/CN scheme with N T = 64, 128, 256 and SAV/CN+SDC method with N T = 32. With the increase of the temporal discretization N T , the energy difference decreases. However, the energy value obtained by the fourthorder SAV/CN+SDC scheme with N T = 32 is closer to the reference value than that of the SAV/CN scheme. It is demonstrated that the higher-order method shows more accurate results with less time discretization points in the dynamic simulation. We also give the morphologies of the crucial moment t = 0.4025 in Fig. 5. The morphology of the reference solution is also set as a reference value to clearly illustrate the differences. And the conclusions from these results are consistent with the energy evolution curves.

4.3.
The influence of multiple length-scales potential. In this subsection, we use the dodecagonal quasiperiodic phases as an example to investigate the influence of multi-length-scale potentials on the stability of aperiodic structures. Concretely, we consider three, four, and five multiple length scale iPFC models. The parameters in the potential function are set as q j = s j−1 , s = 2 cos(π/12), j = 1, · · · , m,   The spatial functions are also discretized in 4dimensional space with 24 4 basis functions. The SAV/CN approach with N T = 256 is adopted to simulate the dynamic process. Fig. 6 lists the corresponding energy evolution plots, initial values, and stationary states. The first row presents the energy evolutions of DDQCs with different length-scale potentials. The second and third rows give the spectral distributions and real morphologies of the m-lengthscale DDQCs, respectively. The last row shows the real morphologies of stationary solutions. From these results, one can see that the three-and four-length-scale potentials both result in 6-fold symmetric periodic crystal, while the five-lengthscale potential can obtain the 12-fold symmetric quasicrystals. Therefore increasing the number of characteristic length-scales in the iPFC model is helpful to stabilize the quasicrystals.

5.
Summary. For the iPFC model, we proposed a second-order SAV/CN scheme which is unconditionally energy stable in the almost periodic function sense and gave the error estimate. Meanwhile, we used the SDC approach to further improve the accuracy of the second-order scheme to the fourth-order method through a onestep correction. The PM was applied to discretize spatial functions for computing aperiodic structures to high accuracy. In numerical simulations, the efficiency of the numerical schemes was demonstrated via the numerical convergence rates and the comparison of the dynamic evolutions. By comparing the dynamic evolutions of the DDQCs with different length-scales, we found that increasing the number of characteristic length-scales in the iPFC model significantly is helpful to stabilize aperiodic structures. The scale parameters are set as q j = s j−1 , j = 1, · · · , m, s = 2 cos(π/12).