Large eddy simulations of magnetized mergers of neutron stars with neutrinos

Neutron star mergers are very violent events involving extreme physical processes: dynamic, strong-field gravity; large magnetic field; very hot, dense matter; and the copious production of neutrinos. Accurate modeling of such a system and its associated multimessenger signals, such as gravitational waves, short gamma ray bursts, and kilonova, requires the inclusion of all these processes and is increasingly important in light of advancements in multimessenger astronomy generally, and in gravitational wave astronomy in particular (such as the development of third-generation detectors). Several general relativistic codes have been incorporating some of these elements with different levels of realism. Here, we extend our code mhduet, which can perform large eddy simulations of magnetohydrodynamics to help capture the magnetic field amplification during the merger, and to allow for realistic equations of state and neutrino cooling via a leakage scheme. We perform several tests involving isolated and binary neutron stars demonstrating the accuracy of the code.

Figure 1

Tests of the eikonal equation. By adopting an explicitly spherical (top) or ellipsoidal (bottom) source, we compare the numerically obtained solution with an analytic solution on the $z=0$ plane. Shown in color map is the analytic function $u(x,y,0)$ of Eq. (62) while the solid contours represent the numerical solution obtained from solving the flat eikonal equation, Eq. (52). The numerical solution agrees very well with the analytic solution and maintains the same symmetry. For comparison, we also include the contours (dashed) obtained with the scheme implemented in had from Paper I, which is largely in agreement despite some irregularities along the diagonals.

Figure 2

Perturbed, cold star with the LS220 EoS. The top panel shows the variations in central density, ${\rho }_0(t)/{\rho }_0(0)$, and in central magnetic field magnitude, $|B_0(t)|/|B_0(0)|$. The bottom panel shows the (normalized) power-spectral density of the quantities in the top panel. The domain of this evolution spans ${[-150\mathrm{km},150\mathrm{km}]}^3$ with finest resolution $\mathrm{\Delta }{x}_{\min }=144\mathrm{m}$. The reference frequencies noted in Table 1 of Paper 1 are shown with vertical, dashed, gray lines. Comparing to Fig. 4 of Paper I, the peak frequencies agree quite well.

Figure 3

Hot, rapidly rotating star. A $2.1M_{\odot }$ (baryonic) star spinning at 730 Hz with an initial temperature of 12 MeV. The maximum density and maximum temperature are shown for all evolutions in the top two panels. The total neutrino luminosity for all species and central magnetic field strength are shown for the evolutions using leakage and with a magnetic field, respectively. With the leakage active, the star cools faster, as expected. An initial magnetization, even very large, has only a very small effect, also as expected. Snapshots of this star at $t=5.26\mathrm{ms}$ are shown in Fig. 4.

Figure 4

Hot, rapidly rotating star at late time ($t=5.26\mathrm{ms}$). From left to right are shown: top: Optical depths for ${\nu }_e$, ${\overline{\nu }}_e$, and ${\nu }_x$; middle: The neutrino luminosities $Q_e$, $Q_a$, and $Q_x$; and bottom: The magnetic field components, $B_x$, $B_y$, and $B_z$. all along the equatorial plane. The time evolution for this magnetized star with leakage is shown in (red, dotted line) Fig. 3. The contour lines display constant density surfaces at $\log (\rho )=(10.5,11,12,13)\mathrm{g}/{\mathrm{cm}}^3$.

Figure 5

Binary neutron star with HS EoS. Snapshots of the density at various times $t=(0,1.92,3.84,5.76,7.68,9.6,11.52,13.44,18.24)\mathrm{ms}$ during the coalescence. Notice that the variations both due to the neutrino dynamics and by the magnetic field occur only from the merger onward. The contours display constant density surfaces at $\log \rho =(10.5,11,12,13)\mathrm{g}/{\mathrm{cm}}^3$. The stars first make contact around the time $t=9.5\mathrm{ms}$ (i.e., close to the middle-right panel), and the remnant fluid has largely circularized by the latest time shown (almost 9 ms after merger).

Figure 6

Binary neutron star with HS EoS. Snapshots of the magnetic field strength, and the same constant density iso-surfaces as in Fig. 5, after the merger at times $t=(11.52,13.44,18.24)\mathrm{ms}$. The top row corresponds to the standard simulation while the bottom row shows the one with LES. Both simulations incorporate the leakage scheme.

Figure 7

Binary neutron star with HS EoS. Average magnetic field strength as a function of time, starting approximately at the merger, for the standard simulation and the LES. Clearly, the magnetic field grows faster and reaches higher values with the LES, even though these simulations employ only medium resolution (see, for instance, Fig. 5 in Ref. [24] to see the effect of the resolution on LES).

Figure 8

Binary neutron star with HS EoS. Snapshots of the temperature (left), electron fraction (middle), and neutrino emission rates (right) at the final time of the simulation $t=18.24\mathrm{ms}$, approximately 9 ms after the merger. The top row corresponds to the standard simulation, while the bottom row shows the LES case. Both of them include magnetic field, although with LES it is much stronger. Notice that the main difference is a small dephasing between these two simulations, possibly due to the stronger magnetic field.

Figure 9

Binary neutron star with HS EoS. Luminosities of the different neutrino species as functions of time, starting approximately at the merger for both the standard simulation and the LES. Again, both of them are for stars with magnetic fields. The much stronger magnetic field (roughly 2 orders of magnitude larger) of the LES simulation arising from its amplification during the turbulent phase of the merger produces only small deviations in the neutrino dynamics compared to the standard, magnetized simulation.

Figure 10

Binary neutron star with HS EoS. Main mode of the gravitational waveform as a function of the retarded time (i.e., subtracting the traveling time of the wave to the surface where it is computed), starting approximately at the merger, for the standard simulation and the LES. Again, no significant differences are observed due to the presence of strong magnetic fields.

Figure 11

Convergence test of the binary GW signal. The primary mode $C_{2,2}$ for three different resolutions of the binary evolution (top). The medium and low resolutions differ from the high resolution by factors of 1.25 and ${(1.25)}^2$. The differences in the phase of the signals (middle) and the absolute differences in the signals (bottom), both measures of the error, are shown, as is the rescaled difference expected between the higher two resolutions if the code converges to third order. The phase appears to converge better than third order while the simple differences in $C_{2,2}$ appear convergent at third order. The medium resolution shown here is the run whose results are presented in the previous figures.