 # A stabilized second order exponential time differencing multistep method for thin film growth model without slope selection

In this paper, a stabilized second order in time accurate linear exponential time differencing (ETD) scheme for the no-slope-selection thin film growth model is presented. An artificial stabilizing term Aτ^2∂Δ^2 u/∂ t is added to the physical model to achieve energy stability, with ETD-based multi-step approximations and Fourier collocation spectral method applied in the time integral and spatial discretization of the evolution equation, respectively. Long time energy stability and detailed ℓ^∞(0,T; ℓ^2) error analysis are provided based on the energy method, with a careful estimate of the aliasing error. In addition, numerical experiments are presented to demonstrate the energy decay and convergence rate.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Consider the continuum model of no-slope-selection thin film epitaxial growth, which is the gradient flow of the following energy functional:

 E(u)=∫Ωε22|Δu|2−12ln(1+|∇u|2) dx, (1.1)

where , is a scaled height function of thin film with periodic boundary condition, and is a constant. Here the first term in represents the surface diffusion, and the logarithm term models the Ehrlich-Schowoebel (ES) effect which describes the effect of kinetic asymmetry in the adatom attachment-detachment; see [13, 26, 27, 37]. Consequently, the growth equation is the gradient flow of (1.1):

 ∂u∂t=−ε2Δ2u−∇⋅(∇u1+|∇u|2), x∈Ω, t∈(0,T]. (1.2)

Taking the inner product with (1.2) by , we obtain the energy decay property for the continuous system:

 dEdt=−∥ut∥2L2.

The equation (1.2) is referred to as the no-slope-selection (NSS) equation, since (1.2) predicts an unbounded mound slope [16, 27]. On the other hand, the slope-selection (SS) case refers to the case when the logarithm term in (1.2) is replaced by , which predicts that mound-like or pyramid structures in the surface profile tend to have a uniform, constant mound slope .

The solution is mass conservative, i.e.,

 ∫Ωu(x,t) dx=∫Ωu(x,0) dx, t>0, (1.3)

due to the periodic boundary condition. For simplicity of presentation, herein we assume to have zero mean-value. Scaling laws that characterize the surface coarsening during the film growth have been a physically interesting problem; see [16, 25, 26, 28, 33]. As for the continuum model, the well-posedness of the initial-boundary-value problem for both SS and NSS equation has been given by [26, 27] through the perturbation analysis and Galerkin spectral approximations. Li and Liu  proved that the minimum energy scales as for small , and the average energy is bounded below by for large . This implies that long time energy stability is required to simulate the coarsening process.

One popular way to construct energy stable numerical scheme is to split the energy functional into convex and concave parts, then apply implicit and explicit time stepping algorithms to the corresponding terms , respectively; see the first such numerical scheme in  for the molecular epitaxy growth model. Since then, there have been various works dedicated to deriving high order and energy stable schemes under this idea, such as [8, 9, 10, 32, 35, 36, 39]. In particular, a first order in time linear scheme was proposed in  to solve the NSS equation, with the energy stability established. On the other hand, the concave term turns out to be nonlinear in this approach, and there has been a well-known difficulty of the convex-splitting approach to construct unconditionally stable higher-order schemes for a nonlinear concave term. Many efforts have been devoted to overcome this subtle difficulty, such as introducing an artificial stabilizing term in the growth equation to balance the explicit treatment of the nonlinear term; see [15, 30, 31, 43]. In addition, there have been some other interesting approaches for the stability of a numerically modified energy functional, such as the invariant energy quadratization (IEQ)  and the scalar auxiliary variable (SAV) methods .

Other than these direct approaches to establish the energy stability, some other ideas have been reported in recent years to obtain higher order temporal accuracy with explicit treatment of nonlinear terms, such as the time exponential time differencing (ETD) approach. In general, an exact integration of the linear part of the NSS equation is involved in the ETD-based scheme, followed by multi-step explicit approximation of the temporal integral of the nonlinear term [3, 4, 12, 19, 20]. An application of such an idea to various gradient models has been reported in recent works [11, 21, 22, 23, 24, 42, 45], with the high order accuracy and preservation of the exponential behavior observed in the numerical experiments. At the theoretical side, some related convergence analyses have also been reported, while a rigorous energy stability analysis has only been provided for the first order scheme . For the second and higher order numerical schemes, a theoretical justification of the energy bound has not been available.

In this paper, we work on a second order in time accurate ETD-based scheme, with Fourier pseudo-spectral approximation in space, and the energy stability analysis is established at a theoretical level. We use an alternate approach to overcome the above-mentioned difficulty, instead of applying convex splitting to the energy functional (1.1). In more details, an artificial stabilizing term is added in the growth equation, where is a positive constant and is the time step size. Also, we apply a linear Lagrange approximation to the nonlinear term as in . This approach enables us to perform a careful energy estimate, so that a decay property for a modified discrete energy functional is proved. This in turn leads to a uniform in time bound for the original energy functional. In our knowledge, this is the first such result for a second order ETD-based scheme, with the primary difficulty focused on the explicit extrapolation for the nonlinear terms. Moreover, we provide a novel convergence analysis for the proposed scheme. Instead of analyzing the operator form of the numerical error function, here we start from the continuous in time ODE system satisfied by the error function, which is derived from the corresponding equations of the numerical solution. With a careful treatment of the aliasing error and estimate of the numerical solution, we are able to derive an error estimate for the proposed scheme.

There have been quite a few physically interesting quantities that may be obtained from the solutions of the NSS equation (1.2), such as the energy, average surface roughness and the average slope. A theoretical analysis in  has implied a lower bound for the energy dissipation as the order of , and an upper bound for the average roughness as the order of , which in turn implies an order for the average slope. In particular, the energy dissipation scale leads to a long time scale nature of the NSS equation (1.2), in comparison with the scale for the slope-selection version and the standard Cahn-Hilliard model. In addition, a lower bound for the energy, estimated as as derived in [8, 9, 41], indicates an intuitive law for the saturation time scale for , because of the energy dissipation scale. All these facts have demonstrated the necessity of energy stable numerical approach to accurately capture the physically interesting quantities, since the energy stability turns out to be a global in time nature for the numerical scheme. In the extensive numerical experiments, the proposed stabilized ETD scheme has produced highly accurate solutions, which obtain the long time asymptotic growth rate of the surface roughness and average slope with relative accuracy within .

The rest of this article is organized as follows. In Section 2 we present the fully discrete numerical scheme. The numerical energy stability is proved in Section 3, followed by and bounds of the numerical solution. Subsequently, an error analysis is given by Section 4, consisting of two lemmas concerning the error of the nonlinear term. In addition, numerical experiments are provided in Section 5, including temporal convergence test and simulation results of the scaling laws for energy, average surface roughness and average slope. Finally, some concluding remarks are given by Section 6.

## 2 The stabilized second order in time ETD mulstistep scheme (sEDTMs2)

Some definitions in  are recalled. Consider . We define , where is a -tuple of non-negative integers with , and . For simplicity, we denote as , and the related semi-norm as . In particular, if , we set as , and as .

In this paper, we focus on the case when and follow the notations used in collocation spectral method; see [5, 6, 17, 18, 21, 38]. Let be a positive integer, be a uniform mesh on , with grid points , where , with , . Denote as the exact solution of (1.2), be the restriction of the exact solution on , be the time step size , and for . Let , be the space of 2D periodic grid functions on , and be the space of trigonometric polynomials in and of degree up to . In this paper, we denote one generic constant which may depend on , the solution , the initial value and time , but is independent of the mesh size and time step size .

Let us firstly recall the following Berstein inequality introduced in [38, p. 33, Lemma 2.5]:

###### Lemma 2.1

For any and , we have

 ∥∂mxu∥Lp(Ω)≤CNm∥u∥Lp(Ω),m≥1. (2.1)

Now we introduce an interpolation operator

onto that reserves the function value on grid points, i.e., for :

 (INf)(x,y)=N∑k,l=−N+1(^fc)k,lexp(2πiL(kx+ly)),with i=√−1, (2.2)

where the coefficients

are given by the discrete Fourier transform of the

grid points:

 (^fc)k,l=14N22N∑i,j=1f(xi,yj)exp(−2πiL(kxi+lyj)). (2.3)

For any , denote as the continuous extension of . When and () are continuous and periodic on , has the following approximation property ([7, Theorem 1.2, p. 72]):

 ∥∂kf−∂kINf∥L2≤Chm−k∥f∥Hm, for 0≤k≤m, m>d2, (2.4)

for dimension . Also, the following bound of is excerpted in [18, Lemma 1]:

###### Lemma 2.2

For any in dimension , we have

 ∥INφ∥Hk≤(√2)d∥φ∥Hk, k≥0. (2.5)

Given , the discrete spatial partial derivatives can be defined as:

 (Dxf)i,j =N∑k,l=−N+12kπiL(^fc)k,lexp(2πiL(kxi+lyj)), (D2xf)i,j =N∑k,l=−N+1−4k2π2L2(^fc)k,lexp(2πiL(kxi+lyj)).

Similarly we have and . For any , and , we can define the discrete gradient, divergence and Laplace operators:

 ∇Nf=(DxfDyf),∇N⋅f=Dxf1+Dyf2,ΔNf=D2xf+D2yf.

Also, to measure the discrete differentiation operators defined above, we introduce the discrete (denoted as ) inner product and norm :

 (f,g)N=h22N∑i,j=1fijgij,∥f∥N=√(f,f)N, (f,g)N=h22N∑i,j=1(f1ijg1ij+f2ijg2ij),∥g∥N=√(f,f)N.

Similarly, we can define the discrete Sobolev norm and the discrete Sobolev semi-norm :

 ∥f∥H2h=(∑|α|≤2∥Dαf∥2N)12,|f|H2h=(∑|α|=2∥Dαf∥2N)12,

where as above is a -tuple of nonnegative integers with , and Furthermore, the following summation by parts formulas are available in [21, proposition 2.2]:

###### Lemma 2.3

For any , and , we have the following summation by parts formula:

 (f,∇N⋅g)N=−(∇Nf,g)N,(f,ΔNg)N=−(∇Nf,∇Ng)N=(ΔNf,g)N.

Recall that the exact solution in (1.2) is assumed to be mean value free, thus we consider the subset of zero-mean grid functions:

 MN0={v∈MN∣(v,1)N=0}={v∈MN∣^v00=0}.

Define , which is symmetric positive definite on .

In turn, the spatial discretization of (1.2) becomes: Given , find such that

 d^udt=−LN^u−fN(^u),LN=ε2Δ2N,t∈(0,T], (2.6) ^u(0)=u(0),

where . Multiplying both sides of (2.6) by yields

 deLNt^udt=−eLNtfN(^u). (2.7)

Integrating (2.7) from to gives

 ^u(tn+1)=e−LNτ^u(tn)−∫tn+1tne−LN(tn+1−t)fN(^u(t))dt. (2.8)

By [21, Lemma 4.1], the error between the exact solution and the solution of (2.8) is of order , given that is smooth enough.

### 2.1 The ETD1 and ETDMs2 schemes

Two exponential time differencing schemes (ETD1 and ETDMs2) have been proposed in , using the energy convex splitting method. Ju et al used a similar spatial discretization as (2.8) is derived, with modified and :

 Lc =LN−κΔN,κ>0, (2.9) fe(^u(t)) =fN(^u(t))+κΔN^u(t). (2.10)

For the ETD1, the term in is simply approximated by . For the ETDMs2, is approximated by the linear Lagrange interpolation:

 fe(^u(tn))+t−tnτ[fe(^u(tn))−fe(^u(tn−1))], t∈[tn,tn+1].

Let be the numerical solution of ETD1 and ETDMs2. Denote as for . Integrating from to , the following expressions are obtained in :

• ETD1:

 un+1h=e−τLcunh−ϕ0(Lc)fe(unh).
• ETDMs2:

 un+1h=e−τLcunh−ϕ0(Lc)fe(unh)−ϕ1(Lc)[fe(unh)−fe(un−1h)],

in which

 ϕ0(x)=x−1(I−e−xτ), ϕ1(x)=x−1[1−(xτ)−1(I−e−xτ)].

The convergence analysis for both ETD1 and ETDMs2 was given by . However, a theoretical proof of the long time energy stability of ETDMs2 was not provided.

### 2.2 The sETDMs2 scheme

In this section we propose a second order in time stabilized exponential time differencing multistep scheme (sETDMs2). In order to guarantee the long time energy stability, a stabilizing term is introduced, in which is a positive constant. Also, we apply the same Lagrange approximation for as in . Define a continuous in time function for :

 fN,L(s,^u(tn),^u(tn−1)):=fN(^u(tn))+sτ[fN(^u(tn))−fN(^u(tn−1))].

We present the following sETDMs2 scheme:

• sETDMs2: For , find such that for any ,

 dun+1s(t)dt+Aτ2% dΔ2Nun+1s(t)dt=−LNun+1s(t)−fN,L(t−tn,uns,un−1s), (2.11)

with for .

The initial step is obtained as follow: Find such that

 du1s(t)dt+Aτ2dΔ2Nu1s(t)dt=−LNu1s(t)−fN(u0s), t∈(0,t1], (2.12)

where . In fact, as for the initial step, the stabilizing term can also be replaced by . Integrating (2.11) and (2.12) from to , one obtains explicit expressions of and for sETDMs2:

• sETDMs2 (matrix form): For ,

 un+1s=e−KNτuns−ϕ0(KN)GN(uns)−ϕ1(KN)[GN(uns)−GN(un−1s)], (2.13)

and the intial step is computed by

 u1s=e−KNτu0s−ϕ0(KN)GN(u0s), (2.14)

in which

 KN=ε2(I+Aτ2Δ2N)−1Δ2N, GN(v)=(I+Aτ2Δ2N)−1fN(v), ϕ0(KN)=K−1N(I−e−KNτ), ϕ1(KN)=K−1N[I−(KNτ)−1(I−e−KNτ)].

Note that sETDMs2 and ETDMs2 have similar matrix forms for solutions at , and they share similar computational complexity.

## 3 Long time energy stability of the sETDMs2 scheme

In later analysis we shall make frequent use of the Gagliardo-Nirenberg inequality and Agmon’s inequality (see [1, 2, 34]): for function in ,

 |v|Wj,p(Ω) (3.1) ∥v∥L∞(Ω) ≤C∥v∥1−n2mL2(Ω)∥v∥n2mHm(Ω),m>d2, v∈Hm(Ω), (3.2)

with .

Consider the following numerical energy functional:

 ~E(un+1s(t)) =EN(un+1s(t))+√3−14∥∥ ∥∥dun+1s(t)dt∥∥ ∥∥2L2(tn,tn+1;ℓ2) +(√3+1)τ224∥∥ ∥∥dΔNun+1s(t)dt∥∥ ∥∥2L2(tn,tn+1;ℓ2), (3.3)

where is the discrete version of the continuous energy functional in (1.1):

 EN(v)=(−12ln(1+|∇Nv|2),1)N+ε22∥ΔNv∥2N,   ∀v∈MN. (3.4)

The main theoretical result of this section is the following theorem.

###### Theorem 3.1

Assume that . The numerical system (2.11) is energy stable in the sense that for any ,

 ~E(un+1s)+(A−2+√36)τ2∥∥ ∥∥dΔNun+1s(t)dt∥∥ ∥∥2L2(tn,tn+1;ℓ2)≤~E(uns).

We define and denote

 ^fN,L(t−tn,uns,un−1s):=β(∇Nuns)+t−tnτ[β(∇Nuns)−β(∇Nun−1s)].

By definition we have: For any and ,

 fN(v)=∇N⋅β(∇Nv), fN,L(ζ,v,w)=∇N⋅^fN,L(ζ,v,w).

Taking an inner product with (2.11) by gives that: For ,

 = (^fN,L(t−tn,uns,un−1s),d∇Nun+1s(t)dt)N. (3.5)

Recall that , thus

 ddtEN(un+1s(t))=(−β(∇Nun+1s(t)),d∇Nun+1s(t)dt)N+ε22d∥ΔNun+1s(t)∥2Ndt.

Therefore, we can rewrite (3.5) as below:

 = (^fN,L(t−tn,uns,un−1s)−β(∇Nun+1s(t)),d∇Nun+1s(t)dt)N. (3.6)

Integrating (3.6) from to , we have

 ∥∥ ∥∥dun+1s(t)dt∥∥ ∥∥2L2(tn,tn+1;ℓ2)+Aτ2∥∥ ∥∥dΔNun+1s(t)dt∥∥ ∥∥2L2(tn,tn+1;ℓ2) +EN(un+1s)−EN(uns) = ∫tn+1tn(β(∇Nuns)−β(∇Nun+1s(t)),d∇Nun+1s(t)dt)N dt +∫tn+1tnt−tnτ(β(∇Nuns)−β(∇Nun−1s),d∇Nun+1s(t)dt)N dt := (I)+(II). (3.7)

As for , an application of Hlder’s inequality gives

 (I) ≤∫tn+1tn∥β(∇Nun+1s(t))−β(∇Nuns)∥N∥∥ ∥∥d∇Nun+1s(t)dt∥∥ ∥∥N dt. (3.8)

Using the inequality in [21, Lemma 3.5] and reversing the order of integration and summation, we obtain

 ∥β(∇Nun+1s(t))−β(∇Nuns)∥N ≤∥∇Nun+1s(t)−∇Nuns∥N =∥∥ ∥∥∫ttnd∇Nun+1s(l)dt dl∥∥ ∥∥N ≤τ12∥∥ ∥∥d∇Nun+1s(t)dt∥∥ ∥∥L2(tn,tn+1;ℓ2). (3.9)

Substituting (3.9) into (3.8) and applying Hlder’s inequality again, we get

 (I) ≤τ∥∥ ∥∥d∇Nun+1s(t)dt∥∥ ∥∥2L2(tn,tn+1;ℓ2). (3.10)

Also notice that by the interpolation inequality, we have

 ∥∥ ∥∥d∇Nun+1s(t)dt∥∥ ∥∥2N≤12τλ∥∥ ∥∥dun+1s(t)dt∥∥ ∥∥2N+τλ2∥∥ ∥∥% dΔNun+1s(t)dt∥∥ ∥∥2N,

where is a positive constant to be decided later. Thus we arrive:

 (I) ≤12λ∥∥ ∥∥dun+1s(t)% dt∥∥ ∥∥2L2(tn,tn+1;ℓ2)+τ2λ2∥∥ ∥∥dΔNun+1s(t)dt∥∥ ∥∥2L2(tn,tn+1;ℓ2). (3.11)

The term could be estimated in a similar way:

 (II) ≤∫tn+1tnt−tnτ∥∇Nuns−∇Nun−1s∥N∥∥ ∥∥d∇Nun+1s(t)dt∥∥ ∥∥N dt ≤τ12∥∥∥d∇Nuns(t)dt∥∥∥L2(tn−1,tn;ℓ2)∫tn+1tnt−tnτ∥∥ ∥∥d∇Nun+1s(t)dt∥∥ ∥∥N dt ≤τ√3∥∥∥d∇Nuns(t)dt∥∥∥L2(tn−1,tn;ℓ2)∥∥ ∥∥d∇Nun+1s(t)dt∥∥ ∥∥L2(tn,tn+1;ℓ2) ≤τ2√3∥∥∥d∇Nuns(t)dt∥∥∥2L2(tn−1,tn;ℓ2)+τ2√3∥∥ ∥∥