GR-Athena++: Binary Neutron Star Merger Simulations with Neutrino Transport

TL;DR

This work presents GR-Athena++ with a grey M1+N0 neutrino transport coupled to GRMHD on adaptive mesh refinement, validated through a comprehensive test suite and cross-code comparisons. It demonstrates a split-step numerical strategy, a robust semi-implicit treatment of stiff radiation-matter coupling, and a novel tapering excision technique that preserves stability through horizon formation. The authors apply the framework to rotating neutron-star collapse and to equal-mass binary neutron star mergers with two different equations of state, highlighting the impact of solvers (LLF vs HLLE) and magnetic fields on remnant structure, ejecta, neutrino emission, and gravitational-wave signals. While the grey M1+N0 approach captures essential microphysics and remains computationally efficient, the paper identifies limitations relative to spectral transport and outlines a path toward multigroup extensions and improved weak rates to refine ejecta composition and kilonova predictions.

Abstract

We present general-relativistic radiation magnetohydrodynamics simulations of binary neutron star mergers performed with GR-Athena++. Neutrino transport is treated using a moment-based, energy-integrated scheme (M1), augmented by neutrino number density evolution (N0). Our implementation is validated through an extensive suite of standard tests and demonstrated to perform robustly under adaptive mesh refinement. As a first application, we simulate the gravitational collapse of a uniformly rotating, magnetized neutron star, demonstrating stable radiation evolution through apparent-horizon formation using a novel excision technique based on the tapering of state vector evolution inside the horizon. To further test robustness in highly dynamic environments, we apply our code to two demanding binary neutron star merger scenarios. We investigate a long-lived remnant with the DD2 equation of state, evolved with full general-relativistic magnetohydrodynamics and M1 neutrino transport. Following this, a gravitational collapse scenario with the SFHo equation of state is explored. We showcase long-term stable evolution on neutrino cooling time-scales, demonstrating robust handling of excision and stable evolution of the post-collapse accretion phase in three-dimensional mergers with magnetic fields and neutrino radiation.

Figures (27)

• Figure 1: Optically thin advection test: Euler-frame energy density $\widetilde{E}$ at time $t=2$, for a variety of grid resolutions. Note the absence of any profile artifacting from the the fluid velocity interface (at $x=0$). The Mesh sampling is indicated in the inset. For each run the MeshBlock sampling (along the $x$-direction) is $N_B=10$, with boundaries demarcated in gray for the $N_M=400$ run.
• Figure 2: Diffusion test: Euler-frame energy density $\widetilde{E}$ at time $t=10$, for a variety of grid resolutions. Inset indicates differing number of Mesh samples for distinct simulations. We have $x\in[-4,\,4]$ during calculation; for this plot, however, we suppress (symmetric) data for $x<0$ to better illustrate differences in the runs. For each run, the MeshBlock sampling is $N_B=10$, with boundaries demarcated in gray for the $N_M=400$ run.
• Figure 3: Diffusion test, moving-medium: Euler-frame energy density $\widetilde{E}$ at time $t=4$, at a variety of grid resolutions. (Upper) Observe clear convergent trend as number of Mesh samples is increased. A run at $N_M=200$ featuring two additional levels of adaptive refinement, is also shown. The AMR criterion is dynamical, and traces the field local maximum. Absence of refinement boundary artifacts indicates level-to-level transfer of data, and refluxing are robust. (Lower) Selected runs in $\log$-scale comparing (de)-activation of FC of the M1-sector in (solid) dashed lines respectively. Observe absence of trailing-edge artifacting with FC. (Common) MeshBlock sampling is $N_B=10$ with boundaries demarcated in gray for the AMR run in the upper panel, and the $N_M=400$ run in the lower panel.
• Figure 4: Shadow casting: Euler-frame energy density $\widetilde{E}$ at time $t=10$. An incident beam is injected at the left boundary, and impinges upon an absorbing cylinder $\mathcal{C}$. Observe that the field remains regular as it propagates towards the two levels of refinement that contain $\mathcal{C}$, and exits towards the right. In gray we indicate MeshBlock boundaries, where the interior of each has $N_B=(20,\,20,\,4)$ samples along each axis. Observe the static refinement structure of the grid towards the central feature.
• Figure 5: Radiating sphere: on-axis profile of Euler-frame energy density $\widetilde{E}$ at time $t=10$. Two choices for the opacities are selected through $\mathcal{K}$. The numerical M1 profiles compare favourably with the full analytical solution based on the Boltzmann equation. Vertical dashed black line indicates location of the spherical surface. In gray MeshBlock boundaries are depicted. Observe the static refinement structure of the grid towards the central feature.