Monte-Carlo closure for moment-based transport schemes in general relativistic radiation hydrodynamics simulations

General relativistic radiation hydrodynamics simulations are necessary to accurately model a number of astrophysical systems involving black holes and neutron stars. Photon transport plays a crucial role in radiatively dominated accretion disks, while neutrino transport is critical to core-collapse supernovae and to the modeling of electromagnetic transients and nucleosynthesis in neutron star mergers. However, evolving the full Boltzmann equations of radiative transport is extremely expensive. Here, we describe the implementation in the general relativistic SpEC code of a cheaper radiation hydrodynamics method which theoretically converges to a solution of Boltzmann's equation in the limit of infinite numerical resources. The algorithm is based on a gray two-moment scheme, in which we evolve the energy density and momentum density of the radiation. Two-moment schemes require a closure which fills in missing information about the energy spectrum and higher-order moments of the radiation. Instead of the approximate analytical closure currently used in core-collapse and merger simulations, we complement the two-moment scheme with a low-accuracy Monte-Carlo evolution. The Monte-Carlo results can provide any or all of the missing information in the evolution of the moments, as desired by the user. As a first test of our methods, we study a set of idealized problems demonstrating that our algorithm performs significantly better than existing analytical closures. We also discuss the current limitations of our method, in particular open questions regarding the stability of the fully coupled scheme.