Consider the continuum model of no-slope-selection thin film epitaxial growth, which is the gradient flow of the following energy functional:
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):
Taking the inner product with (1.2) by , we obtain the energy decay property for the continuous system:
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.,
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]:
For any and , we have
Now we introduce an interpolation operatoronto that reserves the function value on grid points, i.e., for :
where the coefficients
are given by the discrete Fourier transform of thegrid points:
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]):
for dimension . Also, the following bound of is excerpted in [18, Lemma 1]:
For any in dimension , we have
Given , the discrete spatial partial derivatives can be defined as:
Similarly we have and . For any , and , we can define the discrete gradient, divergence and Laplace operators:
Also, to measure the discrete differentiation operators defined above, we introduce the discrete (denoted as ) inner product and norm :
Similarly, we can define the discrete Sobolev norm and the discrete Sobolev semi-norm :
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]:
For any , and , we have the following summation by parts formula:
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:
Define , which is symmetric positive definite on .
In turn, the spatial discretization of (1.2) becomes: Given , find such that
where . Multiplying both sides of (2.6) by yields
Integrating (2.7) from to gives
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 :
For the ETD1, the term in is simply approximated by . For the ETDMs2, is approximated by the linear Lagrange interpolation:
Let be the numerical solution of ETD1 and ETDMs2. Denote as for . Integrating from to , the following expressions are obtained in :
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 :
We present the following sETDMs2 scheme:
sETDMs2: For , find such that for any ,
with for .
The initial step is obtained as follow: Find such that
sETDMs2 (matrix form): For ,
and the intial step is computed by
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
Consider the following numerical energy functional:
where is the discrete version of the continuous energy functional in (1.1):
The main theoretical result of this section is the following theorem.
Assume that . The numerical system (2.11) is energy stable in the sense that for any ,
We define and denote
By definition we have: For any and ,
Taking an inner product with (2.11) by gives that: For ,
Recall that , thus
Therefore, we can rewrite (3.5) as below:
Integrating (3.6) from to , we have
As for , an application of Hlder’s inequality gives
Using the inequality in [21, Lemma 3.5] and reversing the order of integration and summation, we obtain
Also notice that by the interpolation inequality, we have
where is a positive constant to be decided later. Thus we arrive:
The term could be estimated in a similar way: