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

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 two-length-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 ($\geq 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\in\mathbb{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 second-order 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.   Problem formulation

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( \mathit{\boldsymbol{r}})$ be a real-valued or complex-valued function defined on $\mathbb{R}^d$ and let $\epsilon > 0$. We say that $\zeta\in \mathbb{R}^d$ is an $\epsilon$-almost period of $f$ if

A function $f$ is almost periodic on $\mathbb{R}^d$ if it is continuous and if for every $\epsilon$ there exists a number $L = L(\epsilon, f)$ such that any cube with the side length of $L$ on $\mathbb{R}^d$ contains an $\epsilon$-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 $\left\|{\cdot}\right\|^{2}_{AP} = \left<{\cdot, \cdot}\right>_{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.

Proposition 1. (1). An almost periodic function is uniformly continuous and bounded.

(2). If $f( \mathit{\boldsymbol{r}})$ and $g( \mathit{\boldsymbol{r}})$ are almost periodic functions, then $f( \mathit{\boldsymbol{r}}) + g( \mathit{\boldsymbol{r}})$ and $f( \mathit{\boldsymbol{r}}) \cdot g( \mathit{\boldsymbol{r}})$ are almost periodic functions.

These properties can be easily proven from the one-dimensional results [2].

Theorem 2.2. If $f( \mathit{\boldsymbol{r}})$ and $g( \mathit{\boldsymbol{r}})$ are almost periodic functions and differentiable, then the Green's identity holds in the almost periodic sense, i.e.,

Proof. Since $f( \mathit{\boldsymbol{r}})$ and $g( \mathit{\boldsymbol{r}})$ are almost periodic functions, then there exists a constant $M$ such that $\sup\limits_{ \mathit{\boldsymbol{r}}} \left\{ \left|f( \mathit{\boldsymbol{r}})\right|, \left|g( \mathit{\boldsymbol{r}})\right| \right\} \leq M$. Denote

where $\mathit{\boldsymbol{n}}$ is the outward normal of $\partial Q(R)$. In the $d$-dimensional space, we have

Therefore, we obtain the desired conclusion by

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

where $\psi( \mathit{\boldsymbol{r}})$, $\mathit{\boldsymbol{r}}\in\mathbb{R}^d$ $(d = 1, 2, 3)$, 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. $\varepsilon$ and $\alpha$ 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 defining $F = c^{2} \mathcal{F}$, $\psi( \mathit{\boldsymbol{r}}) = \sqrt{c} \phi( \mathit{\boldsymbol{r}})$, $\tilde{\varepsilon} = c \varepsilon$, and $\tilde{\alpha} = \sqrt{c} \alpha$. The rescaled functional is

where

To solve the iPFC free energy functional, we consider the following Allen-Cahn dynamic equation

The initial value is $\phi|_{t = 0} = \phi_{0}$. From Theorem 2.2, we can prove that the system (11) 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 $F_{1}(\phi) = \left<{ \mathcal{N}(\phi), 1}\right>_{AP} + C_{1} \geq 0$, where $C_{1}$ is a positive constant. We introduce a scalar auxiliary variable $\mathcal{R} = \sqrt{F_{1}(\phi)}$ to transform (1) into an equivalent system as

By taking the almost periodic inner products of (3a) with $\mathcal{W}$, and (3b) with $-\phi_{t}$, multiplying (3c) with $2 \mathcal{R}$ and adding them together, we have the following energy dissipation property

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}$ non-overlapping subintervals by the partition $0 = t^{0}<t^{1}<\cdots<t^{n}<\cdots<t^{N_{T}} = T$. The time step size is $\tau^{n} = t^{n+1} - t^{n}$. We denote $t^{n}+\tau^{n}/2$ as $t^{n+1/2}$. Let $\phi^{n} = \phi(t^{n})$ and $\mathcal{R}^{n} = \mathcal{R}(t^{n})$.

Scheme 3.1 (SAV/CN). For $n \geq 1$, given $\phi^{n}$, $\mathcal{R}^{n}$ and $\phi^{n-1}$, $\mathcal{R}^{n-1}$, we update $\phi^{n+1}$ and $\mathcal{R}^{n+1}$ by

where $\mathcal{R}^{n+1/2} = ( \mathcal{R}^{n+1}+ \mathcal{R}^{n})/2$, $\phi^{n+1/2} = (\phi^{n+1}+\phi^{n})/2$ and $\bar{\phi}^{n+1/2} = (3\phi^{n}-\phi^{n-1})/2$.

Theorem 3.1. The SAV/CN scheme satisfies the following energy dissipation mechanism

where

Proof. We take the almost periodic inner products of (4a) with $\tau^{n} \mathcal{W}^{n+1/2}$, and (4b) with $-(\phi^{n+1}-\phi^{n})$.

Multiplying (4c) with $2 \mathcal{R}^{n+1/2}$ yields

Adding (5), (6) and (7) together, we obtain

Since $\mathcal{R}^{n+1/2} = ( \mathcal{R}^{n+1}+ \mathcal{R}^{n})/2$ and $\phi^{n+1/2} = (\phi^{n+1}+\phi^{n})/2$, (8) can be recast as

The desired conclusion is obtained from the above equation.

Remark 3.1. It is noted that the modified energy $\mathcal{F}_{SAV/CN}^{n}$ is different from the original energy $\mathcal{F}(\phi^{n})$ since $\mathcal{R}^{n}$ is obtained from the iteration process.

Remark 3.2(The implementation of the SAV/CN scheme). Denote

Substituting (4b) and (4c) into (4a), we obtain

Eqn. (9) can be rewritten as

Taking the almost periodic inner product with $(I+\frac{1}{2}\tau^{n} \mathcal{G}^2)^{-1} u^{n+1/2}$ leads to

where

From (11), we get

Then we can directly calculate $\phi^{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} = \phi^{n} - \phi(t^{n})$, $w^{n+1} = \mathcal{W}^{n+1} - \mathcal{W}(t^{n+1})$, and $r^{n} = \mathcal{R}^{n} - \mathcal{R}(t^{n})$, we have

Theorem 3.2. For the Allen-Cahn dynamic equation, assume that $\phi^{0}$ is almost periodic and $\left\|{ \phi_{t} }\right\|^{2}_{AP}$ is bounded. Considering a uniform time partition, i.e., $\tau = \tau^{k}, k \leq N_{T}$, we have

where the constant $C$ is independent on $\tau$.

Proof. Let's subtract (3) from (4) at $t^{n+1/2}$

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

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 $| \mathcal{R}(t)| \leq C$ and

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

and

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, \cdots, k-1$, and using the Gronwall inequality.

3.3.   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 $\phi_{[0]}(t)$ has been calculated by some numerical schemes. Then we define the residual $R_{[0]}(t)$ and the error $\epsilon_{[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 $\tau^{n} = t^{n+1} - t^{n}$. To obtain the approximation solution $\phi_{[0]}(t)$, we rewrite the SAV/CN scheme as

where $\phi_{[0]}^{0} = \phi^{0}$, $\mathcal{R}_{[0]}^{n+1/2} = ( \mathcal{R}_{[0]}^{n+1} + \mathcal{R}_{[0]}^{n})/2$, $\phi_{[0]}^{n+1/2} = (\phi_{[0]}^{n+1} + \phi_{[0]}^{n})/2$ and $\bar{\phi}_{[0]}^{n+1/2} = (3\phi_{[0]}^{n} - \phi_{[0]}^{n-1})/2$. We adopt a similar strategy to discretize (24) as follows

where $\epsilon_{[0]}^{0} = \phi(0) - \phi_{[0]}^{0} = 0$, $\phi_{[0]}^{n+1/2, \epsilon} = \phi_{[0]}^{n+1/2} + \epsilon_{[0]}^{n+1/2}$, and $\bar{\phi}_{[0]}^{n+1/2, \epsilon} = \bar{\phi}_{[0]}^{n+1/2} + \bar{\epsilon}_{[0]}^{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\geq 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 $\phi( \mathit{\boldsymbol{r}})$ can be expanded as

where $\mathbf{B} \in \mathbb{R}^{n\times n}$ is an invertible matrix related to the $n$-dimensional primitive reciprocal lattice. The projection matrix $\mathcal{P} \in \mathbb{R}^{d\times 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 $\hat{\phi}( \mathit{\boldsymbol{h}})$ satisfies

In practice, let $\mathbf{N} = (N_1, N_2, \dots, N_n)\in \mathbb{N}^n$, and

The number of elements in the set is $N = (N_1+1)(N_2+1) \cdots (N_n+1)$. Using the PM, the SAV/CN scheme of full discretization reads

where

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, \cdots, 0)^{T} \in \mathbb{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} = \sqrt{2}$, $q_{2} = \sqrt{3}$, $\tilde{\varepsilon} = 10$ and $\tilde{\alpha} = 4$. The computational domain is $[0, 2\pi]$. Correspondingly, the projection matrix $\mathcal{P}$ and $\mathbf{B}$ in the PM both are $1$. The initial data is chosen as $\phi(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.

4.2.   Dynamic evolution

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(\pi/12)$, $\tilde{\varepsilon} = -2$ and $\tilde{\alpha} = 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

and the $\mathbf{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 fourth-order 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(\pi/12)$, $j = 1, \cdots, m$, $m = 3, 4, 5$. The other model parameters are $\tilde{\varepsilon} = -2$ and $\tilde{\alpha} = 2$. In the PM, the $\mathcal{P}$ and $\mathbf{B}$ are consistent with Subsec. 4.2. The spatial functions are also discretized in $4$-dimensional 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$-length-scale 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-length-scale 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 one-step 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.

Acknowledgments

This work is supported by National Natural Science Foundation of China (11771368), Hunan Science Foundation of China (2018JJ2376), Youth Project (18B057) and Key Project (19A500) of Education Department of Hunan Province of China, and Hunan Provincial Innovation Foundation for Postgraduate (0943/431000232).