  The well-known backward difference formulas (BDF) of the third, the fourth
and the fifth orders are investigated for time integration of the phase field
crystal model. By building up novel discrete gradient structures of the
BDF-$\rmk$ ($\rmk=3,4,5$) formulas, we establish the energy dissipation laws at
the discrete levels and then obtain the priori solution estimates for the
associated numerical schemes (however, we can not build any discrete energy
dissipation law for the corresponding BDF-6 scheme because the BDF-6 formula
itself does not have any discrete gradient structures). With the help of the
discrete orthogonal convolution kernels and Young-type convolution
inequalities, some concise $L^2$ norm error estimates (with respect to the
starting data in the $L^2$ norm) are established via the discrete energy
technique. To the best of our knowledge, this is the first time such type $L^2$
norm error estimates of non-A-stable BDF schemes are obtained for nonlinear
parabolic equations. Numerical examples are presented to verify and support the
theoretical analysis.
