An Energy Stable BDF2 Fourier Pseudo-Spectral Numerical Scheme for the Square Phase Field Crystal Equation

TL;DR

This paper introduces an energy stable, second-order BDF2 Fourier pseudo-spectral scheme for the square phase field crystal equation, ensuring stability, accuracy, and efficient solvability for modeling atomic-scale crystal dynamics.

Contribution

It develops a novel energy stable numerical scheme combining BDF2 time discretization with Fourier pseudo-spectral spatial approximation for the complex SPFC equation.

Findings

Proves energy stability and unique solvability of the scheme.

Establishes optimal convergence rates in relevant norms.

Numerical experiments confirm robustness and high accuracy.

Abstract

In this paper we propose and analyze an energy stable numerical scheme for the square phase field crystal (SPFC) equation, a gradient flow modeling crystal dynamics at the atomic scale in space but on diffusive scales in time. In particular, a modification of the free energy potential to the standard phase field crystal model leads to a composition of the 4-Laplacian and the regular Laplacian operators. To overcome the difficulties associated with this highly nonlinear operator, we design numerical algorithms based on the structures of the individual energy terms. A Fourier pseudo-spectral approximation is taken in space, in such a way that the energy structure is respected, and summation-by-parts formulae enable us to study the discrete energy stability for such a high-order spatial discretization. In the temporal approximation, a second order BDF stencil is applied, combined with an appropriate extrapolation for the concave diffusion term(s). A second order artificial Douglas-Dupont-type regularization term is added to ensure energy stability, and a careful analysis leads to the artificial linear diffusion coming at an order lower that that of surface diffusion term. Such a choice leads to reduced numerical dissipation. At a theoretical level, the unique solvability, energy stability are established, and an optimal rate convergence analysis is derived in the $\ell^\infty (0,T; \ell^2) \cap \ell^2 (0,T; H_N^3)$ norm. In the numerical implementation, the preconditioned steepest descent (PSD) iteration is applied to solve for the composition of the highly nonlinear 4-Laplacian term and the standard Laplacian term, and a geometric convergence is assured for such an iteration. Finally, a few numerical experiments are presented, which confirm the robustness and accuracy of the proposed scheme.

Equations

$$E_{\mathrm{pfc}}(\varphi )={\int }_{\Omega }\left\{\frac{1}{4}{\varphi }^4+\frac{1-\epsilon }{2}{\varphi }^2-{\left|\nabla \varphi \right|}^2+\frac{1}{2}(\Delta \varphi {)}^2\right\}d\mathbf{x},0<\epsilon ,$$

$$E_{\mathrm{pfc}}(\varphi )={\int }_{\Omega }\left\{\frac{1}{4}{\varphi }^4-\frac{\epsilon }{2}{\varphi }^2+\frac{1}{2}\varphi {\left(1+\Delta \right)}^2\varphi \right\}d\mathbf{x},$$

$${\left(1+\Delta \right)}^2\varphi =\left(1+2\Delta +{\Delta }^2\right)\varphi =\varphi +2\Delta \varphi +{\Delta }^2\varphi .$$

$${\partial }_t\varphi =\Delta \mu ,\mu :={\delta }_{\varphi }{E}_{\mathrm{pfc}}={\varphi }^3-\epsilon \varphi +(1+\Delta {)}^2\varphi .$$

$$E_{\mathrm{spfc}}(\varphi )$$

$$={\int }_{\Omega }\left\{\frac{1}{4}{\left|\nabla \varphi \right|}^4-\frac{\epsilon }{2}{\varphi }^2+\frac{1}{2}\varphi {\left(1+\Delta \right)}^2\varphi \right\}d\mathbf{x}.$$

$${\partial }_t\varphi =\Delta \mu ,\mu :={\delta }_{\varphi }{E}_{\mathrm{spfc}}=-\nabla \cdot \left(|\nabla \varphi {|}^2\nabla \varphi \right)-\epsilon \varphi +{\left(1+\Delta \right)}^2\varphi .$$

$${\partial }_{\tau }\varphi ={\Delta }_z\left(-{\sigma }^4{\nabla }_z\cdot \left(|{\nabla }_z\varphi {|}^2{\nabla }_z\varphi \right)-\epsilon \varphi +{\left(1+{\sigma }^2{\Delta }_z\right)}^2\varphi \right),$$

$$A\ge \frac{{\epsilon }^2}{16}.$$

$$N=2K+1,\text{mbox}forsomeK\in \mathbb{N}.$$

$$\mathcal{G}_{N}:=\left\{f:\mathbb{Z}^{3}\to\mathbb{R}\ \middle|\ f\ \mbox{is $\Omega_{N}$-periodic}\right\}.$$

$$⟨f,g⟩:=h^3\sum _{i,j,k=0}^{N-1}{f}_{i,j,k}\cdot g_{i,j,k},{\Vert f\Vert }_2:=\sqrt{⟨f,f⟩}.$$

$$f_{i,j,k}=\sum _{\ell ,m,n=-K}^K{\hat{f}}_{\ell ,m,n}^N\exp \left(2\pi \mathrm{i}(\ell x_i+m{y}_j+n{z}_k)\right),$$

$${\hat{f}}_{\ell ,m,n}^N:=h^3\sum _{i,j,k=0}^{N-1}{f}_{i,j,k}\exp \left(-2\pi \mathrm{i}\left(\ell x_i+m{x}_j+n{x}_k\right)\right).$$

$${\mathcal{D}}_x{f}_{i,j,k}:=\sum _{\ell ,m,n=-K}^K\left(2\pi \mathrm{i}\ell \right){\hat{f}}_{\ell ,m,n}^N\exp \left(2\pi \mathrm{i}(\ell x_i+m{y}_j+n{z}_k)\right),$$

$${\mathcal{D}}_x^2{f}_{i,j,k}:=\sum _{\ell ,m,n=-K}^K\left(-4{\pi }^2{\ell }^2\right){\hat{f}}_{\ell ,m,n}^N\exp \left(2\pi \mathrm{i}(\ell x_i+m{y}_j+n{z}_k)\right).$$

$$\Delta_{N}f:=\left({\cal D}_{x}^{2}+{\cal D}_{y}^{2}+{\cal D}_{z}^{2}\right)f,\quad\nabla_{N}f:=\left(\begin{array}[]{c}{\cal D}_{x}f\\ {\cal D}_{y}f\\ {\cal D}_{z}f\end{array}\right),\quad\nabla_{N}\cdot\left(\begin{array}[]{c}f_{1}\\ f_{2}\\ f_{3}\end{array}\right):={\cal D}_{x}f_{1}+{\cal D}_{y}f_{2}+{\cal D}_{z}f_{3},$$

$${\nabla }_N\cdot {\nabla }_{N}f={\Delta }_{N}f.$$

$$f_S(x,y,z)=\sum _{\ell ,m,n=-K}^K{\hat{f}}_{\ell ,m,n}^N\exp \left(2\pi \mathrm{i}(\ell x+my+nz)\right).$$

$$Q_N(g{)}_{i,j,k}:=g(x_i,y_j,z_k),$$

$$Q_N(g{)}_{i,j,k}=\sum _{\ell ,m,n=-K}^K{\widehat{Q_N(g)}}_{\ell ,m,n}^N\exp \left(2\pi \mathrm{i}(\ell x_i+m{y}_j+n{z}_k)\right).$$

$$R_N(g)(x,y,z)=\sum _{\ell ,m,n=-K}^K{\widehat{Q_N(g)}}_{\ell ,m,n}^N\exp \left(2\pi \mathrm{i}(\ell x+my+nz)\right).$$

$$P_N(g)(x,y,z)=\sum _{\ell ,m,n=-K}^K{\hat{g}}_{\ell ,m,n}\exp \left(2\pi \mathrm{i}(\ell x+my+nz)\right),$$

$$\hat{g}_{\ell,m,n}=\int_{\Omega}g(x,y,z)\exp\left(-2\pi{\rm i}\left(\ell x+my+nz\right)\right)d\mbox{\boldmath$x$},$$

$${\hat{g}}_{\ell ,m,n}\ne {\widehat{Q_N(g)}}_{\ell ,m,n}^N.$$

$${\Vert R_N(\phi )\Vert }_{H^r}\le m^{\frac{d}{2}}{\Vert \phi \Vert }_{H^r},$$

$$\Vert \varphi -R_N(\varphi ){\Vert }_{H^k}\le C\Vert \varphi {\Vert }_{H^m}{h}^{m-k},\text{mbox}for0\le k\le m.$$

$${\mathcal{D}}_x{Q}_N(\varphi )=Q_N({\partial }_x{R}_N(\varphi )).$$

$$(-{\Delta }_N{)}^{\gamma }{f}_{i,j,k}:=\sum _{\ell ,m,n=-K}^K{\left(4{\pi }^2({\ell }^2+m^2+n^2)\right)}^{\gamma }{\hat{f}}_{\ell ,m,n}^N\exp \left(2\pi \mathrm{i}(\ell x_i+m{y}_j+n{z}_k)\right),$$

$$(-{\Delta }_N{)}^{-\gamma }{f}_{i,j,k}:=\sum _{\begin{array}{c}\ell ,m,n=-K\\ (\ell ,m,n)\ne 0\end{array}}^K{\left(4{\pi }^2({\ell }^2+m^2+n^2)\right)}^{-\gamma }{\hat{f}}_{\ell ,m,n}^N\exp \left(2\pi \mathrm{i}(\ell x_i+m{y}_j+n{z}_k)\right).$$

$$⟨f,{\Delta }_{N}g⟩=-⟨{\nabla }_{N}f,{\nabla }_{N}g⟩,⟨f,{\Delta }_N^{2}g⟩=⟨{\Delta }_{N}f,{\Delta }_{N}g⟩,⟨f,{\Delta }_N^{3}g⟩=-⟨{\nabla }_N{\Delta }_{N}f,{\nabla }_N{\Delta }_{N}g⟩.$$

Full text

Key words.   square phase field crystal equation, Fourier pseudo-spectral approximation, second order BDF stencil, energy stability, optimal rate convergence analysis, preconditioned steepest descent iteration

AMS Subject Classification   35K30, 35K55, 65K10, 65M12, 65M70

1 Introduction

The phase field crystal (PFC) model was proposed in [21] as a new approach to simulating crystal dynamics at the atomic scale in space but on diffusive scales in time. The model naturally incorporates elastic and plastic deformations, multiple crystal orientations and defects and has already been used to simulate a wide variety of microstructures, such as epitaxial thin film growth [22], grain growth [47], eutectic solidification [23], and dislocation formation and motion [47], to name a few. See the related review [41]. The key idea is that the phase variable describes a coarse-grained temporal average of the number density of atoms and the approach can be related to dynamic density functional theory [2, 39]. This method represents a significant advantage over other atomistic methods, such as molecular dynamics methods where the time steps are constrained by atomic-vibration time scales. In more detail, the dimensionless energy is given by the following form [21, 22, 48]

[TABLE]

where $\Omega\subset\mathbb{R}^{D}$, $D=2$ or 3, $\phi:\Omega\rightarrow\mathbb{R}$ is the atom density field. Typically the parameter $\varepsilon$ which represents a deviation from the melting temperature, satisfies $0<\varepsilon<1$, though it may be possible that $\varepsilon>1$. In this paper, $\phi$ is assumed to be periodic on the rectangular domain $\Omega$. Quite often, in the physics literature especially, the energy is rewritten as

[TABLE]

where

[TABLE]

The two views of the energy allow us to analyze the convexity structure from different perspectives. In (1.1), we view the quadratic term $-\left|\nabla\phi\right|^{2}$ as destabilizing (concave), and $\frac{1-\varepsilon}{2}\phi^{2}$ (and all other terms) as stabilizing (convex). This view is valid only when $a:=1-\varepsilon$ is positive, which, as we have said, may be slightly restrictive. By contrast, in (1.2), we view $-\frac{\varepsilon}{2}\phi^{2}$ as destabilizing and all other terms as stabilizing. In either case, the phase field crystal (PFC) equation is defined as

[TABLE]

When $\varepsilon<0$, the PFC equation can have solutions that exhibit spatial oscillations, and, typically in 2D, the peaks and valleys of $\phi$ are arranged in a hexagonal pattern. Such solutions are considered to represent “solid phase” solutions in the model. It is also possible to have “liquid phase” solutions, which are spatially uniform and constant; and these can even be in coexistence with the solid phase solutions to describe a crystal in equilibrium with its melt. See, for example, the book by Provatas and Elder [42].

Alternate lattice structures, such as “square” symmetry crystal lattices, are possible in 2D solutions. As mentioned in [22], and motivated by the work of [27], a different choice of nonlinear term in the PFC model is needed to obtain a square symmetry crystal lattice rather than the usual hexagonal structure. Specifically, such symmetries can be obtained [27] by replacing $\phi^{4}$ in (1.1) with $\left|\nabla\phi\right|^{4}$. (See also [53] for a related method.) Doing so one obtains an energy of the form

[TABLE]

We observe that there are essential similarities between this energy and the Aviles-Giga-type energy [1]. The square phase field crystal (SPFC) equation is given by the following dynamics

[TABLE]

We will assume for simplicity that $a=1-\varepsilon>0$. In this case, the energy will be bounded from below, a fact that we will prove later. For the standard PFC model and its modified version, there have been extensive numerical works [3, 4, 17, 36, 50, 51, 52, 57], et cetera, in the existing literature. Of course, because of its generality, the new SAV approach of Shen et al. [45] could be applied to the PFC and SPFC problems. For the SPFC equation, few if any simulation results exist, to our knowledge, though a closely related equation is solved in [27].

Regarding scaling and parameters, we first note that $\epsilon$ – which represents the deviation of the temperature from the freezing temperature – in the present setting is not necessarily a small parameter. However, in the standard PFC literature, the characteristic feature size, that is the unit cell size of the pattern, is generally taken to be $O(1)$ and the domain size, $L$, is typically large, perhaps several hundred unit cell sizes wide. If one rescales the domain so that $L=O(1)$, then the feature size is a small fraction of 1, and a small parameter appears in the equation:

[TABLE]

where $\tau$ and $z$ are rescaled time and space coordinates and $\sigma>0$ is a potentially small parameter. In particular, we see that the highest order diffusion term $\sigma^{4}\Delta_{z}^{2}\phi$ and the nonlinear term $-\sigma^{4}\nabla_{z}\cdot\left(|\nabla_{z}\phi|^{2}\nabla_{z}\phi\right)$ will always be of the same magnitude, so that the stiffness in the nonlinear term could never be removed via scaling. Implicit treatment of the highest order surface diffusion as well as the nonlinear term is necessary to ensure an energy stability at a theoretical level.

Regarding the nonlinearity, although the only difference between the standard PFC and SPFC equations is the replacement of $\phi^{4}$ by $|\nabla\phi|^{4}$ in the free energy functional, the analysis and numerical approximation of the later are much more challenging, especially when using pseudo-spectral approximations of spatial derivatives. For example, a nonlinear multigrid numerical solver has been very successfully applied to the standard PFC model in the framework of finite differences [3, 36, 52]. However, such a solver leads to a fairly poor numerical performance in the computation of the 4-Laplacian problem, due to its highly nonlinear nature. In addition, the $H^{-1}$ gradient flow pattern for the 4-Laplacian term makes the overall system even more difficult to solve. All these well-known difficulties have made a rigorous numerical analysis for the SPFC model an open problem until now.

In this article we consider two energy stable schemes for the SPFC equation, both with second order temporal accuracy and Fourier pseudo-spectral spatial discretization. They will differ in construction based only on whether one views $-\left|\nabla\phi\right|^{2}$ as destabilizing or one views $-\frac{\varepsilon}{2}\phi^{2}$ as destabilizing. Since the schemes are closely related, we will conduct the detailed convergence analysis only for one and just comment here and there on the differences one would encounter in the analysis of the other.

The SPFC equation (1.4) is structurally very different from the PFC model, with a much higher nonlinearity. However, its energy (1.3) can still be decomposed into purely convex and concave parts, so that a first order convex splitting scheme could be appropriately derived. On the other hand, first order temporal accuracy – especially in the convex splitting setting which is highly dissipative – is not satisfactory in practical computations, in particular due to the large time scale involved for the SPFC model. Therefore, a second-order-accurate in time, energy stable numerical algorithm is highly desired. Instead of the modified Crank-Nicolson approach for the gradient structure, which has been successfully applied to the Cahn-Hilliard [8, 11, 14, 15, 16, 18, 31, 32, 33, 34] and epitaxial thin film equation [44], we make use of a modified backward differentiation formula (BDF) approach. In more details, we apply the second order BDF concept to derive second order temporal accuracy, but modified so that the concave diffusion term is treated by an explicit extrapolation. Such an explicit treatment for the concave part of the chemical potential ensures the unique solvability of the scheme without sacrificing energy stability. An additional term $A{\Delta t}\Delta(\phi^{k+1}-\phi^{k})$ is added, which represents a second order Douglas-Dupont-type regularization. Moreover, a careful analysis shows that energy stability is guaranteed for the proposed numerical schemes, provided a mild condition is enforced:

[TABLE]

See also a related discussion in [12] for the epitaxial thin film growth model.

In the spatial discretization, we use Fourier pseudo-spectral approximation for its ability to capture more detailed structure with a reduced computational cost. Summation-by-parts formulae and aliasing error control techniques enable us to derive unique solvability and energy stability for the fully discrete numerical scheme. As a result of this discrete energy stability, a uniform-in-time discrete $H^{2}$ bound for the numerical solution becomes available. With such an $H_{N}^{2}$ bound at hand, we are able to derive a discrete $W_{N}^{1,6}$ bound for the numerical solution, uniform-in-time, with the help of discrete Sobolev embedding in the Fourier pseudo-spectral space. Such an embedding analysis cannot be derived from a straightforward calculation. Instead, detailed discrete Fourier analyses, combined with certain non-trivial aliasing error estimates, are involved in the derivation. In turn, these preliminary estimates enable one to obtain an optimal rate ($O({\Delta t}^{2}+h^{m})$) convergence analysis for the proposed numerical scheme, in the $\ell^{\infty}(0,T;\ell^{2})\cap\ell^{2}(0,T;H_{N}^{3})$ norm.

There have been related works of spectral/pseudo-spectral differentiation to various gradient flow models; see the related references [14, 37], et cetera. On the other hand, such a differentiation turns out to be a global operator in space, and this feature leads to great challenges in the numerical implementations, especially in the case of an implicit treatment of nonlinear terms. Moreover, due to the composition of the highly nonlinear 4-Laplacian and the regular Laplacian operators, the numerical implementation becomes even more challenging. As a result, an efficient solver for a regularized p-Laplacian equation, in the discrete $H^{-1}$ space, has to be utilized. In a recent work [25], a preconditioned steepest descent (PSD) algorithm was proposed for such problems. At each iteration stage, only a purely linear elliptic equation needs to be solved to obtain a search direction, and the numerical efficiency for such an elliptic equation could be greatly improved with the use of FFT-based solvers. In turn, an optimization in the given search direction becomes one-dimensional, with its well-posedness assured by convexity arguments. Moreover, a geometric convergence of such an iteration could be theoretically derived, for both the $L^{2}$ and $H^{-1}$ gradient flow structures, so that a great improvement of the numerical efficiency is justified. For the $L^{2}$ gradient flow, the numerical comparison has been reported in [26] for the epitaxial thin film model, with an application of the Polak-Ribiére variant of NCG (nonlinear conjugate gradient) method [40], reported in [44, 49]. For the 4-Laplacian problem in an $H^{-1}$ gradient flow, the PSD solver in the finite difference version has been reported in [24] for the functionalized Cahn-Hilliard/Willmore model, while an efficient spectral/pseudo-spectral PSD solver for the composition of 4-Laplacian and regular Laplacian operators has not been available in the existing literature.

The outline of the paper is given as follows. In Section 2 we present the numerical scheme. First we review the Fourier pseudo-spectral approximation in space and recall an aliasing error control technique. Then we formulate the proposed numerical scheme. Subsequently, the unique solvability and energy stability analyses are provided in Section 3, and an optimal rate convergence analysis is established in Section 4. In addition, the details of the preconditioned steepest descent (PSD) solver are outlined in Section 5. Some numerical results are presented in Section 6. Finally, some concluding remarks are made in Section 7.

2 The numerical scheme

2.1 Fourier pseudo-spectral approximations

The Fourier pseudo-spectral method is also referred as the Fourier collocation spectral method. It is closely related to the Fourier spectral method, but complements the basis by an additional pseudo-spectral basis, which allows to represent functions on a quadrature grid. This simplifies the evaluation of certain operators, and can considerably speed up the calculation when using fast algorithms such as the fast Fourier transform (FFT); see the related descriptions in [5, 10, 13, 14, 29, 30, 35, 55, 56].

To simplify the notation in our pseudo-spectral analysis, we assume that the domain is given by $\Omega=(0,1)^{3}$, $N_{x}=N_{y}=N_{z}=:N\in\mathbb{N}$ and $N\cdot h=1$. We further assume that $N$ is odd:

[TABLE]

The analyses for more general cases are a bit more tedious, but can be carried out without essential difficulty. The spatial variables are evaluated on the standard 3D numerical grid $\Omega_{N}$, which is defined by grid points $(x_{i},y_{j},z_{k})$, with $x_{i}=ih$, $y_{j}=jh$, $z_{k}=kh$, $0\leq i,j,k\leq 2K+1$. This description for three-dimensional mesh ($d=3$) can here and elsewhere be trivially modified for the two-dimensional case ($d=2$).

We define the grid function space

[TABLE]

Given any periodic grid functions $f,g\in\mathcal{G}_{N}$, the $\ell^{2}$ inner product and norm are defined as

[TABLE]

The zero-mean grid function subspace is denoted $\mathring{\mathcal{G}}_{N}:=\left\{f\in\mathcal{G}_{N}\ \middle|\ \langle f,1\rangle=:\overline{f}=0\right\}$. For $f\in\mathcal{G}_{N}$, we have the discrete Fourier expansion

[TABLE]

where the discrete Fourier coefficients are given by

[TABLE]

The collocation Fourier spectral first and second order derivatives of $f$ are defined as

[TABLE]

The differentiation operators in the $y$ and $z$ directions, ${\cal D}_{y}$, ${\cal D}_{y}^{2}$, ${\cal D}_{z}$ and ${\cal D}_{z}^{2}$ can be defined in the same fashion. In turn, the discrete Laplacian, gradient and divergence operators are given by

[TABLE]

at the point-wise level. It is straightforward to verify that

[TABLE]

See the derivations in the related references [5, 6, 28].

Definition 2.1**.**

Suppose that the grid function $f\in\mathcal{G}_{N}$ has the discrete Fourier expansion (2.3). Its spectral extension into the trigonometric polynomial space $\mathcal{P}_{K}$ (the space of trigonometric polynomials of degree at most $K$) is defined as

[TABLE]

We write $S_{N}(f)=f_{S}$ and call $S_{N}:\mathcal{G}_{N}\to\mathcal{P}_{K}$ the spectral interpolation operator. Suppose $g\in C_{\rm per}(\Omega,\mathbb{R})$. We define the grid projection $Q_{N}:C_{\rm per}(\Omega,\mathbb{R})\to\mathcal{G}_{N}$ via

[TABLE]

The resultant grid function may, of course, be expressed as a discrete Fourier expansion:

[TABLE]

We define the de-aliasing operator $R_{N}:C_{\rm per}(\Omega,\mathbb{R})\to\mathcal{P}_{K}$ via $R_{N}:=S_{N}(Q_{N})$. In other words,

[TABLE]

Finally, for any $g\in L^{2}(\Omega,\mathbb{R})$, we define the (standard) Fourier projection operator $P_{N}:L^{2}(\Omega,\mathbb{R})\to{\mathcal{P}}_{K}$ via

[TABLE]

where

[TABLE]

are the (standard) Fourier coefficients.

Remark 2.2**.**

Note that, in general, for $g\in C_{\rm per}(\Omega,\mathbb{R})$, $P_{N}(g)\neq R_{N}(g)$, and, in particular,

[TABLE]

However, if $g\in\mathcal{P}_{K}$ to begin with, then $\hat{g}_{\ell,m,n}=\widehat{Q_{N}(g)}_{\ell,m,n}^{N}$. In other words, $R_{N}:\mathcal{P}_{K}\to\mathcal{P}_{K}$ is the identity operator.

To overcome a key difficulty associated with the $H^{m}$ bound of the nonlinear term obtained by collocation interpolation, the following lemma is introduced.

Lemma 2.3**.**

Suppose that $m$ and $K$ are non-negative integers, and, as before, assume that $N=2K+1$. For any $\varphi\in{\cal P}_{mK}$ in $\mathbb{R}^{d}$, we have the estimate

[TABLE]

for any non-negative integer $r$.

The case of $r=0$ was proven in Weinan E’s earlier papers [19, 20]. The case of $r\geq 1$ was analyzed in a recent article by Gottlieb and Wang [30]. The proof of the following estimate can be found in [6].

Proposition 2.4**.**

Suppose that $\phi\in H_{\rm per}^{m}(\Omega)$ and $m>\frac{d}{2}$. Then, there is some constant $C>0$, such that

[TABLE]

We also note the following identity:

Proposition 2.5**.**

Suppose that $\phi\in C_{\rm per}(\Omega)$. Then

[TABLE]

Similar identities are available to the higher order derivatives.

In addition, we introduce the discrete fractional operator $(-\Delta_{N})^{\gamma}$ (with $\gamma>0$):

[TABLE]

for a grid function $f$ with the discrete Fourier expansion as (2.3). Similarly, for a grid function $f\in\mathring{\mathcal{G}}_{N}$ of (discrete) mean zero, a discrete version of the operator $(-\Delta)^{-\gamma}$ may be defined as

[TABLE]

Observe that, in this way of defining the inverse operator, the result is a periodic grid function of zero mean, i.e, $(-\Delta_{N})^{-\gamma}f\in\mathring{\mathcal{G}}_{N}$.

Detailed calculations show that the following summation-by-parts formulas are valid (see the related discussions in [7, 9, 29, 30]): for any periodic grid functions $f,g\in\mathcal{G}_{N}$,

[TABLE]

Since the SPFC equation (1.4) is an $H^{-1}$ gradient flow, we need a discrete version of the norm $\|\cdot\|_{H^{-1}}$ defined on $\mathring{\mathcal{G}}_{N}$. For any $f,g\in\mathring{\mathcal{G}}_{N}$, we define

[TABLE]

so that the $\|\cdot\|_{-1,N}$ norm could be introduced as

[TABLE]

In addition to the standard $\ell^{2}$ norm, we also introduce the $\ell^{p}$, $1\leq p<\infty$, and $\ell^{\infty}$ norms for a grid function $f\in\mathcal{G}_{N}$:

[TABLE]

The discrete $H^{1}$ and $H^{2}$ norms are introduced as

[TABLE]

For any periodic grid function $\phi\in\mathcal{G}_{N}$, the discrete SPFC energy is defined as

[TABLE]

The inequalities in the next lemma will play an important role in the convergence analysis.

Lemma 2.6**.**

For any $f\in\mathcal{G}_{N}$, we have

[TABLE]

Proof.

An application of the summation by parts formula (2.17) gives

[TABLE]

Meanwhile, another summation by parts formula reveals that

[TABLE]

Therefore, a combination of (2.24) and (2.25) leads to

[TABLE]

which in turn results in

[TABLE]

This finishes the proof of the first inequality in (2.23). For the second part, we make use of the preliminary estimate (2.25), and substitute the derived inequality in (2.27):

[TABLE]

so that the second inequality in (2.23) is valid. This completes the proof for Lemma 2.6. ∎

The following result corresponds to a discrete Sobolev embedding from $H_{N}^{2}$ to $W_{N}^{1,6}$ in the pseudo-spectral space. Similar discrete embedding estimates, in the lower order ones, could be found in Lemma 2.1 of [14]; also see the related results [25, 26] in the finite difference version. A direct calculation is not able to derive these inequalities; instead, a discrete Fourier analysis has to be applied in the derivation; the details will be left in Appendix A.

Proposition 2.7**.**

For any periodic grid function $f$, we have

[TABLE]

for some constant $C$ only dependent on $\Omega$.

3 Solvability and stability of the fully discrete schemes

A modified second order BDF temporal discretization is applied to the SPFC equation, combined with Fourier pseudo-spectral approximation in space:

[TABLE]

The discrete chemical potential is chosen in two different ways: scheme (1) ($-|\nabla\phi|^{2}$ is destabilizing)

[TABLE]

and, scheme (2) ($-\frac{\varepsilon}{2}\phi^{2}$ is destabilizing)

[TABLE]

Comparing with the standard BDF2 algorithm, the concave diffusion terms have been explicitly updated, both for the sake of unique solvability and stability. Furthermore, a second order Douglas-Dupont-type regularization term has to be added in the chemical potential to ensure energy stability. Similar ideas could be found for the epitaxial thin film growth model [26] and the Cahn-Hilliard model [54]. On the other hand, a careful analysis will reveal a subtle difference in the artificial regularization between this work and the two earlier ones, as will be demonstrated in the next section.

Remark 3.1**.**

We point out that the second scheme (3.3) is inspired by the ideas in the paper [46], where they construct an energy stable numerical scheme for the original PFC equation based on a modified Crank-Nicolson time discretization and a new convex splitting, one that views $-\frac{\varepsilon}{2}\phi^{2}$ as the destabilizing (concave) term in the energy. They show in some numerical tests that this new approach yields more accurate numerical results compared to the scheme in [36]. The scheme here is new in that we use the BDF2 approach, and it is applied to a different energy/gradient flow than the one considered in [46].

We will always assume that

[TABLE]

By $\Phi$ we denote the exact PDE solution for (1.4), and the initial value for the numerical methods is taken as

[TABLE]

We observe that (3.2) and (3.3) are two-step numerical methods, so that a “ghost” point extrapolation for $\phi^{-1}$ is useful. To preserve the second order accuracy in time, we apply the following approximation:

[TABLE]

A careful Taylor expansion indicates an $O({\Delta t}^{2}+h^{m})$ accuracy for such an approximation:

[TABLE]

Theorem 3.2**.**

Given $\phi^{k},\phi^{k-1}\in\mathcal{G}_{N}$, with $\overline{\phi^{k}}=\overline{\phi^{k-1}}$, there exist unique solutions – both labeled $\phi^{k+1}\in\mathcal{G}_{N}$, for simplicity, though they are in general distinct – for the numerical schemes (3.2) and (3.3). The schemes are mass conservative, i.e., $\overline{\phi^{k}}\equiv\overline{\phi^{0}}:=\beta_{0}$, for any $k\geq 0$.

Proof.

We prove the result only for (3.2); the proof for the other is similar. By taking a discrete summation of (3.2), and making use of the fact that $\overline{\Delta_{N}\mu_{i}^{k+1}}=0$, $i=1,2$, as well as the mass conservation of the previous time steps: $\overline{\phi^{k}}=\overline{\phi^{k-1}}=\beta_{0}$, we are able to conclude that $\overline{\phi^{k+1}}=\beta_{0}$, for any $k\geq 0$, provided a solution exists.

Next, observe that (3.2) can be rewritten as

[TABLE]

where

[TABLE]

The nonlinear equation (3.6) can be recast as a minimization problem for the following discrete energy functional:

[TABLE]

for any $\phi\in\mathcal{G}_{N}$. In turn, the strict convexity of $F_{N}$ (in terms of $\phi$), over the hyperplane of discrete functions satisfying the mass condition $\bar{\phi}=\beta_{0}$, implies a unique numerical solution for (3.2). ∎

Theorem 3.3**.**

For $k\geq 1$, define the discrete modified energies

[TABLE]

and

[TABLE]

Provided that

[TABLE]

solutions of the numerical schemes (3.2) ($i=1$) and (3.3) ($i=2$) satisfy the respective dissipation properties

[TABLE]

Proof.

Since $\phi^{k+1}-\phi^{k}\in\mathring{\mathcal{G}}_{N}$, we take a discrete inner product of (3.2) with $(-\Delta_{N})^{-1}(\phi^{k+1}-\phi^{k})$. The following identities can be derived:

[TABLE]

where the non-negative remainder $R_{4}\geq 0$ depends upon $\phi^{k+1}$, $\phi^{k}$. Furthermore,

[TABLE]

Meanwhile, an application of Cauchy inequality indicates the following estimate:

[TABLE]

A combination of (3.13) – (3.19) and (3.20) yields

[TABLE]

In addition, under the condition that

[TABLE]

we have

[TABLE]

Therefore, with an introduction of the modified discrete energy given by (3.9), we arrive at $\mathcal{E}_{N,1}^{k+1}\leq\mathcal{E}_{N,1}^{k}$, under the condition (3.11).

For the second scheme (3.3), we need the additional identities

[TABLE]

A combination of (3.13), (3.14), (3.24) – (3.26) and (3.20) yields

[TABLE]

Thus, under the condition $A\geq\frac{\varepsilon^{2}}{16}$, which is the same as for the first scheme, we have $\mathcal{E}_{N,2}^{k+1}\leq\mathcal{E}_{N,2}^{k}$. ∎

Remark 3.4**.**

We can also get a stability involving the chemical potential, if it is needed. It is observed that

[TABLE]

Thus, for example,

[TABLE]

The first term on the right-hand side was completely utilized in our stability analyses. But, at the expense of a larger splitting parameter $A$, some of the term can be spared, and a norm stability of the chemical potential can be derived.

Remark 3.5**.**

The energy stability of numerical methods for gradient flow PDE have attracted a lot of attentions over the years. For the second order numerical scheme using the BDF temporal stencil, an artificial Douglas-Dupont regularization term has to be added to ensure the energy stability, as demonstrated in recent works [26, 38, 54], for the epitaxial thin film growth and Cahn-Hilliard equations, respectively. In particular, a careful comparison reveals that, the energy growth term coming from the concave diffusion process has to be compensated by the artificial diffusion term, as well as the additional stability in the temporal stencil. This in turn requires an artificial diffusion term to have the same diffusion power as the surface diffusion term, which becomes acceptable in the practical computations.

For our first scheme (3.2), on the other hand, if we carry out the same analysis for the SPFC equation (1.4), an artificial diffusion term with the power index as $-A{\Delta t}\Delta_{N}^{3}(\phi^{k+1}-\phi^{k})$, is needed to balance the energy growth part coming from the concave diffusion. This artificial diffusion term has an even higher diffusion power than the surface diffusion term, which in turn may lead to an extra numerical dissipation in the long time simulation. To avoid such a numerical artifact, we make use of the additional stability estimate coming from the surface diffusion part, so that only an artificial diffusion power index as $-A{\Delta t}\Delta_{N}(\phi^{k+1}-\phi^{k})$ is needed in the numerical algorithm, and the energy stability could be derived using a more involved analysis. This choice avoids a higher order artificial diffusion term, therefore a reduced numerical dissipation is expected for the numerical effect.

As a direct consequence of the energy stability, a uniform in time $H_{N}^{2}$ bound for the numerical solution is given as follows.

Corollary 3.6**.**

Suppose that the initial data are sufficiently regular so that for the first scheme (3.2)

[TABLE]

and for the second scheme (3.3),

[TABLE]

for some $\tilde{C}_{0}$ that is independent of $h$, and $A\geq\frac{\varepsilon^{2}}{16}$. Then we have the following uniform (in time) $H_{h}^{2}$ bound for the numerical solution:

[TABLE]

where $\tilde{C}_{1}>0$ depends on $\Omega$ and $\tilde{C}_{0}$, but is independent of $h$, ${\Delta t}$ and final time.

Proof.

As a result of (3.12), for the first scheme (3.2), the following energy bound is available:

[TABLE]

For the second scheme (3.3)

[TABLE]

On the other hand, the point-wise quadratic inequality, $\frac{1}{4}|\nabla_{N}\phi|^{4}-2|\nabla_{N}\phi|^{2}\geq-4$, implies that

[TABLE]

Its substitution into either (3.29) or (3.30) yields

[TABLE]

so that, since $a>0$,

[TABLE]

which implies

[TABLE]

This completes the proof of Corollary 3.6. ∎

Remark 3.7**.**

As a combination of the uniform in time $H_{N}^{2}$ bound (3.28) and the discrete Sobolev embedding inequality (2.29), we arrive at a uniform in time $W_{N}^{1,6}$ estimate for the numerical solution:

[TABLE]

This estimate will be useful in the convergence analysis presented below.

4 Convergence analysis

Now we proceed into the $\ell^{\infty}(0,T;\ell^{2})\cap\ell^{2}(0,T;H_{N}^{3})$ convergence analysis for the proposed numerical scheme. With an initial data with sufficient regularity, we could assume that the exact solution has regularity of class $\mathcal{R}$:

[TABLE]

We prove convergence only for the first scheme (3.2); the details of the proof for the second are similar.

Theorem 4.1**.**

Given initial data $\Phi_{0}\in H_{\rm per}^{m+6}(\Omega)$, suppose the exact solution for SPFC equation (1.4) is of regularity class $\mathcal{R}$. Then, provided $A\geq\frac{\varepsilon^{2}}{16}$ and ${\Delta t}$ and $h$ are sufficiently small, we have

[TABLE]

where $C>0$ is independent of ${\Delta t}$ and $h$, and ${\Delta t}=T/M$.

Proof.

For $\Phi\in\mathcal{R}$, a careful consistency analysis indicates the following truncation error estimate:

[TABLE]

with $\|\tau^{k+1}\|_{2}\leq C({\Delta t}^{2}+h^{m})$. The derivation of (4.3) is accomplished with the help of the spectral approximation estimate (2.13) and other related estimates; the details are left to interested readers.

The numerical error function is defined at a point-wise level:

[TABLE]

In turn, subtracting the numerical scheme (3.2) from (4.3) gives

[TABLE]

where

[TABLE]

Taking a discrete inner product of (4.5)-(4.6) with $e^{k+1}$, with a repeated application of summation by parts, we get

[TABLE]

The temporal stencil term could be analyzed in a standard way:

[TABLE]

The artificial diffusion on the left hand side of (4.7) could be handled as follows:

[TABLE]

The term associated with the local truncation error could be bounded with the help of Cauchy inequality:

[TABLE]

The term associated with the concave diffusion could be bounded in a similar way:

[TABLE]

On the other hand, an application of the preliminary inequality (2.23) reveals that

[TABLE]

for $\ell=k-1,k,k+1$, in which the Young’s inequality has been applied at the second step. Its substitution into (4.11) yields

[TABLE]

For the nonlinear error term, we begin with the first term. An application of discrete Hölder inequality indicates that

[TABLE]

in which the preliminary estimate (3.33) (which comes from the uniform in time $H_{N}^{2}$ bound (3.28)) and (2.23) have been applied in the second and third steps, respectively. A bound for the second nonlinear error term could be derived in a similar way:

[TABLE]

in which the constant $C^{*}$ corresponds to the preliminary estimate, $\|\nabla_{N}\Phi^{k+1}\|_{6}\leq C^{*}$, for the exact solution $\Phi$. As a result, we get

[TABLE]

with $\tilde{C}_{2}=C((C^{*})^{2}+\tilde{C}_{1}^{2})$, in which the preliminary inequality (2.23) has been applied in the second step. In turn, a bound for the nonlinear error inner product term becomes available

[TABLE]

with the Young’s inequality applied in the last step.

Subsequently, a substitution of (4.8)-(4.10), (4.13) and (4.17) into (4.7) results in

[TABLE]

where

[TABLE]

Therefore, with an application of discrete Gronwall inequality, and making use of the fact that $\|\tau^{k+1}\|_{2}\leq C({\Delta t}^{2}+h^{m})$, we arrive at

[TABLE]

with $\hat{C}$ independent on ${\Delta t}$ and $h$. In turn, the desired convergence estimate is available

[TABLE]

This completes the proof of Theorem 4.1. ∎

5 Preconditioned steepest descent solver

In this section we describe a preconditioned steepest descent (PSD) algorithm following the practical and theoretical framework in [25]. We give the details for the first proposed BDF2 scheme (3.2); the details for the second (3.3) and related schemes will be similar. We first note that (3.2) can be recast as a minimization problem for (3.8). One observes that the fully discrete scheme (3.2) is the discrete variation of the strictly convex energy (3.8) set equal to zero. The nonlinear scheme at a fixed time level may be expressed as (3.6)-(3.7).

The essential idea of the PSD solver is to use a linearized version of the nonlinear operator as a pre-conditioner, or in other words, as a metric for choosing the search direction. In the following, we use the notation introduced in Theorem 3.2. A linearized version of the nonlinear operator $\mathcal{N}_{N}$, denoted as $\mathcal{L}_{N}:\mathring{\mathcal{G}}_{N}\to\mathring{\mathcal{G}}_{N}$, is defined as follows:

[TABLE]

In turn, this positive, symmetric operator could be used as a pre-conditioner for the numerical iteration. Specifically, this “metric” is used to find an appropriate search direction for the steepest descent solver [25]. Given the current iterate $\phi_{n}\in\mathcal{G}_{N}$, we define the following search direction problem: find $d_{n}\in\mathcal{G}_{N}$ such that

[TABLE]

where $r_{n}$ is the nonlinear residual of the $n^{\rm th}$ iterate $\phi_{n}$. Of course, this equation can be efficiently solved using the Fast Fourier Transform (FFT).

Subsequently, the next iterate is obtained as

[TABLE]

where $\alpha_{n}\in\mathbb{R}$ is the unique solution to the steepest descent line minimization problem

[TABLE]

Following similar techniques reported in [25] for the finite difference numerical method, a theoretical analysis ensures a geometric convergence of the iteration sequence:

[TABLE]

where $\beta$ is independent of $N$ and $\phi^{k+1}$ stands for the exact numerical solution to (3.2) at time level $k+1$, i.e., $\mathcal{N}_{N}[\phi^{k+1}]=f$.

Remark 5.1**.**

We note that our PSD method can be viewed as a quasi-Newton method, with an orthogonalization (line search) step. Indeed, $\mathcal{L}_{N}$ may be viewed as an approximation of the Jacobian. To fit more neatly into the framework of a traditional quasi-Newton method, one could just take step size equal to 1, so that the correction is just $\phi_{n+1}:=\phi_{n}+d_{n}$. Alternatively, one can just use quadratic line search methods to obtain an approximation of $\alpha_{n}$, call it $\alpha_{n}^{\rm q}$, to obtain a good-enough approximation that one can still prove a geometric convergence rate that is independent of $N$.

Remark 5.2**.**

There have been quite a few existing works [26, 43, 44, 49] on the the slope-selection (SS) model for epitaxial thin film growth, including both numerical analysis and implementation issues. These existing works only deal with the 4-Laplacian term in an $L^{2}$ gradient flow. This work is, to our knowledge, the first attempt to address the energy stability for a 4-Laplacian energy term in an $H^{-1}$ gradient flow. Because of its highly nonlinear nature, a 4-Laplacian term in an $H^{-1}$ flow turns out to be much more challenging than in an $L^{2}$ flow, at both the theoretical and numerical levels. Furthermore, most existing works for the 4-Laplacian problem have been focused on either the finite difference and finite element spatial approximations, while this work reports a Fourier pseudo-spectral approach for the first time. Such a spectral approach has the advantage of avoiding complicated staggered mesh points in the evaluation of the gradient variables, which may lead to much reduced numerical accuracy for a complicated physical energy (see a related report [24] for the functionalized Cahn-Hilliard model). Meanwhile, the analysis for such a global spatial discretization becomes much more involved, as described in this work.

6 Numerical results

6.1 Convergence test for the numerical scheme

In this subsection we perform some numerical experiments to verify the accuracy order of the proposed numerical scheme. In particular, it is observed that the search direction and Poisson-like equations can also be efficiently solved by using the Fourier pseudo-spectral method (see the related discussions in [5, 10, 28, 35]) and Fast Fourier Transform (FFT).

To test the convergence rate, we choose the following exact solution for (1.4) on the square domain $\Omega=(0,1)^{2}$:

[TABLE]

We take $a=0.975$, and we choose the artificial diffusion coefficient as $A=0.25$. The final time is taken as $T=0.16$.

To investigate the accuracy in space, we fix ${\Delta t}=10^{-4}$, and absorb the temporal discretization errors into the external source term, so that the spatial approximation error dominates the overall numerical error. We compute solutions with grid sizes $N=6$ to $N=20$ in increments of 2, and we solve up to time $T=0.16$. The $\ell^{2}$ numerical errors, computed by the proposed numerical schemes (3.2) and (3.3), are displayed in Fig. 1. The spatial spectral accuracy is apparently observed for the phase variable. Due to the round-off errors, a saturation of spectral accuracy appears with an increasing $N$, for both schemes.

To explore the temporal accuracy, we fix the spatial resolution as $N=128$ so that the numerical error is dominated by the temporal ones. We compute solutions with a sequence of time step sizes, ${\Delta t}=\frac{T}{N_{k}}$, with $N_{k}=100$ to $N_{k}=800$ in increments of 100, and the same final time $T=0.16$. Fig. 2 shows the discrete $\ell^{2}$ norms of the errors between the numerical and exact solutions, computed by the proposed numerical schemes (3.2) (BDF2-ES-1) and (3.3) (BDF2-ES-2), respectively. A clear second order accuracy has been demonstrated for both schemes. Moreover, while both schemes preserve second order temporal accuracy, it is interesting to observe that, the scale of the numerical error for the second scheme (3.3) is lower than the first scheme (3.2). Such a phenomenon comes from a subtle fact that, the second scheme (3.3) takes a closer form than the standard BDF2 scheme. On the other hand, the standard BDF2 scheme gives smaller truncation errors than the ones with extrapolation and regularization terms, although a theoretical justification of its energy stability is not available.

6.2 Numerical simulation of square symmetry patterns

The $4$-Laplacian term in (1.4) gives preference to rotationally invariant patterns with square symmetry. In this subsection, we perform two-dimensional numerical simulations showing the emergence of these patterns. The rest of the parameters are given by $a=0.5$ and $\Omega=(0,L)^{2}$, with $L=100$. The initial data for the simulations are given by

[TABLE]

where the $r_{i,j}$ are uniformly distributed random numbers in $[0,1]$. In addition, we add nucleation sites at specific locations in the domain, with magnitude 10, at (50,50) as an example of one nucleation site, at (25,25), (25,75), (75,25), and (75,75) as another example of four nucleation sites. For the temporal step size ${\Delta t}$, we use increasing values of ${\Delta t}$ in the time evolution: ${\Delta t}=0.05$ on the time interval $[0,1000]$ and ${\Delta t}=0.1$ on the time interval $[1000,9000]$. Whenever a new time step size is applied, we initiate the two-step numerical scheme by taking $\phi^{-1}=\phi^{0}$, with the initial data $\phi^{0}$ given by the final time output of the last time period. The time snapshots of the evolution by using the given parameters are presented in Figures 3 (one nucleation site) and 4 (four nucleation sites). These tests confirm the emergence of the rotationally invariant square-symmetry patterns in the density field.

To illustrate the energy stability property of the proposed numerical scheme, we display the energy evolution of the one nucleation site example, up to $t=1000$, in Figure 5. Three different time step sizes, ${\Delta t}_{1}=0.1$, ${\Delta t}_{2}=0.05$ and ${\Delta t}_{3}=0.025$, have been taken in the comparison. The energy dissipation property is satisfied in the numerical simulation. Moreover, the numerical results have a perfect agreement between different time step sizes; this is an amazing fact, due to the large time scale of the numerical simulations.

7 Concluding remarks

In this article, we proposed and analyzed two energy stable Fourier pseudo-spectral schemes for the square phase field crystal (SPFC) equation, a gradient flow to model the crystal growth. The schemes exhibit second order temporal accuracy and spectral accuracy in space. In the energy functional expansion and the corresponding PDE formulations, a composition of the 4-Laplacian and the regular Laplacian operators makes the physical system very challenging, at both the theoretical and numerical levels. To overcome this well-known difficulty, we come up with a modified BDF scheme, with the second order BDF stencil applied in the time direction, combined with an appropriate extrapolation for the concave diffusion term, to ensure the unique solvability and energy stability. In particular, a second order artificial Douglas-Dupont regularization term is added to ensure the energy stability, and a careful treatment leads to the regularization diffusion at a lower order than the surface diffusion term. Such a subtle analysis avoids a higher order artificial diffusion term, therefore a reduced numerical dissipation is expected for the numerical effect. At the theoretical side, the unique solvability, energy stability could be proved with the help of the summation-by-parts formulas in the Fourier pseudo-spectral space. In addition, the energy stability yields a uniform in time $H_{N}^{2}$ bound for the numerical solution, and an aliasing error control technique enables us to derive a discrete Sobolev embedding from $H^{2}$ into $W^{1,6}$. As a result of such a nonlinear estimate, we are able to derive an optimal rate error analysis in the $\ell^{\infty}(0,T;\ell^{2})\cap\ell^{2}(0,T;H_{N}^{3})$ norm. In the numerical implementation, the preconditioned steepest descent (PSD) iteration is needed to deal with the composition of the highly nonlinear 4-Laplacian term and the standard Laplacian term, and a geometric convergence could be proved for this iteration. It is the first such result for a 4-Laplacian solver in an $H^{-1}$ gradient flow. A few numerical experiments are presented to demonstrate the efficiency and accuracy of the proposed scheme, including the numerical accuracy test and numerical simulations of square symmetry patterns, with one nucleation site and four nucleation sites, respectively.

Acknowledgements

This work is supported in part by the Longshan Talent Project of SWUST 18LZX529 (K. Cheng), the grants NSF DMS-1418689 (C. Wang), NSF DMS-1418692 and NSF DMS-1719854 (S. Wise).

Appendix A Proof of Proposition 2.7

For simplicity of presentation, in the analysis of $\|\nabla_{N}f\|_{6}$, we are focused on the estimate of $\|{\cal D}_{x}f\|_{6}$. Due to the periodic boundary condition for $f$ and its cell-centered representation, it has a corresponding discrete Fourier transformation, as the form given by (2.3):

[TABLE]

Then we make its extension to a continuous function:

[TABLE]

The following result is excerpted as Lemma A.2 in [26]; similar analyses have also been reported in recent works [14, 25], et cetera.

Lemma A.1**.**

[26]** For $g\in\mathcal{G}_{N}$, we have

[TABLE]

Then we proceed into the proof of Proposition 2.7.

Proof.

We denote a discrete grid function, $g:={\cal D}_{x}f$, at a point-wise level. Since $f$ corresponds to $f_{N}\in{\cal B}^{K}$ (the space of trigonometric polynomials of degree at most $K$), we see that $g_{N}=\partial_{x}f_{N}$. As a result, an application of the preliminary estimate (A.3) reveals that

[TABLE]

in which the fourth step is based on the Sobolev embedding in the continuous space: $\|\partial_{x}f_{N}\|_{L^{6}}\leq C\|\Delta f_{N}\|_{L^{2}}$, and the last step comes from the fact that $f$ uniquely corresponds to $f_{N}$. Similar estimates in the other directions could be derived at the same manner:

[TABLE]

This completes the proof of Proposition 2.7. ∎