The operator splitting method with semi-implicit spectral deferred correction for molecular beam epitaxial growth models

Abstract

This paper presents a high order time discretization method by combining the time splitting method with semi-implicit spectral deferred correction method to simulate the molecular beam epitaxial growth models with and without slope selection. The original problem is split into linear and nonlinear subproblems. The Fourier spectral method is adopted for the linear part, and a second-order SSP-RK method together with the Fourier spectral method is used for the nonlinear part. The scheme takes advantage of avoiding nonlinear iteration. However, the temporal error is dominated by the operator splitting error, which is second order for Strang splitting. In order to achieve higher order numerical algorithm in time, we consider a semi-implicit spectral deferred correction (SDC) method to reduce the splitting error. Specifically, the temporal order is increased by one with each correction loop in the SDC framework. Numerical results are given to illustrate that the high order temporal algorithm is a practical, accurate and efficient simulation tool for molecular beam epitaxial growth models.

Introduction

The molecular beam epitaxy (MBE) equation is a fourth-order nonlinear parabolic equation. Epitaxy is referred to the deposition of a crystalline overlayer on a substrate. MBE model is a high vacuum technique for the growth of epitaxial layers, usually semiconductors, that utilizes thermal beams of source atoms or molecules impinging on a single crystal substrate. It can produce many complex structures of varying layers that are further processed to produce a range of electronic and optoelectronic devices, including high speed transistors, high-efficiency solar cells, light-emitting diodes and solid state lasers. In the recent years, MBE equation becomes an important and challenging research topic in the material science which motivates people to develop mathematical models as well as simulation tools to understand the primary mechanisms behind the technology.

In this work, we consider a broadly-used continuum model for the MBE, which is derived via a L²-gradient flow associated with the energy

$$E(u)={\displaystyle {\int }_{\Omega }\left(\frac{{\mathrm{ϵ}}^2}{2}|\Delta u{|}^2+F(\nabla u)\right)}\text{d}\mathbf{x}$$

(1)

where Ω is a smooth domain in R^d (d = 1, 2, 3), u: $\Omega \to R^d$ is a scaled height function of a thin film in a co-moving frame, $F(\nabla u)$ is a smooth function of its argument $\nabla u$ , and ϵ is a constant.

There are two popular choices of the nonlinear bulk potential:

(I) the logarithmic potential for the case without slope selection model

$$F(\nabla u)=-\frac{1}{2}\ln (|\nabla u{|}^2+1)$$

(2)

(II) the double well potential for the case of slope selection model

$$F(\nabla u)=\frac{1}{4}{(|\nabla u{|}^2-1)}^2$$

(3)

Due to the presence of the negative term, the energy $E(\varphi )$ with (2) is expected to behave badly, and is in fact poorly understood mathematically. The $E(\varphi )$ with (3) appears in several areas of material modeling.^1–3

The $L^2-$ gradient flow of (1) with (2) and (3) for the MBE models are respectively given without slope selection model

$$u_t=-\frac{\delta E}{\delta u}=-{\mathrm{ϵ}}^2{\Delta }^{2}u-\nabla ·\left(\frac{\nabla u}{1+|\nabla u{|}^2}\right)$$

(4)

and with slope selection model

$$u_t=-\frac{\delta E}{\delta u}=-{\mathrm{ϵ}}^2{\Delta }^{2}u+\nabla ·\left((|\nabla u{|}^2-1)\nabla u\right)$$

(5)

The initial condition of both models is

$$u(\mathbf{x},0)=u_0(\mathbf{x}),\mathbf{x}\in \Omega$$

The periodic boundary condition is adopted here since it is usually used to remove all boundary integrals. As for the well-posedness and solution regularity of the initial-boundary-value problems of the above models, see literature.^4–6

In the context of MBE, equation (3) selects the slope of the MBE surface and thus equation (5) is called the MBE model with slope selection, and correspondingly equation (4) is called the MBE model without slope selection. The interfacial dynamics governed by equations (4) and (5) are different. Solutions to equation (5) exhibit pyramidal structures, where the faces of the pyramids have slopes $|\nabla u|\hspace{0.17em}\approx \hspace{0.17em}1$ (see, e.g., literature ⁷ ) However, solutions to model (4) exhibit mound-like structures, the slopes of which (on an infinite domain) may grow unbounded (see, e.g., literature ⁸ )

The models (4) and (5) satisfy energy dissipative law and mass conservation

$$\frac{d}{dt}E(u)=-{\displaystyle {\int }_{\Omega }{(u_t)}^2}\text{d}\mathbf{x}\le 0,\frac{d}{dt}{\displaystyle {\int }_{\Omega }u}(\mathbf{x},t)\mathrm{dx}=0$$

There have been many numerical studies on the MBE model: the convex splitting method Ju et al. 2017,^3,9 the stabilization method Xu and Tang 2006, ⁶ the invariant energy quadratization (IEQ) method, ¹⁰ the scalar auxiliary variable (SAV) method11 adaptive time-stepping algorithms,^12,13 exponential time difference (ETD) scheme ¹⁴ and so on. The operator splitting strategy is also a very effective numerical method to solve MBE model.^15,16 The basic idea of the operator splitting method ¹⁷ is that the original problem is divied into subproblems which are simpler than the original problem. Then the approximate solution of the original problem is composed using the exact or approximate solutions of the subproblems in a given sequential order.

The spectral deferred correction (SDC) was first introduced by Dutt et al. ¹⁸ to solve the Cauchy problem for ODEs. The key idea of the SDC method is to convert the original ODEs into the corresponding Picard equation and then apply a deferred correction procedure in the integral formulation, driven by either the explicit or the implicit Euler marching scheme. An advantage of SDC method is that it is a one step method and possess favorable accuracy and stability properties even for versions with very high order of accuracy.

The motivation of this study is to an second-order operator splitting Fourier method combined with SDC method for MBE model. The MBE equation equation with or without slope selection was first divided into linear fourth order parabolic equation and nonlinear equation. The linear part was solved using Fourier spectral spatial discretization, and the nonlinear part was solved by the spectral method combined with a second-order SSP-RK method in time. Then, a semi-implicit spectral deferred correction (SDC) method is employed to improve the temporal accuracy. It is obvious that we can achieve high order accuracy in both time and space. Finally, we present numerical experiments demonstrating the convergence order and numerical properties of the proposed methods.

The outline of the paper is as follows. In Section 2, a fast explicit operator splitting spectral method for problem (1) is developed. In Section 3, a high order semi-implicit SDC method combining with the operator splitting spectral algorithm is introduced. In Section 4, the numerical results confirming the accuracy and efficiency of the obtained method are presented. Finally, conclusions are drawn in Section 5.

Fast explicit operator splitting method for MBE model

In this section, we present a fast explicit operator splitting method for MBE models (4) and (5). The proposed method is based on the operator splitting method that combines spatial discretization by Fourier spectral method with temporal discretization by the SSP-RK method. We restrict the description to two space dimension (2 D) for notational simplicity.

The operator splitting method for MBE models (4) and (5) are constructed as follows. Writing both equations in the following abstract initial value problem

$$u_t=L(u)+N(\nabla u)$$

(6)

where

$$\begin{array}{l}L(u)=-{\mathrm{\epsilon }}^2{\Delta }^{2}u\\ N(\nabla u)=\{\begin{array}{l}-\nabla ·\left(\frac{\nabla u}{1+|\nabla u{|}^2}\right),\left\{\frac{\text{Eq}.(\text{4})\text{without}}{\text{slope}\text{selection}.}\right\}\\ \nabla ·\left((|\nabla u{|}^2-1)\nabla u\right),\left\{\frac{\text{Eq}.(\text{5})\text{ with }}{\text{slope}\text{selection}.}\right\}\end{array}\end{array}$$

Splitting the original models into two subequations

$$u_t=L(u)$$

(7)

$$u_t=N(\nabla u)$$

(8)

Then given a time step τ, the second-order Strang splitting method ¹⁷ reads as

$$u(x,y,t+\tau )={\mathcal{S}}^L\left(\frac{\tau }{2}\right){\mathcal{S}}^N(\tau ){\mathcal{S}}^L\left(\frac{\tau }{2}\right)u(x,y,t)+O({\tau }^3)$$

(9)

where ${\mathcal{S}}^L$ and ${\mathcal{S}}^N$ denote the exact operators associated with equations (7) and (8), respectively.

Introducing the following spatial grid ${\Omega }_h^{\mathit{per}}=\left\{(x_j,y_k)=\left(a+jh,a+kh\right),0\le j,k\le N-1)\right\}$ for the periodic boundary condition, where $h=(b-a)/N$ . All of the 2 D periodic grid functions defined on ${\Omega }_h^{\mathit{per}}$ are denoted by ${\mathcal{M}}^N$ . Let $u_{jk}^m\approx u(x_j,y_k,m\tau )$ , where $\tau =T/M$ is the time step and M is the total number of time steps. For a function $u\in {\mathcal{M}}^N$ , the 2D discrete Fourier transform $\widehat{u}=Pu$ is defined componentwise ¹⁹ by

$$\begin{array}{l}{\mathcal{F}}_{\mathcal{N}}[u(t)](p,q):={\widehat{u}}_{pq}(t)\\ =\frac{h^2}{C_p{C}_q{(b-a)}^2}{\displaystyle \sum _{j=0}^{N-1}{\displaystyle \sum _{k=0}^{N-1}{u}_{jk}}(t){\phi }_{pq}(x_j,y_k)}\end{array}$$

(10)

The function u can be reconstructed via the corresponding inverse transform $u=P^{-1}\widehat{u}$ with components given by

$${\mathcal{F}}_N^{-1}[\widehat{u}(t)](j,k):=u_{j,k}(t)={\displaystyle \sum _{p=-N/2}^{N/2}{\displaystyle \sum _{q=-N/2}^{N/2}{\widehat{u}}_{pq}}(t){\phi }_{pq}(x_j,y_k)}$$

(11)

where c_p, c_q and ${\phi }_{p,q}$ are defined as

$$\begin{array}{l}{\phi }_{p,q}=\exp \left(i\frac{2p\pi (x-a)}{b-a}+i\frac{2q\pi (y-a)}{b-a}\right)\\ C_r=\{\begin{array}{l}2,\left|r\right|=\frac{N}{2},\hfill \\ 1,\left|r\right|<\frac{N}{2},\hfill \end{array}r=p\text{or}q\end{array}$$

Let $\widehat{{\mathcal{M}}^N}=\left\{Pu\right|u\in {\mathcal{M}}^N\}$ and define the operators $\widehat{D_x}$ and $\widehat{D_y}$ on $\widehat{{\mathcal{M}}^N}$ as

$${\left({\widehat{D}}_x\widehat{u}\right)}_{pq}=\frac{2p\pi i}{b-a}{\widehat{u}}_{pq},{\left({\widehat{D}}_y\widehat{u}\right)}_{pq}=\frac{2q\pi i}{b-a}{\widehat{u}}_{pq}$$

then the Fourier spectral approximations to the first order partial derivatives can be represented as

$$D_x=P^{-1}{\widehat{D}}_{x}P,D_y=P^{-1}{\widehat{D}}_{y}P$$

For any $u^1,u^2\in {\mathcal{M}}^N,\mathbf{u}=(u^1,u^2)\in {\mathcal{M}}^N\times {\mathcal{M}}^N$ , the discrete gradient and divergence are given respectively by

$${\nabla }_N{u}^1={(D_x{u}^1,D_y{u}^1)}^T,{\nabla }_N·\mathbf{u}=D_x{u}^1+D_y{u}^2$$

(12)

Let $X_N:=\text{span}\{{\phi }_{pq}(x,y):-N/2\le p,q\le N/2\}$ and $P_N:L^2(\Omega )\to X_N$ be the L²-orthogonal projection, defined by

$$(P_{N}u-u,v)=0,\forall v\in X_N$$

Thus

$$\left(P_{N}u\right)(x,y,t)={\displaystyle \sum _{p=-N/2}^{N/2}{\displaystyle \sum _{q=-N/2}^{N/2}{\widehat{u}}_{pq}}}(t){\phi }_{pq}(x,y)$$

(13)

is the truncated series. The L² estimate of the projection errors between $P_{N}u$ and u in fractional Sobolev spaces is stated as follows:

Theorem 1. ¹⁹ For any $u\in H_{\mathit{Per}}^s(\Omega )=\{u\in H^s(\Omega )|u\text{is}\Omega \text{-periodic}\}$ with s > 1, it holds that

$$\parallel P_{N}u-u{\parallel }_2\le N^{-s}|u{|}_s$$

where $\parallel ·{\parallel }_m$ and $|·{|}_m$ represents the seminorm of $H_{\mathit{Per}}^m(\Omega ) (m\le s)$ .

Now we give the numerical approximations ${\mathcal{S}}_h^L(\tau )$ and ${\mathcal{S}}_h^N(\tau )$ .

Step 1: ${\mathcal{S}}^L(\tau )$ → ${\mathcal{S}}_h^L(\tau )$ .

We first consider equation (7) subject to periodic boundary condition from t_m to $t_{m+1}$ with an initial condition u^m. Based on Fourier spectral method, equation (7) can be solved in the following analytical formula

$$u^{m+1,(1)}={\mathcal{F}}_N^{-1}\left\{e^{-\tau {\mathrm{ϵ}}^2{\lambda }_{pq}^2}{\mathcal{F}}_N[u^m]\left(p,q\right)\right\}$$

(14)

where ${\lambda }_{pq}={\left(\frac{2p\pi }{b-a}\right)}^2+{\left(\frac{2q\pi }{b-a}\right)}^2,p,q=0,\pm 1,\pm 2,\pm \frac{N}{2}.$

By Parseval’s formula, it yields

$$\parallel u^{m+1}\parallel \le \parallel u^m\parallel$$

(15)

where $\parallel ·\parallel$ represents the discrete L²-norm, that is

$$\parallel u^m\parallel =\sqrt{h^2{\displaystyle \sum _{i=0}^{N-1}{\displaystyle \sum _{j=0}^{N-1}{u}_{ij}^m}}}$$

Hence, we have the following lemma.

Lemma 1. For any grid function $u=\left\{u_{ij}\right|0\le i,j\le N-1\}$ , it holds that

$$\parallel {\mathcal{S}}_h^L(\tau )u\parallel \le \parallel u\parallel$$

Step 2: ${\mathcal{S}}^N(\tau )$ → ${\mathcal{S}}_h^N(\tau )$ .

We now consider the nonlinear equation (8) subject to periodic boundary conditions. Based on the second order SSP-RK method, ²⁰ equation (8) can be solved by

$$\{\begin{array}{l}{u}^{*}=u^m+\tau N(\nabla u^m)\hfill \\ u^{m+1}=\frac{1}{2}{u}^m+\frac{1}{2}{u}^{*}+\frac{\tau }{2}N(\nabla u^{*})\hfill \end{array}$$

(16)

Define a mapping $\beta :{ℝ}^2\to {ℝ}^2$ as

$$\beta (\upsilon )=\{\begin{array}{l}-\frac{\upsilon }{1+|\upsilon {|}^2},\text{equation}(\text{4})\text{without}\text{slope}\text{selection}\hfill \\ (|\upsilon {|}^2-1)\upsilon ,\text{equation}(\text{5})\text{ with }\text{slope}\text{selection}\hfill \end{array}$$

(17)

Taking Fourier transform in both sides of the scheme (16) and combining with equation (12), we have the fully discrete scheme

$$\{\begin{array}{l}{u}^{*}=u^m+\tau {\nabla }_N·\beta ({\nabla }_N{u}^m)\hfill \\ u^{m+1}=\frac{1}{2}{u}^m+\frac{1}{2}{u}^{*}+\frac{\tau }{2}{\nabla }_N·\beta ({\nabla }_N{u}^{*})\hfill \end{array}$$

(18)

Now we are in the position to give a second-order fast explicit operator splitting method. The whole procedures from t_m to $t_{m+1}$ are stated as follows:

$$\mathbf{Algorithms}1\{\begin{array}{l}{u}^{m+1,(1)}={\mathcal{F}}_N^{-1}\left\{e^{-\frac{\tau }{2}{\mathrm{ϵ}}^2{\lambda }_{pq}^2}{\mathcal{F}}_N[u^m](p,q)\right\}\hfill \\ u^{*}=u^m+\tau {\nabla }_N·\beta ({\nabla }_N{u}^m)\hfill \\ u^{m+1,(2)}=\frac{1}{2}{u}^m+\frac{1}{2}{u}^{*}+\frac{\tau }{2}{\nabla }_N·\beta ({\nabla }_N{u}^{*})\hfill \\ u^{m+1}={\mathcal{F}}_N^{-1}\left\{e^{-\frac{\tau }{2}{\mathrm{ϵ}}^2{\lambda }_{pq}^2}{\mathcal{F}}_N[u^{m+1,(2)}](p,q)\right\}\hfill \end{array}$$

SDC method for improving the temporal order of accuracy

The Algorithms 1 takes advantage of avoiding nonlinear iteration and reaching spectral accuracy in space. However, the temporal error is dominated by the operator splitting error. In order to achieve higher order numerical algorithm in time, we consider a semi-implicit SDC method to reduce the splitting error.

Suppose the time interval $[0,T]$ has been divided into M non-overlapping intervals by the partition $0=t_0<t_1<\cdots <t_n<\cdots <t_M=T$ . We shall describe below the semi-implicit SDC method which will be used over one time step $[t_n,t_{n+1}]$ with length $\Delta t_n=t_{n+1}-t_n$ . $u_n$ denotes the numerical approximation of $u(t_n)$ , with $u_0=u(t_0)$ .

A single time step of an SDC method begins by first dividing the time step $[t_n,t_{n+1}]$ into a set of intermediate sub-steps defined by P + 1 Legendre-Gauss-Radau nodes (including the two endpoints) with $t_n=t_{n,0}<t_{n,1}<\cdots <t_{n,m}\cdots <t_{n,P}=t_{n+1}$ . Let $\Delta t_{n,m}=t_{n,m+1}-t_{n,m}$ and $u_{n,m}^k$ denotes the k-th order approximation to $u(t_{n,m})$ . The procedure of an SDC method with P + 1 distributed quadrature nodes and with K correction loops is the following:

Algorithms 2

step 1. (prediction step) Let $u_{n,0}^{(0)}=u^n$ , then use Algorithm 1 to obtain a numerical solution $u_{n,1}^{(0)},u_{n,1}^{(0)},\dots ,u_{n,P}^{(0)}$ , which is a low order approximation to the exact solution at Legendre-Gauss-Radau nodes (including the two endpoints).

step 2. (correction loops) For $k=0,\dots ,K-1$

$$u_{n,0}^{k+1}=u_{n,0}^{(0)}$$

For $m=0,1,\dots ,P-1$

$$\begin{array}{l}{u}_{n,m+1}^{k+1}=u_{n,m}^{k+1}+\Delta t_{n,m}\left(L(u_{n,m+1}^{k+1})-L(u_{n,m+1}^k)\right)\\ +{\Pi }_{n,m}^{n,m+1}\left(L(u^k)+N(u^k)\right)\end{array}$$

where ${\Pi }_{n,m}^{n,m+1}\left(L(u^k)+N(u^k)\right)$ is the integral of the P-th degree interpolating polynomial on the P + 1 points $\left(t_{n,m},L(u_{n,m}^k)+N(u_{n,m}^k)\right)$ over the subinterval $[t_{n,m},t_{n,m+1}]$ .

step 3. Finally we have $u^{n+1}=u_{n,P}^K$ .

Remark 1 ¹⁸ The local truncation error obtained with the Algorithms 2 is $O({\tau }^{\min (K+2,P+1)})$ .

Numerical experiments

In this section, we perform a series of numerical experiments to demonstrate the accuracy and efficiency of the obtained numerical algorithm.

Problem 1: Convergence test

In order to show the order of convergence of the obtained method, we start with the study of the 2D case on $\Omega =[-5,5]\times [-5,5]$ with $\mathrm{ϵ}=0.5$ and the initial condition is as follows

$$u_0(x,y)=0.01\exp \left(-5(x^2+y^2)\right)$$

Since both models (4) and (5) exhibit the same behavior in a short time, we only consider the convergence of our method for model (4) in this problem.

We first consider the convergence order both in terms of L²- and $L^{\infty }$ -norm. However, the exact solution for MBE model is unknown, to this end, we define the error and convergence rate in the temporal direction as follows

$$\begin{array}{ll}{\text{Err}}_{L^{\infty }}^M=\parallel U^{M,N}-U^{2M,N}{\parallel }_{\infty },& {\text{order}}_{L^{\infty }}={\log }_2({\text{Err}}_{L^{\infty }}^M/{\text{Err}}_{L^{\infty }}^{2M})\\ {\text{Err}}_{L^2}^M=\parallel U^{M,N}-U^{2M,N}\parallel ,& {\text{order}}_{L^2}={\log }_2({\text{Err}}_{L^2}^M/{\text{Err}}_{L^2}^{2M})\end{array}$$

where $U^{M,N}=\left\{u_{i,j}^{M,N}\right|1\le i,j\le N\}$ as the numerical solution of time grid $\tau =\frac{T}{M}$ and space grid $h=\frac{b-a}{N}$ at time T.

The computational results are presented in Table 1 for the MBE model without slope selection. One may see that the numerical orders of time accuracy is close to the optimal $(K+2)$ -order.

Table 1.

Temporal accuracy for MBE without slope selection at T = 5 with N = 128 and P = 3 for Problem 1.

K	$\tau =$	$\frac{1}{5}$	$\frac{1}{10}$	$\frac{1}{20}$	$\frac{1}{40}$	$\frac{1}{80}$
0	${\text{Err}}_{L^{\infty }}^{\tau }$	1.25e-2	4.32e-3	1.25e-3	3.34e-4	8.60e-5
	${\text{order}}_{L^{\infty }}$	–	1.53	1.79	1.90	1.96
	${\text{Err}}_{L^2}^{\tau }$	2.63e-2	9.04e-3	2.61e-3	6.95e-4	1.79e-4
	${\text{order}}_{L^2}$	–	1.54	1.79	1.91	1.96
1	${\text{Err}}_{L^{\infty }}^{\tau }$	8.66e-4	1.34e-4	1.89e-5	2.55e-6	3.34e-7
	${\text{order}}_{L^{\infty }}$	–	2.69	2.83	2.89	2.93
	${\text{Err}}_{L^2}^{\tau }$	1.84e-3	2.83e-4	3.93e-5	5.19e-6	6.69e-7
	${\text{order}}_{L^2}$	–	2.70	2.85	2.92	2.96
2	${\text{Err}}_{L^{\infty }}^{\tau }$	1.93e-4	1.67e-5	1.37e-6	1.12e-7	8.96e-9
	${\text{order}}_{L^{\infty }}$	–	3.53	3.61	3.61	3.64
	${\text{Err}}_{L^2}^{\tau }$	4.19e-4	3.58e-5	2.71e-6	2.08e-7	1.54e-8
	${\text{order}}_{L^2}$	–	3.55	3.73	3.70	3.76

Secondly, we test the spatial errors by letting N vary and fixing K = 2 and M = 2000, which is large enough to avoid contamination of the temporal error. The error in the spatial direction as follows

$$\begin{array}{l}{\text{Err}}_{L^{\infty }}^N=\parallel U^{M,N}-U^{M,2N}(2:2:2N,2:2:2N){\parallel }_{\infty }\hfill \\ {\text{Err}}_{L^2}^N=\parallel U^{M,N}-U^{M,2N}(2:2:2N,2:2:2N){\parallel }_0\hfill \end{array}$$

The $L^{\infty }$ and L² errors at T = 5 for the MBE model without slope selection are shown in Figure 1 on a loglog scale. It is clear from the figure that our algorithm is exponentially convergent.

Figure 1.

Spatial L² and $L^{\infty }$ errors at time T = 5 for MBE model without slope selection for Problem 1.

Problem 2: Growth dynamics

This example is to study the long time dynamical behavior of the MBE models (4) and (5). We take the same initial condition and boundary condition as in Problem 1. The simulation parameters are chosen as follows

$$\mathrm{ϵ}=0.1,N=128,T=400,\tau =0.01$$

The numerical results are shown in Figures 2 and 3. We see from Figure 2 that solutions to model (4) exhibit mound-like structures. There are many small hills (red part) and valleys (blue part) at the early period, while at the final time t = 400, the system saturates to a one-hill-one-valley structure. However, the interfacial dynamics in governed by equation (5) is different, see the Figure 3. The morphology of the growing interface is characterized by the development of pyramidal structures. These are induced by the presence of the slope selection in the model. The pyramidal structures are made of nearly flat facets that meet at sharp edges.

Figure 2.

The numerical solutions at t = 10,50,100,200,400 for the MBE models (4).

Figure 3.

The numerical solutions at t = 10,50,100,200,400 for the MBE models (5). (a) Mass for MBE without slope. (b) Energy for MBE without slope. (c) Mass for MBE with slope. (d) Energy for MBE with slope.

The time evolution of the energy E(t), and mass M(t) are presented in Figure 4, which show that the proposed method is mass conserving and energy decreasing when solving the MBE models.

Figure 4.

The curves of mass and energy for Problem 2.

Problem 3: Coarsening dynamics

For the no-slope-selection MBE model (4), it is shown in Golubovic ⁸ that

$$E(t)\hspace{0.17em}\sim \hspace{0.17em}O(-\ln t),R(t)\hspace{0.17em}\sim \hspace{0.17em}O(t^{1/2})$$

However, for the slope-selection MBE model (5), it is shown in Moldovan and Golubović ⁷ that

$$E(t)\hspace{0.17em}\sim \hspace{0.17em}O(t^{-1/3}),R(t)\hspace{0.17em}\sim \hspace{0.17em}O(t^{1/3})$$

Here the surface roughness R(t) is defined by

$$R(t)=\sqrt{\frac{1}{\Omega }{\displaystyle {\int }_{\Omega }{\left(\varphi (\mathbf{x},t)-\overline{\varphi }(t)\right)}^2}\text{d}\mathbf{x}},\overline{\varphi }=\frac{1}{\Omega }{\displaystyle {\int }_{\Omega }\varphi }(\mathbf{x},t)\text{d}\mathbf{x}$$

In order to study the energy decay rate and the surface roughness growth rate, we take an initial condition as $\varphi (x,y,0)=\mathit{rand}(x,y)$ on the domain $\Omega =[0,12.8]\times [0,12.8]$ . The rand(x, y) is a random number varying from –0.001 to 0.001. The simulation parameters are chosen as follows

$$\mathrm{ϵ}=0.1,N=128,T=600,\tau =0.01$$

In Figures 5 and 6, the time evolution of the energy and roughness for models (4) and (5) are plotted respectively. For the MBE model (4), Figure 5(b) shows the linear fitting of the semi-log energy data up to t = 600, where the fitting line is $E(t)=-40\ln (t)-20$ . Figure 5(d) shows the linear fitting of the log-log surface roughness data up to t = 600, where the fitting line is $R(t)=0.4t^{1/2}$ . For the MBE model (5), Figure 6(b) and (d) show the linear fitting of the log-log surface energy and roughness data up to t = 600, where the fitting lines are $E(t)=25t^{-1/3}$ and $R(t)=0.15t^{1/3}$ , respectively.

Figure 5.

Time evolution of the energy (first line) and roughness (second line) for the MBE model without slope selection.

Figure 6.

Time evolution of the energy (first line) and roughness (second line) for the MBE model with slope selection.

Conclusions

In this work, by a combination of the operator splitting method and SDC method, a high order operator splitting Fourier spectral scheme for the MBE model with and without slope selection were developed. The semi-implicit SDC method was employed to improve the temporal accuracy. Numerical experiments were presented to confirm the accuracy and efficiency of the method. Additionally, during the simulations for the coarsening process, it is clearly observed that the energy decay and roughness growth rates. In future work, we plan to to analyze the convergence rate of the method theoretically.