Smooth equations of state for high-accuracy simulations of neutron star binaries

High-accuracy numerical simulations of merging neutron stars play an important role in testing and calibrating the waveform models used by gravitational wave observatories. Obtaining high-accuracy waveforms at a reasonable computational cost, however, remains a significant challenge. One issue is that high-order convergence of the solution requires the use of smooth evolution variables, while many of the equations of state used to model the neutron star matter have discontinuities, typically in the first derivative of the pressure. Spectral formulations of the equation of state have been proposed as a potential solution to this problem. Here, we report on the numerical implementation of spectral equations of state in the spectral Einstein code. We show that, in our code, spectral equations of state allow for high-accuracy simulations at a lower computational cost than commonly used “piecewise polytrope” equations state. We also demonstrate that not all spectral equations of state are equally useful: different choices for the low-density part of the equation of state can significantly impact the cost and accuracy of simulations. As a result, simulations of neutron star mergers present us with a trade-off between the cost of simulations and the physical realism of the chosen equation of state.

Figure 1

Mass-radius relationships for the grid of spectral equations of state with ${\mathrm{\Gamma }}_0=1.35692$ (black) and ${\mathrm{\Gamma }}_0=2$ (red) provided in the appendix, Tables 3 and 4. Both choices of ${\mathrm{\Gamma }}_0$ are compatible with the most likely range of mass and radii, but the lower ${\mathrm{\Gamma }}_0$ choice allows for the construction of stiffer equations of state.

Figure 2

Mass-radius relationships for the equations of state from Table 1, and for the piecewise-polytropic equation SLyPP.

Figure 3

Density profile of a $1.35M_{\odot }$ neutron star for our SLy-like equations of state. Both the density and radius are in $G=c=M_{\odot }=1$ units, and we plot the density against the radius in isotropic coordinates (i.e., for a conformally flat metric). We note that the isotropic radius differs from the Schwarzschild radius quoted elsewhere in this paper and more commonly used in the literature. We use it here as isotropic coordinates are closer to the coordinates used in our simulations.

Figure 4

Error in the orbital phase of an equal mass neutron star-neutron star (NSNS) binary evolution as a function of time, using the $\mathrm{SLy}\mathrm{\Gamma }135$ equation of state and the SLyPP equation of state. Here and in the following figures, abrupt changes in the slope of our error estimates occur when the conservative error estimate switches between using the estimate obtained from comparing the low and high resolution results and the estimate obtained from comparing the medium and high resolution results.

Figure 5

Error in the orbital phase of an equal mass NSNS binary evolution as a function of time, using the $\mathrm{SLy}\mathrm{\Gamma }2$ equation of state and a fast or slow trigger of the spectral adaptive mesh refinement algorithm. The plot shows $\sim 6.5$ orbits of evolution. Merger would occur after $\sim 10$ orbits, at $t\sim 7000$ in these units.

Figure 6

Error in the orbital phase of an equal mass NSNS binary evolution as a function of time, using the $\mathrm{SLy}\mathrm{\Gamma }2$ and $\mathrm{SLy}\mathrm{\Gamma }135$ equations of state.

Figure 7

Number of time steps taken as a function of time for equal mass NSNS binary evolutions using the $\mathrm{SLy}\mathrm{\Gamma }2$ and $\mathrm{SLy}\mathrm{\Gamma }135$ equations of state (high resolution).