© Wuhan University 2024

0 Introduction

The epitaxial thin film growth models, as a class of the famous gradient flow type partial differential equations, have been applied to obtain high quality crystal materials. In this paper, we consider the unconditionally energy stable scheme with the variable time steps by combining the second-order backward differential formula and the scalar auxiliary variable (also called BDF2-SAV), for computing the epitaxial thin film growth models ^[1–3] in the form of

${\partial }_t\Phi =-\kappa \mu ,\mu ={ϵ}^{\mathrm{2}}{∆}^{\mathrm{2}}\Phi +\nabla \cdot g\left(\nabla \Phi \right),\mathit{x}\in \mathrm{\Omega },$(1)

in which $g\left(\mathit{u}\right)=\left(\mathrm{1}-{\left|\mathit{u}\right|}^{\mathrm{2}}\right)\mathit{u}$ for the model with slope selection, and $g\left(\mathit{u}\right)=\frac{\mathit{u}}{\mathrm{1}+{\left|\mathit{u}\right|}^{\mathrm{2}}}$ for the model without slope selection. The model parameters $\kappa$ and ${ϵ}^{\mathrm{2}}$ are the mobility constant and diffusion coefficient respectively. Meanwhile, the models can be regarded as the gradient flow of the energy functional in the form of

$E\left(\Phi \right)={\int }_{\Omega }^{}\left[\frac{{ϵ}^{\mathrm{2}}}{\mathrm{2}}{\left|∆\Phi \right|}^{\mathrm{2}}+G\left(\nabla \Phi \right)\right]dx,$(2)

where $G\left(\mathit{u}\right)=\frac{\mathrm{1}}{\mathrm{4}}{\left({\left|\mathit{u}\right|}^{\mathrm{2}}-\mathrm{1}\right)}^{\mathrm{2}}$ for the model with slope selection, and $G\left(\mathit{u}\right)=-\frac{\mathrm{1}}{\mathrm{2}}\mathrm{l}\mathrm{n}\left(\mathrm{1}+{\left|\mathit{u}\right|}^{\mathrm{2}}\right)$ for the model without slope selection.

It is well known that the epitaxial thin film growth models exist multi-scale behaviors in the time evolution, which means that the solution changes quickly in certain time intervals and slowly in others. As a highly effective technique for capturing multi-scale behavior, the adaptive time-stepping strategy is suitable for simulating the epitaxial thin film growth models. Two second order finite difference schemes ^[4] were proposed for solving the molecular beam epitaxial (MBE) model with slope selection, and proved to be unconditionally energy stable under the uniform time mesh. Based on the numerical schemes, the adaptive time-stepping algorithms were designed in Refs.[4,5] . The nonlinear BDF2 scheme with variable time steps^[6] was designed for the MBE model without slope selection, and presented the energy dissipation law under the restriction ${\tau }_n/{\tau }_{n-\mathrm{1}}<$3.561. In Ref.[7], the variable time step BDF2 scheme was also proposed for the same model, and proved to be energy stable in the modified version under the condition ${\tau }_n/{\tau }_{n-\mathrm{1}}<$ 4.864 5. The above mentioned numerical methods are all nonlinear and need the inner iteration in the practical computation, and are proved to be energy stable under certain time step constraints.

The unconditionally energy stable linear scheme is more efficient during long time evolution compared with the nonlinear scheme. The scalar auxiliary variable (SAV) approach firstly proposed in Ref.[8] is very serviceable for tackling the nonlinear term in the gradient flows. The SAV type schemes combining BDF and Crank-Nicolson were proposed for the epitaxial thin film growth models on the uniform mesh^[9]. The numerical schemes are more efficient by combining the SAV approach and the BDF formulas or adaptive time-stepping strategies^[10]. Inspired by the idea in Ref.[11], we use a first-order approximation to discretize the dynamical equation of the auxiliary variable and develop the variable time step BDF2-SAV scheme for simulating the epitaxial thin film growth models, and present the unconditional energy stability for the proposed scheme with regardless time step size.

The outline of the paper is organized as follows. Next section presents the BDF2-SAV scheme with the spatial discretization by the Fourier pseudo-spectral method. In Section 2, the unconditional energy stability is obtained by applying the discrete gradient decomposition of BDF kernels. In Section 3, the numerical experiments are carried out to test the convergence rate and effectiveness of the proposed method. At the final section, we present a brief summary.

1 The Numerical Scheme

Consider the nonuniform time levels $\mathrm{0}=t_{\mathrm{0}}<t_{\mathrm{1}}<t_{\mathrm{2}}<\cdots <t_N<T$. Denote the temporal step ${\tau }_n=t_n-t_{n-\mathrm{1}}$ and the maximum step $\tau =\underset{\mathrm{1}\le n\le N}{\mathrm{m}\mathrm{a}\mathrm{x}}{\tau }_n$. Let the adjacent step-ratio $r_n={\tau }_n/{\tau }_{n-\mathrm{1}}$ for $\mathrm{2}\le n\le N$ and $r_n=\mathrm{0}$ for $n=\mathrm{1}$. For the arbitrary real sequence $\{v^n|n=\mathrm{0,1},\mathrm{2},\cdots ,N\}$, denote the difference operator ${\nabla }_{\tau }{v}^n=v^n-v^{n-\mathrm{1}}$ and the explicit extrapolation formula ${\widehat{v}}^n=\left(\mathrm{1}+r_n\right)v^{n-\mathrm{1}}-r_n{v}^{n-\mathrm{2}}$ for $n\ge \mathrm{2}$, while ${\widehat{v}}^n=v^{\mathrm{0}}$ for $n=\mathrm{1}$. The BDF2 formula can be expressed as the convolution summation

$D_{\mathrm{2}}{v}^n={\displaystyle \sum _{i=\mathrm{1}}^n}{b}_{n-i}^{(n)}{\nabla }_{\tau }{v}^i,n\ge \mathrm{1},$(3)

in which the coefficients $b_{\mathrm{0}}^{\left(n\right)}=\frac{\mathrm{1}+\mathrm{2}{r}_n}{{\tau }_n\left(\mathrm{1}+r_n\right)}$,$b_{\mathrm{1}}^{\left(n\right)}=-\frac{r_n^{\mathrm{2}}}{{\tau }_n\left(\mathrm{1}+r_n\right)}$, $b_j^{\left(n\right)}=\mathrm{0}$, $\mathrm{2}\le j\le n-\mathrm{1}$, are also called the discrete BDF kernels. Specifically, $D_{\mathrm{2}}{v}^{\mathrm{1}}=\frac{\mathrm{1}}{{\tau }_{\mathrm{1}}}{\nabla }_{\tau }{v}^{\mathrm{1}}=D_{\mathrm{1}}{v}^{\mathrm{1}}$ is the BDF1 formula.

We consider the problem (1) equipped with the periodic boundary conditions. Let the domain $\mathrm{\Omega }={\left(\mathrm{0},L\right)}^{\mathrm{2}}$ and divide it by ${\mathrm{\Omega }}_h=\{{\mathit{x}}_h=\left(ih,jh\right)|\mathrm{1}\le i,j\le M(\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{g}\mathrm{e}\mathrm{r})\}$ with the space step $h=L/M$. Define the space of grid functions ${{V}}_h=\{v|v=\left(v_h\right)$ is $\left(L,L\right)$-period for ${\mathit{x}}_h\in {\overline{\mathrm{\Omega }}}_h={\mathrm{\Omega }}_h\cup \partial {\mathrm{\Omega }}_h\}$ and the function $\phi \left(\mathit{x}\right)={\mathrm{e}}^{\omega \mathrm{i}\left(mx+ny\right)}$ with $\omega =\mathrm{2}\mathrm{\pi }/L$. Suppose ${\mathrm{ℱ}}_M$ is the space containing all trigonometric polynomials of degree up to $M/\mathrm{2}$. Let $I_M:L^{\mathrm{2}}\left(\mathrm{\Omega }\right)\to {\mathrm{ℱ}}_M$ be the trigonometric interpolation operator, namely, $\left(I_{M}u\right)\left(\mathit{x}\right)={\displaystyle \sum _{m,n=-M/\mathrm{2}}^{M/\mathrm{2}-\mathrm{1}}}{\tilde{u}}_{m,n}\phi \left(\mathit{x}\right)$, in which the coefficients ${\tilde{u}}_{m,n}$ are determined by $\left(I_{M}u\right)\left({\mathit{x}}_{\mathit{h}}\right)=u_h$. The $p$-th order pseudo-spectral derivatives of $u_h$ are given as

$\begin{array}{l}{D}{}_x^p{u}_h={\displaystyle \sum _{m,n=-M/\mathrm{2}}^{M/\mathrm{2}-\mathrm{1}}}{\left(\omega m\mathrm{i}\right)}^p{\tilde{u}}_{m,n}\phi \left({\mathit{x}}_{\mathit{h}}\right),\\ {D}{}_y^p{u}_h={\displaystyle \sum _{m,n=-M/\mathrm{2}}^{M/\mathrm{2}-\mathrm{1}}}{\left(\omega n\mathrm{i}\right)}^p{\tilde{u}}_{m,n}\phi \left({\mathit{x}}_{\mathit{h}}\right).\end{array}$

The discrete gradient and Laplacian operators are defined by ${\nabla }_h{u}_h={\left({D}{}_x^{\mathrm{1}}{u}_h,{D}{}_y^{\mathrm{1}}{u}_h\right)}^{\mathrm{T}}$, ${∆}_h{u}_h={\nabla }_h\cdot \left({\nabla }_h{u}_h\right)={D}{}_x^{\mathrm{2}}{u}_h+{D}{}_y^{\mathrm{2}}{u}_h$. For any grid functions $v,w\in {{V}}_h$, one defines the discrete inner product $〈v,w〉=h^{\mathrm{2}}{\displaystyle \sum _{{\mathit{x}}_h\in {\mathrm{\Omega }}_h}^{}}{v}_h{w}_h$ and the corresponding $L^{\mathrm{2}}$ norm $\Vert v\Vert =\sqrt[]{〈v,v〉}$. Furthermore, the discrete Green's formula ^[12] gives $〈{∆}_h^{\mathrm{2}}u,v〉=〈{∆}_{h}u,{∆}_{h}v〉$, for any $u,v\in {{V}}_h$.

We denote $E_{\mathrm{1}}\left(\Phi \right)={\int }_{\Omega }^{}\left[\frac{\alpha }{\mathrm{2}}{\left|∆\Phi \right|}^{\mathrm{2}}+G\left(\nabla \Phi \right)\right]d\mathit{x}$ with $\mathrm{0}\le \alpha <{ϵ}^{\mathrm{2}}$. It follows from Ref. [9] that $E_{\mathrm{1}}\left(\Phi \right)$ has a lower bound. One selects a scalar auxiliary variable defined by $\gamma \left(t\right)=\sqrt[]{E_{\mathrm{1}}\left(\Phi \right)+C_{\mathrm{0}}}$, where $C_{\mathrm{0}}$ is a positive constant to guarantee that the term under the root sign is positive. Then, the total energy of the models becomes

$\overline{E}\left(\Phi ,\gamma \right)={\int }_{\Omega }^{}\frac{{ϵ}^{\mathrm{2}}-\alpha }{\mathrm{2}}{\left|∆\Phi \right|}^{\mathrm{2}}d\mathit{x}+{\gamma }^{\mathrm{2}}-C_{\mathrm{0}}.$(4)

Denote $\mathrm{\mathscr{H}}\left(\Phi \right)=\alpha {∆}^{\mathrm{2}}\Phi +\nabla \cdot g\left(\nabla \Phi \right)$ and $\xi =\frac{\gamma \left(t\right)}{\sqrt[]{E_{\mathrm{1}}\left(\Phi \right)+C_{\mathrm{0}}}}$. The models (1) are rewritten as the expanding system, namely

${\partial }_t\Phi =-\kappa \mu ,\mu =\left({ϵ}^{\mathrm{2}}-\alpha \right){∆}^{\mathrm{2}}\Phi +\left(\mathrm{2}-\xi \right)\xi \mathrm{\mathscr{H}}\left(\Phi \right),$(5)

${\gamma }^{\text{'}}=\frac{\mathrm{2}-\xi }{\mathrm{2}\sqrt[]{E_{\mathrm{1}}\left(\Phi \right)+C_{\mathrm{0}}}}{\int }_{\Omega }^{}\mathit{\mathscr{H}}\left(\Phi \right){\partial }_t\Phi d\mathit{x}.$(6)

It is natural to check that the dissipative law holds, that is $\frac{\mathrm{d}}{\mathrm{d}t}\overline{E}\left(\Phi ,\gamma \right)=-\kappa {\Vert \nabla \mu \Vert }_{L^{\mathrm{2}}}^{\mathrm{2}}\le \mathrm{0}$.

Based on the SAV approach combining the BDF formula, the variable time step energy stable numerical scheme is designed for solving the epitaxial thin film growth models. That is, finding ${\varphi }^n\in {{V}}_h$ and ${\gamma }^n\in \mathbb{R}$ such that

$D_{\mathrm{2}}{\varphi }^n=-\kappa \left[\left({ϵ}^{\mathrm{2}}-\alpha \right){\mathrm{\Delta }}_h^{\mathrm{2}}{\varphi }^n+\left(\mathrm{2}-{\xi }^n\right){\xi }^n\mathrm{\mathscr{H}}\left({\widehat{\varphi }}^n\right)\right],\mathrm{1}\le n\le N,$(7)

$\{\begin{array}{c}{\nabla }_{\tau }{\gamma }^n=\frac{\mathrm{2}-{\xi }^n}{\mathrm{2}\sqrt[]{E_{\mathrm{1}}\left({\varphi }^{n-\mathrm{1}}\right)+C_{\mathrm{0}}}}〈\mathrm{\mathscr{H}}\left({\widehat{\varphi }}^n\right),{\nabla }_{\tau }{\varphi }^n〉,\\ {\xi }^n=\frac{{\gamma }^n}{\sqrt[]{E_{\mathrm{1}}\left({\varphi }^{n-\mathrm{1}}\right)+C_{\mathrm{0}}}},\mathrm{1}\le n\le N.\end{array}$(8)

Remark 1   For the proposed numerical scheme, the scalar auxiliary variable $\gamma$ with the first-order approximation is different from that in Ref. [9]. Meanwhile, it holds $\left|\mathrm{1}-\left(\mathrm{2}-{\xi }^n\right){\xi }^n\right|=O\left({\tau }^{\mathrm{2}}\right)$, which ensures that the convergence order of ${\varphi }^n$ can arrive at second order.

Next, we show how to carry out the new proposed numerical scheme. It follows from (7) that

$\begin{array}{l}{b}_{\mathrm{0}}^{\left(n\right)}{\varphi }^n+\kappa \left({ϵ}^{\mathrm{2}}-\alpha \right){\mathrm{\Delta }}_h^{\mathrm{2}}{\varphi }^n\\ =\left(b_{\mathrm{0}}^{\left(n\right)}-b_{\mathrm{1}}^{\left(n\right)}\right){\varphi }^{n-\mathrm{1}}+b_{\mathrm{1}}^{\left(n\right)}{\varphi }^{n-\mathrm{2}}-\kappa \left(\mathrm{2}-{\xi }^n\right){\xi }^n\mathrm{\mathscr{H}}\left({\widehat{\varphi }}^n\right).\end{array}$(9)

Let ${\varphi }_{\mathrm{1}}^n$ and ${\varphi }_{\mathrm{2}}^n$ be the solution of the following two linear systems respectively,

$\left[b_{\mathrm{0}}^{\left(n\right)}+\kappa \left({ϵ}^{\mathrm{2}}-\alpha \right){\mathrm{\Delta }}_h^{\mathrm{2}}\right]{\varphi }_{\mathrm{1}}^n=\left(b_{\mathrm{0}}^{\left(n\right)}-b_{\mathrm{1}}^{\left(n\right)}\right){\varphi }^{n-\mathrm{1}}+b_{\mathrm{1}}^{\left(n\right)}{\varphi }^{n-\mathrm{2}},$(10)

$\left[b_{\mathrm{0}}^{\left(n\right)}+\kappa \left({ϵ}^{\mathrm{2}}-\alpha \right){\mathrm{\Delta }}_h^{\mathrm{2}}\right]{\varphi }_{\mathrm{2}}^n=-\kappa \mathrm{\mathscr{H}}\left({\widehat{\varphi }}^n\right).$(11)

Define ${\varphi }^n={\varphi }_{\mathrm{1}}^n+\left(\mathrm{2}-{\xi }^n\right){\xi }^n{\varphi }_{\mathrm{2}}^n$. One can check that it is the solution of the linear system (9). In fact, one solves (10) and (11) to get ${\varphi }_{\mathrm{1}}^n$ and ${\varphi }_{\mathrm{2}}^n$ respectively. It remains to calculate the value of ${\xi }^n$. By substituting the definition of ${\varphi }^n$ into (8), one has a nonlinear algebraic equation about ${\xi }^n$:

$\mathrm{2}{\xi }^n\left[E_{\mathrm{1}}\left({\varphi }^{n-\mathrm{1}}\right)+C_{\mathrm{0}}\right]-\mathrm{2}{\gamma }^{n-\mathrm{1}}\sqrt[]{E_{\mathrm{1}}\left({\varphi }^{n-\mathrm{1}}\right)+C_{\mathrm{0}}}$

$-\left(\mathrm{2}-{\xi }^n\right)〈\mathrm{\mathscr{H}}\left({\widehat{\varphi }}^n\right),{\varphi }_{\mathrm{1}}^n〉-{\left(\mathrm{2}-{\xi }^n\right)}^{\mathrm{2}}{\xi }^n〈\mathrm{\mathscr{H}}\left({\widehat{\varphi }}^n\right),{\varphi }_{\mathrm{2}}^n〉=\mathrm{0}.$(12)

One solves the above equation by using the Newton's iteration to get the value of ${\xi }^n$. When ${\varphi }_{\mathrm{1}}^n$, ${\varphi }_{\mathrm{2}}^n$and${\xi }^n$ are known, the solutions of ${\varphi }^n$ and ${\gamma }^n$ will be obtained right now.

2 Energy Stability

We develop the unconditional energy stability for the numerical scheme (7)-(8) by using the discrete gradient decomposition of the BDF2 formula. Actually, for any real sequence $\{w_n{\}}_{n=\mathrm{1}}^N$, a gradient structure was proposed in Ref.[13] as follows

$\mathrm{2}{w}_n{\displaystyle \sum _{i=\mathrm{1}}^n}{b}_{n-i}^{(n)}{w}_i\ge \frac{r_{n+\mathrm{1}}^{\frac{\mathrm{3}}{\mathrm{2}}}}{\mathrm{1}+r_{n+\mathrm{1}}}\frac{w_n^{\mathrm{2}}}{{\tau }_n}-\frac{r_n^{\frac{\mathrm{3}}{\mathrm{2}}}}{\mathrm{1}+r_n}\frac{w_{n-\mathrm{1}}^{\mathrm{2}}}{{\tau }_{n-\mathrm{1}}}+R_L(r_n,r_{n+\mathrm{1}})\frac{w_n^{\mathrm{2}}}{{\tau }_n}$(13)

where the time step ratio is set to be $\mathrm{0}<r_n<$4.864 for 2$\le n\le N$, and the binary function $R_L\left(z,s\right)=\frac{\mathrm{2}+\mathrm{4}z-z^{\frac{\mathrm{3}}{\mathrm{2}}}}{\mathrm{1}+z}-\frac{s^{\frac{\mathrm{3}}{\mathrm{2}}}}{\mathrm{1}+s}>\mathrm{0}$, $\mathrm{0}<z,s<$4.864. Based on the structure, we give the energy decay law under no time step restriction.

Theorem 1   The variable time step BDF2-SAV scheme is unconditionally energy stable in modified version at the discrete levels. In details, define the modified discrete energy formula

$\mathrm{ℰ}\left[{\varphi }^n,{\gamma }^n\right]=\overline{E}\left[{\varphi }^n,{\gamma }^n\right]+\frac{r_{n+\mathrm{1}}^{\frac{\mathrm{3}}{\mathrm{2}}}}{\mathrm{1}+r_{n+\mathrm{1}}}\frac{{\Vert {\nabla }_{\tau }{\varphi }^n\Vert }^{\mathrm{2}}}{\kappa {\tau }_n},$(14)

in which $\overline{E}\left[{\varphi }^n,{\gamma }^n\right]=\frac{{ϵ}^{\mathrm{2}}-\alpha }{\mathrm{2}}{\Vert {∆}_h{\varphi }^n\Vert }^{\mathrm{2}}+{\left({\gamma }^n\right)}^{\mathrm{2}}-C_{\mathrm{0}}$. For $\mathrm{0}<r_n<$4.864 with 2$\le n\le N$, it holds

$\mathrm{ℰ}\left[{\varphi }^n,{\gamma }^n\right]\le \mathrm{ℰ}\left[{\varphi }^{n-\mathrm{1}},{\gamma }^{n-\mathrm{1}}\right].$(15)

Proof   Taking the discrete inner product on both sides of the equality (6) with ${\nabla }_{\tau }{\varphi }^n/\kappa$, multiplying the expansion equation (8) by $\mathrm{2}{\gamma }^n$, it gives that

$\begin{array}{l}\frac{\mathrm{1}}{\kappa }〈D_{\mathrm{2}}{\varphi }^n,{\nabla }_{\tau }{\varphi }^n〉+\left({ϵ}^{\mathrm{2}}-\alpha \right)〈{\mathrm{\Delta }}_h^{\mathrm{2}}{\varphi }^n,{\nabla }_{\tau }{\varphi }^n〉\\ =-〈\left(\mathrm{2}-{\xi }^n\right){\xi }^n\mathrm{\mathscr{H}}\left({\widehat{\varphi }}^n\right),{\nabla }_{\tau }{\varphi }^n〉,\end{array}$(16)

$\mathrm{2}{\gamma }^n{\nabla }_{\tau }{\gamma }^n=\frac{\left(\mathrm{2}-{\xi }^n\right){\gamma }^n}{\sqrt[]{E_{\mathrm{1}}\left({\varphi }^{n-\mathrm{1}}\right)+C_{\mathrm{0}}}}〈\mathrm{\mathscr{H}}\left({\widehat{\varphi }}^n\right),{\nabla }_{\tau }{\varphi }^n〉.$(17)

Summing up the equalities (16) and (17), it follows that

$\frac{\mathrm{1}}{\kappa }〈D_{\mathrm{2}}{\varphi }^n,{\nabla }_{\tau }{\varphi }^n〉+\left({ϵ}^{\mathrm{2}}-\alpha \right)〈{∆}_h{\varphi }^n,{\nabla }_{\tau }{∆}_h{\varphi }^n〉+\mathrm{2}{\gamma }^n{\nabla }_{\tau }{\gamma }^n=\mathrm{0},$(18)

which indicates that there is no need to estimate the nonlinear term in the numerical scheme.

To handle the term related to BDF2 formula, we use the gradient structure (13) to obtain

$\frac{\mathrm{1}}{\kappa }〈D_{\mathrm{2}}{\varphi }^n,{\nabla }_{\tau }{\varphi }^n〉\ge \frac{r_{n+\mathrm{1}}^{\frac{\mathrm{3}}{\mathrm{2}}}}{\mathrm{1}+r_{n+\mathrm{1}}}\frac{{\Vert {\nabla }_{\tau }{\varphi }^n\Vert }^{\mathrm{2}}}{\kappa {\tau }_n}-\frac{r_n^{\frac{\mathrm{3}}{\mathrm{2}}}}{\mathrm{1}+r_n}\frac{{\Vert {\nabla }_{\tau }{\varphi }^{n-\mathrm{1}}\Vert }^{\mathrm{2}}}{\kappa {\tau }_{n-\mathrm{1}}},$(19)

where the time step ratio satisfies $\mathrm{0}<r_n<$4.864. Furthermore, by applying the elementary equality $\mathrm{2}a\left(a-b\right)=a^{\mathrm{2}}-b^{\mathrm{2}}+{\left(a-b\right)}^{\mathrm{2}}$ and the Green's formula, it holds that

$\begin{array}{l}({ϵ}^{\mathrm{2}}-\alpha )〈{∆}_h{\varphi }^n,{\nabla }_{\tau }{∆}_h{\varphi }^n〉\\ =\frac{{ϵ}^{\mathrm{2}}-\alpha }{\mathrm{2}}({\Vert {∆}_h{\varphi }^n\Vert }^{\mathrm{2}}-{\Vert {∆}_h{\varphi }^{n-\mathrm{1}}\Vert }^{\mathrm{2}}+{\Vert {∆}_h{\nabla }_{\tau }{\varphi }^n\Vert }^{\mathrm{2}})\end{array}$(20)

$\mathrm{2}{\gamma }^n{\nabla }_{\tau }{\gamma }^n={\left({\gamma }^n\right)}^{\mathrm{2}}-{\left({\gamma }^{n-\mathrm{1}}\right)}^{\mathrm{2}}+{\left({\nabla }_{\tau }{\gamma }^n\right)}^{\mathrm{2}}.$(21)

Substituting the estimates (19)-(21) into the equality (18), it leads to the claimed energy dissipation law (15) at the discrete levels. This completes the proof.

3 Numerical Experiments

We focus on the temporal accuracy of the numerical method (7)-(8) on random time meshes. Let time step size ${\tau }_n=T{\sigma }_n/S$ for $\mathrm{1}\le n\le N$, in which $S={\displaystyle \sum _{n=\mathrm{1}}^N}{\sigma }_n$ and ${\sigma }_n\in \left(\mathrm{0,1}\right)$ is the uniformly distributed random number. The experiment temporal convergence order is computed by

$\mathrm{O}\mathrm{r}\mathrm{d}\mathrm{e}{\mathrm{r}}_{\varphi }=\mathrm{l}\mathrm{o}{\mathrm{g}}_{\mathrm{2}}\left(e_{\varphi }\left(N\right)/e_{\varphi }\left(\mathrm{2}N\right)\right)/\mathrm{l}\mathrm{o}{\mathrm{g}}_{\mathrm{2}}\left(\tau \left(N\right)/\tau \left(\mathrm{2}N\right)\right),$

where the discrete $L^{\mathrm{2}}$ norm error is recorded by $e_{\varphi }\left(N\right)=\Vert \Phi \left(T\right)-{\varphi }^N\Vert$, and $\tau \left(N\right)$ denotes the maximum time step size for total $N$ subintervals. Also, we test the order of SAV function $\gamma \left(t\right)$ in the way of

$\mathrm{O}\mathrm{r}\mathrm{d}\mathrm{e}{\mathrm{r}}_{\gamma }=\mathrm{l}\mathrm{o}{\mathrm{g}}_{\mathrm{2}}\left(e_{\gamma }\left(N\right)/e_{\gamma }\left(\mathrm{2}N\right)\right)/\mathrm{l}\mathrm{o}{\mathrm{g}}_{\mathrm{2}}\left(\tau \left(N\right)/\tau \left(\mathrm{2}N\right)\right),$

in which the corresponding error is defined by $e_{\gamma }\left(N\right)=\underset{\mathrm{1}\le n\le N}{\mathrm{m}\mathrm{a}\mathrm{x}}\left|\gamma \left(t_n\right)-{\gamma }^n\right|$.

To verify the numerical accuracy, one computes the models with an artificial forcing term $f\left(\mathit{x},t\right)$ and the exact solution $\Phi \left(\mathit{x},t\right)=\mathrm{c}\mathrm{o}\mathrm{s}t\mathrm{s}\mathrm{i}\mathrm{n}x\mathrm{s}\mathrm{i}\mathrm{n}y$, namely

${\partial }_t\Phi +\kappa \left[{ϵ}^{\mathrm{2}}{∆}^{\mathrm{2}}\Phi +\nabla \cdot g\left(\nabla \Phi \right)\right]=f\left(\mathit{x},t\right),\mathit{x}\in {\left(\mathrm{0,2}\mathrm{\pi }\right)}^{\mathrm{2}},$(22)

where the parameters $\kappa =\mathrm{0.5},{ϵ}^{\mathrm{2}}=\mathrm{0.5}$. Consider the above model with slope selection. We solve it until time $T=\mathrm{6}$ by using the numerical scheme (7)-(8) with the parameters $\alpha =\mathrm{0},C_{\mathrm{0}}=\mathrm{1},$$M=\mathrm{128}$. The errors and convergence orders are presented in Table 1. Again, we compute the model (22) with $g\left(\nabla \Phi \right)=\frac{\nabla \Phi }{\mathrm{1}+{\left|\nabla \Phi \right|}^{\mathrm{2}}}$ and $T=\mathrm{1}$ by applying the proposed scheme with the parameters $\alpha =\mathrm{0.1}{ϵ}^{\mathrm{2}},C_{\mathrm{0}}=\mathrm{2},M=\mathrm{256}$. The related numerical results are listed in Table 2. The maximum step ratio ($\mathrm{m}\mathrm{a}\mathrm{x}{r}_n$) and the number ($N_{\mathrm{1}}$) of time levels with the step ratios $r_n\ge \mathrm{4.864}$ are also listed in Tables 1-2. From the tables, we find that the first-order convergence of the SAV function ${\gamma }^n$ does not affect the second-order convergence of the numerical solution ${\varphi }^n$, and the adjacent time step ratio can exceed 4.864 in practical calculation.

We simulate the epitaxial thin film growth models (1) with the initial value $\Phi \left(\mathit{x},\mathrm{0}\right)=\mathrm{0.1}[\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{3}x\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{2}y+\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{5}x$

$\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{5}y$, the domain $\mathrm{\Omega }={\left(\mathrm{0,2}\mathrm{\pi }\right)}^{\mathrm{2}}$ and the parameters $\kappa =\mathrm{1},{ϵ}^{\mathrm{2}}=\mathrm{0.1}$. We run the numerical schemes (7)-(8) with $M=\mathrm{128}$ by applying the following adaptive time-stepping strategy^[14]

${\tau }_{\mathrm{a}\mathrm{d}\mathrm{a}}=\mathrm{m}\mathrm{a}\mathrm{x}\{{\tau }_{\mathrm{m}\mathrm{i}\mathrm{n}},\frac{{\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\sqrt[]{\mathrm{1}+\lambda {\Vert {\partial }_{\tau }{\varphi }^n\Vert }^{\mathrm{2}}}}\},\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{n}{\tau }_{n+\mathrm{1}}=\mathrm{m}\mathrm{i}\mathrm{n}\{{\tau }_{\mathrm{a}\mathrm{d}\mathrm{a}},r_{\mathrm{u}\mathrm{s}\mathrm{e}\mathrm{r}}{\tau }_n\},$(23)

in which parameter $\lambda >\mathrm{0}$ need be selected by the user, ${\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ and ${\tau }_{\mathrm{m}\mathrm{i}\mathrm{n}}$ are the pre-selected maximum and minimum time steps, respectively. Here, we set $\lambda =\mathrm{1}\mathrm{000},{\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}={\mathrm{10}}^{-\mathrm{1}},{\tau }_{\mathrm{m}\mathrm{i}\mathrm{n}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$. The numerical energy curves and the corresponding adaptive time steps are summarized respectively for the epitaxial thin film growth models with or without slope selection, see Figs. 1 and 2. It is found that the time evolutions of the total energy are consistent with that in Ref.[9] and the proposed numerical scheme combining the adaptive time-stepping algorithm performs well for the present numerical simulations.

Fig. 1 The total energy and the related time step sizes for the numerical solution of the model with slope selection by using the scheme (7)-(8) with ${\mathit{C}}_{\mathrm{0}}=\mathrm{2},$$\mathit{\alpha }=\mathrm{0}$

Fig. 2 The total energy and the related time step sizes for the numerical solution of the model without slope selection by using the scheme (7)-(8) with ${\mathit{C}}_{\mathrm{0}}=\mathrm{2.75},$$\mathit{\alpha }=\mathrm{0}$

4 Conclusion

In this article, the BDF2-SAV scheme with the variable time steps is designed for the epitaxial thin film growth models with or without slope selection. The unconditional energy stability in modified version is proved without any time step constraints. The effectiveness and accuracy of the proposed scheme are demonstrated by the numerical experiments.

Note that the SAV scheme (7)-(8) is a coupled system, and the auxiliary variable makes the nonlinearity term more complicated. Therefore, it can be foreseen that the error estimate may be very sophisticated. We will focus on the $L^{\mathrm{2}}$ norm convergence analysis in the future work.

References

All Tables

All Figures

	Fig. 1 The total energy and the related time step sizes for the numerical solution of the model with slope selection by using the scheme (7)-(8) with ${\mathit{C}}_{\mathrm{0}}=\mathrm{2},$$\mathit{\alpha }=\mathrm{0}$
In the text

	Fig. 2 The total energy and the related time step sizes for the numerical solution of the model without slope selection by using the scheme (7)-(8) with ${\mathit{C}}_{\mathrm{0}}=\mathrm{2.75},$$\mathit{\alpha }=\mathrm{0}$
In the text