Binary neutron star mergers with tabulated equations of state in SPHINCS_BSSN

Abstract

The dynamics and observable signatures of neutron star mergers are governed by physics under the most extreme conditions. They are particularly impacted by the high-density equation of state, which for the most sophisticated models is usually available in the form of tables. Numerical relativity codes usually evolve particularly well-behaved numerical ("conservative") variables, but at the price that the physically interesting ("primitive") variables need to be found at every computational element and at every integration sub-step by means of expensive (and not always successful) root-finding algorithms. We have recently developed the Lagrangian numerical relativity code SPHINCS_BSSN which evolves the spacetime on an adaptive mesh with well tested methods, but the fluid is evolved by means of freely moving particles. Since our evolution equations differ from those of conventional numerical relativity, we need to develop new conservative-to-primitive algorithms if we want to use tabulated equations of state. We present here three such algorithms: a 3D and a 2D Newton-Raphson method and a 1D root-finding algorithm based on Ridders' method. We find the 3D method to be very fast and robust with an average failure fraction in a full-blown neutron star merger simulation (with the DD2 equation of state) well below 1%. While we do not find obvious advantages for the 2D method, the 1D Ridders' method is slow, but essentially fail-safe. Therefore, we choose the 3D Newton-Raphson as default and fall back to the 1D Ridders' method as a safe "parachute".

Figures (5)

• Figure 1: Summary of the recovery methods
• Figure 2: Accuracy of con2prim in the $(\rho, T)$ plane for different schemes for the DD2 EOS. The regions colored red are those where con2prim failed i.e. con2prim did not converge or converged with a recovery error less than $10^{-10}$. The regions colored blue are those where con2prim converged with an accuracy of at least $10^{-10}$, and the shade of blue represents the recovery error with the darker shades representing higher recovery errors. The top panels show the results for different schemes for Case A : flat spacetime with velocity $v=(0.9, 0, 0)$ and generalized Lorentz factor $\Theta \simeq 2.29$. The middle panels show the results for Case B : Cartesian Kerr-Schild spacetime with $v=(0.2, 0, 0)$ and $\Theta \simeq 3.47$. The bottom panels show the results for Case C : Cartesian Kerr-Schild spacetime with $v=(0.9, 0, 0)$ and $\Theta \simeq 10.26$. We find that 3D Newton-Raphson and 2D Newton-Raphson show non-uniform convergence behavior where the recovery error is higher ($\lesssim 10^{-10}$) for densities $\gtrsim 10^{14} \, \mathrm{g/cm}^3$ and lower ($\lesssim 10^{-12}$) elsewhere. On the other hand, the recovery error is uniformly low everywhere for 1D Ridders' method with values $\lesssim 10^{-15}$ at $\Theta=2.29$ and $\lesssim 10^{-13}$ at $\Theta=10.26$. Thus, 1D Ridders' method is more accurate, especially in high-density regions.
• Figure 3: Computational costs associated with different con2prim schemes in terms of number of EOS interpolations required. We plot the histograms of frequency vs number of EOS calls using 20 bins between 0 to 230 (280, 4000) for 3D Newton-Raphson (2D Newton-Raphson, 1D Ridders'). We also show the mean and the median number of EOS calls using dashed red and green lines respectively on top of the histograms. We compare the costs of 3D Newton-Raphson, 2D Newton-Raphson and 1D Ridders' method for $\Theta = 2.29$ in the top panels, $\Theta = 3.47$ in the middle panels and $\Theta = 10.26$ in the bottom panels. We find that 1D Ridders' method, even though the most accurate and robust, is also the most expensive among all schemes, with median EOS calls $\gtrsim 650$. On the other hand, 3D Newton-Raphson and 2D Newton-Raphson, though less robust, are much cheaper with median EOS calls $\lesssim 15$.
• Figure 4: 2D cross-sections of density and temperature of the binary neutron star merger simulation just before the merger at $7.06 \, \mathrm{s}$, just after the merger at $7.86 \, \mathrm{s}$ and the final simulation time at $18.96 \, \mathrm{s}$. The cross-sections are in the $xy$-plane at $z=0$ for the simulation with 3 million particles. We show the density in the top panels and the temperature (capped at 20 MeV for visibility) in the bottom panels.
• Figure 5: Failure fraction of our primary 3D Newton-Raphson based con2prim method as a function of time. At every sub-step, we define the failure fraction as the ratio of number of particles where 3D Newton-Raphson fails and backup 1D Ridders' method is needed to the total number of particles. The blue line shows the average failure fraction at every timestep. The red line shows the average failure fraction over every $0.1 \, \mathrm{ms}$ of time. We find that the average failure rate of 3D Newton-Raphson is $\sim0.1\%$ in the pre-merger phase and $\lesssim0.01\%$ in the post-merger phase.