  Classical molecular dynamics simulation is performed mostly using the
established velocity Verlet integrator or other symplectic propagation schemes.
In this work, an alternative formulation of numerical propagators for classical
molecular dynamics is introduced based on an expansion of the time evolution
operator in series of Chebyshev and Newton polynomials. The suggested
propagators have, in principle, arbitrary order of accuracy which can be
controlled by the choice of expansion order after that the series is truncated.
However, the expansion converges only after a minimum number of terms is
included in the expansion and this number increases linearly with the time step
size. Measurements of the energy drift demonstrate the acceptable long-time
stability of the polynomial propagators. It is shown that a system of
interacting Lennard-Jones particles is tractable by the proposed technique and
that the scaling with the expansion order is only polynomial while the scaling
with the number of particles is the same as with the conventional velocity
Verlet. The proposed method is, in principle, extendable for further
interaction force fields and for integration with a thermostat, and can be
parallelized to speed up the computation of every time step.
