  We consider a continuum phase field model for crystal growth via molecular
beam epitaxy, with the goal of determining stable numerical time integration
methods for the dynamics. We parametrize a class of semi-implicit methods that
are linear in the updated field, which allows for efficient implementation with
fast Fourier transforms. We perform unconditional von Neumann stability
analysis to identify the region of stability in parameter space, and then test
these predictions numerically for gradient stability. We find strong agreement
between the approaches.
