Evaluating radiation transport errors in merger simulations using a Monte-Carlo algorithm

Neutrino-matter interactions play an important role in the post-merger evolution of neutron star-neutron star and black hole-neutron star mergers. Most notably, they determine the properties of the bright optical/infrared transients observable after a merger. Unfortunately, Boltzmann's equations of radiation transport remain too costly to be evolved directly in merger simulations. Simulations rely instead on approximate transport algorithms with unquantified modeling errors. In this paper, we use for the first time a time-dependent general relativistic Monte-Carlo (MC) algorithm to solve Boltzmann's equations and estimate important properties of the neutrino distribution function ~10ms after a neutron star merger. We do not fully couple the MC algorithm to the fluid evolution, but use a short evolution of the merger remnant to critically assess errors in our approximate gray two-moment transport scheme. We demonstrate that the analytical closure used by the moment scheme is highly inaccurate in the polar regions, but performs well elsewhere. While the average energy of polar neutrinos is reasonably well captured by the two-moment scheme, estimates for the neutrino energy become less accurate at lower latitudes. The two-moment formalism also overestimates the density of neutrinos in the polar regions by ~50%, and underestimates the neutrino pair-annihilation rate at the poles by factors of 2-3. Although the latter is significantly more accurate than one might have expected before this study, our results indicate that predictions for the properties of polar outflows and for the creation of a baryon-free region at the poles are likely to be affected by errors in the two-moment scheme, thus limiting our ability to reliably model kilonovae and gamma-ray bursts.