Accurate representation of thermal neutron scattering in Monte Carlo transport simulations requires that the molecular vibrations of the target material be accounted for. Historically, this has been achieved by precomputing large multidimensional tables that are a function of temperature and the cosine of the scattering angle, as well as incoming and outgoing neutron energy. Most commonly used sampling techniques for thermal neutron scattering rely on large multidimensional tables, where higher resolution results in an increase in required memory and attempts to reduce memory can result in grid coarseness errors. An alternative sampling method is introduced here that is a significant departure from precomputed tables and instead relies on a more physical model of the scattering behavior. The phonon sampling method classifies neutron scattering events by the number of phonons excited/de-excited during the scattering collision. In doing so, energy exchange may be obtained via rejection sampling, and an analytical representation of the momentum exchange is obtained. This sampling method has been tested on graphite, yttrium hydride, and uranium nitride, and preliminary implementation of the phonon sampling method shows accurate results for angular and energy distributions, though resulting in up to a 40% slowdown in overall calculation time. This notable slowdown is countered, however, by a large reduction in storage (over 99% reduction compared to standard multidimensional tables).