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.
