  We present a systematic numerical relativity study of the impact of different
treatment of microphysics and grid resolution in binary neutron star mergers.
We consider series of simulations at multiple resolutions comparing
hydrodynamics, neutrino leakage scheme, leakage augmented with the M0 scheme
and the more consistent M1 transport scheme. Additionally, we consider the
impact of a sub-grid scheme for turbulent viscosity. We find that viscosity
helps to stabilise the remnant against gravitational collapse but grid
resolution has a larger impact than microphysics on the remnant's stability.
The gravitational wave (GW) energy correlates with the maximum remnant density,
that can be thus inferred from GW observations. M1 simulations shows the
emergence of a neutrino trapped gas that locally decreases the temperature a
few percent when compared to the other simulation series. This
out-of-thermodynamics equilibrium effect does not alter the GW emission at the
typical resolutions considered for mergers. Different microphysics treatments
impact significantly mass, geometry and composition of the remnant's disc and
ejecta. M1 simulations show systematically larger proton fractions. The
different ejecta compositions reflect into the nucleosynthesis yields, that are
robust only if both neutrino emission and absorption are simulated. Synthetic
kilonova light curves calculated by means of spherically-symmetric
radiation-hydrodynamics evolutions up to 15 days post-merger are mostly
sensitive to ejecta's mass and composition; they can be reliably predicted only
including the various ejecta components. We conclude that advanced microphysics
in combination with resolutions higher than current standards appear essential
for robust long-term evolutions and astrophysical predictions.
