  In this paper, a new two-parameter family of regularized kernels is
introduced, suitable for applying high-order time stepping to N-body systems.
These high-order kernels are derived by truncating a Taylor expansion of the
non-regularized kernel about $(r^2+\epsilon^2)$, generating a sequence of
increasingly more accurate kernels. This paper proves the validity of this
two-parameter family of regularized kernels, constructs error estimates, and
illustrates the benefits of using high-order kernels through numerical
experiments.
