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.