Hybrid approach to long-term binary neutron-star simulations

One of the main challenges in the numerical modeling of binary neutron-star mergers are long-term simulations of the postmerger remnant over timescales of the order of seconds. When this modeling includes all the aspects of the complex physics accompanying the remnant, the computational costs can easily become enormous. To address this challenge in part, we have developed a novel hybrid approach in which the solution from a general-relativistic magnetohydrodynamics (GRMHD) code solving the full set of the Einstein equations in Cartesian coordinates is coupled with another GRMHD code in which the Einstein equations are solved under the conformally flat condition (CFC). The latter approximation has a long history and has been shown to provide an accurate description of compact objects in nonvacuum spacetimes. An important aspect of the CFC approximation is that the elliptic equations need to be solved only for a fraction of the steps needed for the underlying hydrodynamical/magnetohydrodynamical evolution, thus allowing for a gain in computational efficiency that can be up to a factor of $\sim 6(230)$ in three-dimensional (two-dimensional) simulations. We present here the basic features of the new code, the strategies necessary to interface it when importing both two- and three-dimensional data, and a novel and robust approach to the recovery of the primitive variables. To validate our new framework, we have carried out a number of tests with various coordinates systems and different numbers of spatial dimensions, involving a variety of astrophysical scenarios, including the evolution of the postmerger remnant of a binary neutron-star merger over a timescale of one second. Overall, our results show that the new code, bhac+, is able to accurately reproduce the evolution of compact objects in nonvacuum spacetimes and that, when compared with the evolution in full general relativity, the CFC approximation reproduces accurately both the gravitational fields and the matter variables at a fraction of the computational costs. This opens the way for the systematic study of the secular matter and electromagnetic emission from binary-merger remnants.

Figure 1

Average relative error ${ϵ}_{\text{recov}}$ (left panel), number of iterations (middle panel), and EOS calls (right panel) required to reach the desired tolerance in recovering the primitive variables from the conserved ones. The example shown employs the LS220 EOS [107] with parameter values $W=2$ and $Y_e=0.1$. The magnetic field is nonzero and set so that the pressure ratio is $b^2/2p={10}^{-3}$, with $b^2≔b_i{b}^i$ being the strength of the magnetic field in the fluid frame.

Figure 2

Number of iterations required to reach the desired tolerance in recovering the primitive variables from the conserved ones when using the LS220 EOS with parameters values $\rho ={10}^{11}\mathrm{g}{\mathrm{cm}}^{-3}$, $T=5\mathrm{MeV}=5.8\times {10}^{10}\mathrm{K}$, and $Y_e=0.1$.

Figure 3

Top panel: Relative difference in the evolution of the central rest-mass density ${\rho }_c(t)$ when normalized to the initial one for two nonrotating stars when evolved in either the Cowling approximation (BU0-cow, red solid line) or with a dynamical spacetime within the CFC approximation (BU0-dyn, blue solid line). Bottom panel: the same as in the top one but for the central value of the lapse function ${\alpha }_c(t)$.

Figure 4

PSD of the evolution of the normalized central rest-mass density ${\rho }_c(t)/{\rho }_c(0)$ of the BU0-dyn test computed over a timescale of 10 ms as in the top panel of Fig. 3. Two clear peaks are visible in the PSD and show a very good match with the expected eigenfrequencies of the fundamental mode ($F$-mode, red dashed line), of its first overtone ($H_1$-mode; blue dashed line) and second overtone ($H_2$-mode; green dashed line) computed from perturbative studies [113].

Figure 5

Evolution of the central rest-mass density normalized to the initial value in the migration test. Reported with solid lines of different color are the evolution by bhac+ (blue line) and by FIL (red line). Note the excellent agreement especially over the first few oscillations; the black dotted line represents the central rest-mass density of the star on the stable branch having the same gravitational mass, which is higher than the asymptotic solutions since it does not account for the matter lost in the nonlinear shocks at the stellar surface.

Figure 6

Left panel: 2D distributions of the rest-mass density (left column) and of the rotational velocity (right column) of a magnetized and differentially rotating star (test magnetized-DRNS in Table 1). The top row refers to the initial time $t=0$, while the bottom row to the final time $t=10\mathrm{ms}$, corresponding to about eight rotational periods, and showing a very good preservation of the axisymmetric equilibrium. Right panel: Indicated with different symbols are the 1D profiles of the rest-mass density, of the angular velocity, and of the toroidal magnetic field, all normalized to their maximum values for the same test in the left panel. The top row displays the profiles at a polar angle $\theta =\pi /4$, while the bottom row shows them on the equatorial plane $\theta =\pi /2$. For all quantities, the dashed lines represent the initial profiles, which are well preserved even after eight rotational periods. The inset in the bottom row reports ${\log }_{10}[\rho ]/{\log }_{10}[{\rho }_c(0)]$ and shows that the surface of the star is captured with a couple of cells only.

Figure 7

Left panel: 1D profiles at $\theta =90°$ and at $t=10\mathrm{ms}$ of the rest-mass density (top row), of the angular velocity (middle row), and of the conformal factor (bottom row) for a rapidly and uniformly rotating neutron star (test DD2RNS-mr in Table 1). Solid lines of different color refer to the evolutions carried out by bhac+ (blue line) or FIL (red line) showing a very good preservation of the equilibrium in the high-density regions of the star. The deviations from the initial profiles at $t=0\mathrm{ms}$ (black dashed lines) are comparable in the two codes and are typical of simulations of rapidly rotating stars (note the logarithmic scale employed here). Right panel: Evolution of the relative difference in the central rest-mass density (top row) and of the conformal factor (bottom row) for the same model reported in the left panel and the same convention for the line types. In addition, each row contains also the data of a high-resolution simulation (DD2RNS-hr, green solid line), highlighting how the differences with FIL can be decreased by a higher spatial resolution.

Figure 8

Evolution of the maximum rest-mass density normalized to the initial value in the 3D head-on collision test. The red line represents the result of FIL, while the blue line shows the evolution in bhac+. Also shown with a black and gray solid line are the evolutions with different efficiency ratios, namely ${\xi }_{\mathrm{eff}}=1/3$ and $1/10$, respectively. The blue solid line is obtained with ${\xi }_{\mathrm{eff}}=1$ and is obviously the closest to the FIL evolution.

Figure 9

2D comparison on the principal planes of the rest-mass density and temperature between in the head-on collision of two neutron stars (test head-on) as computed by bhac+ (left part of each plane) and FIL (right part of each plane). In each column, the top part of each panel reports the rest-mass density, while the bottom part shows the temperature. Finally, the four panels refer to four different times, namely, the instant when the two neutron stars start colliding (top-left panel, $t\simeq 1.89\mathrm{ms}$), that of the maximum compression (top-right panel, $t\simeq 2.26\mathrm{ms}$), that of the minimum compression (bottom-left panel, $t\simeq 2.47\mathrm{ms}$), and that of the late-time evolution (bottom-right panel, $t\simeq 6.47\mathrm{ms}$). Note that the color bars are different in the four panels and the use of negative-$z$ regions is done for visualization purposes only since the simulations actually employ a symmetry across the $z=0$ plane. Finally, note the very good agreement between the two evolutions despite the difference in coordinate systems, truncation order in the solution of the GRMHD equations, and different treatment of the spacetime evolution.

Figure 10

Left panel: 2D slices of the conformal factor (left) and of the rest-mass density (right) for a BNS postmerger remnant at $\overline{t}=50\mathrm{ms}$ as evolved by bhac+ (top part of the panel) and by FIL (bottom part); the data has been handed-off at $t_{\mathrm{HO},1}=20\mathrm{ms}$. Right panel: The same as on the left but at $\overline{t}=100\mathrm{ms}$ and evolved with data handed-off at $t_{\mathrm{HO},2}=50\mathrm{ms}$ (see Fig. 11 for a quantitative comparison in 1D). The white (gray) solid, dashed and dotted contours in the left (right) part of each panel refer to rest-mass densities of ${10}^{13}\mathrm{g}/{\mathrm{cm}}^3$ (solid lines), ${10}^{12}\mathrm{g}/{\mathrm{cm}}^3$ (dashed lines), and ${10}^{11}\mathrm{g}/{\mathrm{cm}}^3$ (dotted lines), respectively. The colormap for the conformal factor is tuned to highlight the location of rest-mass densities of the order of the ${10}^{13}\mathrm{g}/{\mathrm{cm}}^3$, which may be taken as reference for the location of the surface of the HMNS.

Figure 11

1D slices at $z=0$ of the rest-mass density (left column), of the specific internal energy (middle column) and of the conformal factor (right column) for a BNS postmerger remnant at different times, i.e., $\overline{t}=20\mathrm{ms}$ (top row, test DD2BNS-HO@20ms), $\overline{t}=50\mathrm{ms}$ (middle row, test DD2BNS-HO@50ms) and $\overline{t}=100\mathrm{ms}$ (bottom row, test DD2BNS-HO@50ms). While the FIL data is always indicated with a red line, the bhac+ data is shown with different colors depending on the HO time, i.e., with a blue solid line in the top row for $t_{\mathrm{HO},1}$ and with a green solid line in the middle and bottom rows for $t_{\mathrm{HO},2}$. Also reported with an orange solid line in the bottom row is the bhac+ evolution at $t=1000\mathrm{ms}$ (see Fig. 10 for a more qualitative comparison in 2D).

Figure 12

Evolution of the central rest-mass density (top) and of central conformal factor (bottom) as obtained from different evolutions. In particular, the red solid line refers to the FIL simulation carried out till $t=112\mathrm{ms}$, the blue solid line shows the bhac+ evolution with HO at $t_{\mathrm{HO},1}$, while the green solid line shows the bhac+ evolution with HO at $t_{\mathrm{HO},2}$ and continued till $t=1000\mathrm{ms}$. Also shown with a light-blue solid line is the bhac+ evolution with HO at $t_{\mathrm{HO},1}$ but with a rescaling of the conformal factor $\tilde{\psi }=1.00038\psi$ to account for the slightly different spacetimes in FIL and bhac+.

Figure 13

Top panel: Comparison of the evolutions of the central rest-mass density of the BNS remnant when different efficiency ratios are employed. In particular, the top panel reports the evolution with data at $t_{\mathrm{HO},2}$ when ${\xi }_{\mathrm{eff}}=1$ (green solid line; this is the same as in the top panel of Fig. 12), ${\xi }_{\mathrm{eff}}=1/3$ (blue dashed line) and ${\xi }_{\mathrm{eff}}=1/10$ (cyan dotted line). The bottom panel reports instead the relative difference with respect to the ${\xi }_{\mathrm{eff}}=1$ reference evolution (blue and cyan solid lines for ${\xi }_{\mathrm{eff}}=1/3$ and $1/10$, respectively) and shows that the variance is $\lesssim 0.01\%$ even for the most extreme case of ${\xi }_{\mathrm{eff}}=1/10$; the latter simulation was performed with a computational gain of a factor 3.6 with respect to the corresponding ${\xi }_{\mathrm{eff}}=1$ simulation.