The S(α, β) double-differential thermal neutron scattering law tabulated in Evaluated Nuclear Data File (ENDF) File 7 is, by convention, produced theoretically through fundamental scattering physics models. Currently, no published ENDF evaluations contain covariance data for S(α, β) or associated scattering cross sections. Furthermore, no accepted methodology exists for quantifying or representing these covariances. Thermal scattering cross sections depend on the interatomic structure and dynamics of the material. For many solids, the influence of these properties on inelastic scattering cross sections can be adequately described through the phonon energy spectrum. The phonon spectrum can be viewed as a probability density function and is commonly the fundamental input for calculating S(α, β). Probable variation in the shape of the phonon spectrum may be established that characterizes uncertainties in the physics models and methodology employed in its production. Through Monte Carlo sampling of perturbations from the reference phonon spectrum, an S(α, β) covariance matrix may be generated. With appropriate sensitivity information, the S(α, β) covariance matrix can be propagated to generate covariance data for differential and integral cross sections. In this work, hexagonal graphite is used as an example material for demonstrating the proposed procedures for analyzing, calculating, and processing uncertainty information for theoretically generated thermal neutron inelastic scattering data.