  Magnetic fields are expected to play a key role in the dynamics and the
ejection mechanisms that accompany the merger of two neutron stars. General
relativistic magnetohydrodynamic (MHD) simulations offer a unique opportunity
to unravel the details of the ongoing physical processes. Nevertheless, current
numerical studies are severely limited by the fact that any affordable
resolution remains insufficient to fully capture the small-scale dynamo,
initially triggered by the Kelvin-Helmholtz instability, and later sourced by
several MHD processes involving differential rotation. Here, we alleviate this
limitation by using explicit large-eddy simulations, a technique where the
unresolved dynamics occurring at the sub-grid scales (SGS) is modeled by extra
terms, which are functions of the resolved fields and their derivatives. The
combination of high-order numerical schemes, high resolutions, and the gradient
SGS model allow us to capture the small-scale dynamos produced during the
binary neutron star mergers. Here we follow the first 50 milliseconds after the
merger and, for the first time, we find numerical convergence on the magnetic
field amplification, in terms of integrated energy and spectral distribution
over spatial scales. We also find that the average intensity of the magnetic
field in the remnant saturates at $\sim 10^{16}$~G around $5$~ms after the
merger. After $20-30$~ms, both toroidal and poloidal magnetic field components
grow continuously, fed by the winding mechanism that provides a slow inverse
cascade. We find no clear hints for magneto-rotational instabilities, and no
significant impact of the magnetic field on the redistribution of angular
momentum in the remnant in our simulations, probably due to the very turbulent
and dynamical topology of the magnetic field at all stages, with small-scale
components largely dominating over the large-scale ones.
