  In this article, I provide significant mathematical evidence in support of
the existence of short-time approximations of any polynomial order for the
computation of density matrices of physical systems described by arbitrarily
smooth and bounded from below potentials. While for Theorem 2, which is
``experimental'', I only provide a ``physicist's'' proof, I believe the present
development is mathematically sound. As a verification, I explicitly construct
two short-time approximations to the density matrix having convergence orders 3
and 4, respectively. Furthermore, in the Appendix, I derive the convergence
constant for the trapezoidal Trotter path integral technique. The convergence
orders and constants are then verified by numerical simulations. While the two
short-time approximations constructed are of sure interest to physicists and
chemists involved in Monte Carlo path integral simulations, the present article
is also aimed at the mathematical community, who might find the results
interesting and worth exploring. I conclude the paper by discussing the
implications of the present findings with respect to the solvability of the
dynamical sign problem appearing in real-time Feynman path integral
simulations.
