In several relativistic astrophysics scenarios, the understanding of the rich magnetohydrodynamics is hampered by the limitations set by the achievable numerical resolution. In these cases, it is a tremendous challenge to accurately simulate numerically all the relevant scales. We present how to study such systems by using large eddy simulations with a self-consistent subgrid-scale gradient model that we generalized to the special relativistic case in a previous work and now extend to the general relativistic case. Adapted from nonrelativistic scenarios, the so-called gradient model allows us to capture part of the effects of the hidden dynamics on the resolved scales, by means of a physically agnostic, mathematically based Taylor expansion of the nonlinear terms in the conservative evolution equations' fluxes. One of the main applications is the binary neutron star mergers, where the collision excites a nontrivial amplification at small spatial scales. Motivated by this scenario, we assess the validity of this approach in bounding-box simulations of the magnetic Kelvin-Helmholtz instability. Several resolutions and a broad range of scenarios are considered in order to carefully test the performance of the model under three crucial aspects: (i) highly curved backgrounds, (ii) jumps on the fluid density profiles, and (iii) strong shocks. The results suggest that our extension of the gradient subgrid-scale model to general relativistic magnetohydrodynamics is a promising approach for studying binary neutron star mergers and other relevant astrophysical scenarios.

Performing accurate large eddy simulations in compressible, turbulent magnetohydrodynamics is more challenging than in non-magnetized fluids due to the complex interplay between kinetic, magnetic and internal energy at different scales. Here we extend the sub-grid-scale gradient model, so far used in the momentum and induction equations, to account also for the unresolved scales in the energy evolution equation of a compressible ideal MHD fluid with a generic equation of state. We assess the model by considering box simulations of the turbulence triggered across a shear layer by the Kelvin-Helmholtz instability, testing cases where the small-scale dynamics cannot be fully captured by the resolution considered, such that the efficiency of the simulated dynamo effect depends on the resolution employed. This lack of numerical convergence is actually a currently common issue in several astrophysical problems, where the integral and fastest-growing-instability scales are too far apart to be fully covered numerically. We perform a-priori and a-posteriori tests of the extended gradient model. In the former, we find that, for many different initial conditions and resolutions, the gradient model outperforms other commonly used models in terms of correlation with the residuals coming from the filtering of a high-resolution run. In the second test, we show how a low-resolution run with the gradient model is able to quantitatively reproduce the evolution of the magnetic energy (the integrated value and the spectral distribution) coming from higher-resolution runs. This extension is the first step towards the implementation in relativistic magnetohydrodynamics.

The detection of binary neutron star mergers represents one of the most important and complex astrophysical discoveries of the recent years. One of the unclear aspects of the problem is the turbulent magnetic field amplification, initially triggered by the Kelvin-Helmholtz instability at much smaller scales than any reachable numerical resolution nowadays. Here we present numerical simulations of the first 10 milliseconds of a binary neutron star merger. First, we confirm in detail how the simulated amplification depends on the numerical resolution and is distributed on a broad range of scales, as expected from turbulent magnetohydrodynamics theory. We find that an initial large-scale magnetic field of 10 11 G inside each star is amplified in the remnant to root-mean-square values above 10 16 G within the first 5 milliseconds for our highest-resolution run. Then, we run large eddy simulations, exploring the performance of the subgrid-scale gradient model, already tested successfully in previous turbulent box simulations. We show that the addition of this model is especially important in the induction equation, since it leads to an amplification of the magnetic field comparable to a higher-resolution run, but with a greatly reduced computational cost. In the first 10 milliseconds, there is no clear hint for an ordered, large-scale magnetic field, which should indeed occur in longer timescales through magnetic winding and the magnetorotational instability.