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.