In　this　paper　we　propose　and　analyze　a　(temporally)　third　order　accurate　backward　differentiation　formula　(BDF)　numerical　scheme　for　the　no-slope-selection　(NSS)　equation　of　the　epitaxial　thin　film　growth　model,　with　Fourier　pseudo-spectral　discretization　in　space.　The　surface　diffusion　term　is　treated　implicitly,　while　the　nonlinear　chemical　potential　is　approximated　by　a　third　order　explicit　extrapolation　formula　for　the　sake　of　solvability.　In　addition,　a　third　order　accurate　Douglas-Dupont　regularization　term,　in　the　form　of　-A　Delta　t(2)Delta　N-2(un+1-un),　is　added　in　the　numerical　scheme.　A　careful　energy　stability　estimate,　combined　with　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　l(infinity)(0,　T;l(2))　boolean　AND　l(2)(0,T;Hh(2))　norm,　with　the　help　of　a　linearized　estimate　for　the　nonlinear　error　terms.　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　epsilon　=　0.02　(up　to　T　=　3x10(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.