General relativistic simulations of collapsing binary neutron star mergers with Monte Carlo neutrino transport

Recent gravitational wave observations of neutron-star-neutron-star and neutron-star-black-hole binaries appear to indicate that massive neutron stars may not be too uncommon in merging systems. These discoveries have led to an increased interest in the simulation of merging compact binaries involving massive stars. In this paper, we present a first set of evolution of massive neutron star binaries using Monte Carlo radiation transport for the evolution of neutrinos. We study a range of systems, from nearly symmetric binaries that collapse to a black hole before forming a disk or ejecting material, to more asymmetric binaries in which tidal disruption of the lower mass star leads to the production of more interesting postmerger remnants. For the latter type of systems, we additionally study the impact of viscosity on the properties of the outflows, and compare our results to two recent simulations of identical binaries performed with the whiskythc code. We find agreement on the black hole properties, disk mass, and mass and velocity of the outflows within expected numerical uncertainties, and some minor but noticeable differences in the evolution of the electron fraction when using a subgrid viscosity model, with viscosity playing a more minor role in our simulations. The method used to account for r-process heating in the determination of the outflow properties appears to have a larger impact on our result than those differences between numerical codes. We also use the simulation with the most ejected material to verify that our newly implemented Lagrangian tracers provide a reasonable sampling of the matter outflows as they leave the computational grid. We note that, given the lack of production of hot outflows in these mergers, the main role of neutrinos in these systems is to set the composition of the postmerger remnant. One of the main potential uses of our simulations is, thus, as improved initial conditions for longer evolutions of such remnants.

Figure 1

Baryon density in the equatorial plane of the binary just before an apparent horizon is detected (top), and 2 ms later (bottom). From left to right, we show simulations SFHo-161-151, SFHo-180-131, SFHo-178-106, and LS220-178-106. Note the different color scales on the top and bottom panels.

Figure 2

Temperature of the fluid in the equatorial plane 4 ms after black hole formation for simulation LS220-178-106. We clearly see the distinction between the hot disk $T\sim (5–20)\mathrm{MeV}$ and the cold extended tidal tail.

Figure 3

Baryon density and electron fraction of the fluid 2 ms after detection of an apparent horizon for simulations SFHo178-106 (top), and LS220-178-106 (bottom). We show the vertical slice $y=0$.

Figure 4

Vertical slice through the post-merger remnant 10 ms after collapse to a black hole. We show the SFHo-178-106 (top) and LS220-178-106 (bottom) cases and the baryon density (left), electron fraction (center), and temperature (right). The structure of the accretion disk is very similar in both simulations, but the LS220 simulation has a higher electron fraction in the polar regions and a disk wind absent at that time in the SFHo simulation.

Figure 5

Distribution of $Y_e$ for the unbound matter of simulations LS220-178-106-vis (solid black) and LS220-178-106-novis (dashed red). In both cases, very neutron rich material with $Y_e\lesssim 0.05$ heavily dominates. We also show the result of a nonviscous simulation using a lower number of MC packets (${10}^7$, dot-dashed green). Insufficient MC resolution during collapse leads to the ejection of a small amount of high-$Y_e$ matter.

Figure 6

Corner plot of the distribution of unbound matter leaving the grid as a function of time after the formation of an apparent horizon (estimate of the time at which the matter crosses a sphere $\sim 180\mathrm{km}$ from the black hole), temperature, polar angle, and velocity. We also show $(1,2,3)-\sigma$ contours of the distribution function in each 2D plot. The data in this plot is from simulation LS220-178-106. Plot generated using the corner python library from [60].

Figure 7

Left: neutrino luminosity at the boundary of our computational domain in the LS220-178-106-vis simulation for ${\nu }_e$ (solid black), ${\overline{\nu }}_e$ (dashed red), and ${\nu }_x$ (all species combined, dot-dashed green). The drop at the merger time is due to an extension of the computational domain. We note that the neutrino luminosity remains low over the entire evolution, by the standard of neutron star mergers. The other simulations have comparable peaks before apparent horizon formation, and even lower luminosity at later times. Right: energy spectrum of the escaping neutrinos, in 1 MeV bins. For neutrinos emitted close to the surface of the star (low-energy ${\nu }_e$), the structure of the bins used in the simulation remains very visible at the outer boundary. Each packet is weighted by the number of neutrinos that it represents (not by energy) and we show only the neutrinos leaving the grid after collapse of the remnant (including all neutrinos does not meaningfully change these results). We additionally plot for each species a Fermi-Dirac distribution with the same average energy as the observed neutrinos ($[12.3,18.1,23.8]\mathrm{MeV}$).

Figure 8

Distribution of neutrinos as a function of $\cos \theta$, with $\theta$ the angle with respect to the axis of rotation of the remnant.

Figure 9

Density (top) and temperature (bottom) of the matter leaving the computational grid according to Lagrangian tracers following the fluid from the beginning of the inspiral (blue) or from observing directly matter fluxes on the grid boundary. Despite the long evolution of the tracers, the two remain in decent qualitative agreement, with the tracers slightly oversampling higher density regions.

Figure 10

Properties of the fluid according to the Lagrangian tracers. Here, $\theta$ is the angle with respect to the rotation axis of the remnant. We clearly see the increase in density of the core before collapse to a black hole (at $t=0$ in these coordinates), followed by the formation of an accretion disk at radius $(20–30)\mathrm{km}$ at density $\sim {10}^{12}g/{\mathrm{cm}}^3$. As before, we show $(1,2,3)-\sigma$ contours of the various distributions, weighted by mass. We note that tracers that cross the apparent horizon of the black hole are removed from the simulation.