  In this paper we propose and analyze a (temporally) third order accurate
exponential time differencing (ETD) numerical scheme for the no-slope-selection
(NSS) equation of the epitaxial thin film growth model, with Fourier
pseudo-spectral discretization in space. A linear splitting is applied to the
physical model, and an ETD-based multistep approximation is used for time
integration of the corresponding equation. In addition, a third order accurate
Douglas-Dupont regularization term, in the form of $-A \dt^2 \phi_0 (L_N)
\Delta_N^2 ( u^{n+1} - u^n)$, is added in the numerical scheme. A careful
Fourier eigenvalue analysis results in the energy stability in a modified
version, and a theoretical justification of the coefficient $A$ becomes
available. As a result of this energy stability analysis, a uniform in time
bound of the numerical energy is obtained. And also, the optimal rate
convergence analysis and error estimate are derived in details, in the
$\ell^\infty (0,T; H_h^1) \cap \ell^2 (0,T; H_h^3)$ norm, with the help of a
careful eigenvalue bound estimate, combined with the nonlinear analysis for the
NSS model. This convergence estimate is the first such result for a third order
accurate scheme for a gradient flow. Some numerical simulation results are
presented to demonstrate the efficiency of the numerical scheme and the third
order convergence. The long time simulation results for $\varepsilon=0.02$ (up
to $T=3 \times 10^5$) have indicated a logarithm law for the energy decay, as
well as the power laws for growth of the surface roughness and the mound width.
In particular, the power index for the surface roughness and the mound width
growth, created by the third order numerical scheme, is more accurate than
those produced by certain second order energy stable schemes in the existing
literature.
