Gravitational waves and mass ejecta from binary neutron star mergers: Effect of large eccentricities

As current gravitational wave (GW) detectors increase in sensitivity, and particularly as new instruments are being planned, there is the possibility that ground-based GW detectors will observe GWs from highly eccentric neutron star binaries. We present the first detailed study of highly eccentric BNS systems with full $(3+1)\mathrm{D}$ numerical relativity simulations using consistent initial conditions, i.e., setups which are in agreement with the Einstein equations and with the equations of general relativistic hydrodynamics in equilibrium. Overall, our simulations cover two different equations of state (EOSs), two different spin configurations, and three to four different initial eccentricities for each pairing of EOS and spin. We extract from the simulated waveforms the frequency of the $f$-mode oscillations induced during close encounters before the merger of the two stars. The extracted frequency is in good agreement with $f$-mode oscillations of individual stars for the irrotational cases, which allows an independent measure of the supranuclear equation of state not accessible for binaries on quasicircular orbits. The energy stored in these $f$-mode oscillations can be as large as ${10}^{-3}{M}_{\odot }\sim {10}^{51}\mathrm{erg}$, even with a soft EOS. In order to estimate the stored energy, we also examine the effects of mode mixing due to the stars’ offset from the origin on the $f$-mode contribution to the GW signal. While in general (eccentric) neutron star mergers produce bright electromagnetic counterparts, we find that for the considered cases with fixed initial separation the luminosity decreases when the eccentricity becomes too large, due to a decrease of the ejecta mass. Finally, the use of consistent initial configurations also allows us to produce high-quality waveforms for different eccentricities which can be used as a test bed for waveform model development of highly eccentric binary neutron star systems.

Figure 1

Tracks of the star centers for the ${\mathrm{SLy}}^{00}$ cases. Each panel refers to a different input eccentricity: $e=0.40$ (upper left, red), $e=0.45$ (upper right, green), $e=0.50$ (lower left, blue), $e=0.60$ (lower right, orange). The star that originally starts on the left-hand side always has its track colored black.

Figure 2

The evolution of the proper distance between the two non-spinning NSs with SLy EOS and with different initial eccentricities. For the system with less initial eccentricity, the stars do not approach each other as closely, compared to the system with relatively larger initial eccentricity. In other words, the larger the initial eccentricity the smaller is the impact parameter.

Figure 3

NS trajectories computed from the coordinate positions of the lapse minima for the two stars of the $\mathrm{SLy}-e{0.50}^{00}-\text{R}2$ configuration. We also show 2d-plots of the density and velocity field at different times, with the unbound material [computed using the criterion in Eq. (2)] shown in the brown to green color scale, while the bound material is shown in a blue to red color scale. The times corresponding to the inset plots are marked with colored circles on the trajectory. Top left (cyan): Time after the first encounter, where the stars come as close as $\sim 45\mathrm{km}$ (proper distance). As the stars undergo a tightly bound whirl orbit, one can see that some low-density matter is ejected. Top right (orange): After periastron encounter at a later time in the evolution. Here also we see similar features as in the top left plot except that there is also mass transfer at a density $\sim \mathcal{O}({10}^{-9})–\mathcal{O}({10}^{-8})$ [$\mathcal{O}({10}^8)–\mathcal{O}({10}^9)\mathrm{g}/{\mathrm{cm}}^3$]. Bottom right (yellow): Time just before the merger. Unbound matter with densities $\mathcal{O}({10}^{-7})–\mathcal{O}({10}^{-6})$ [$\mathcal{O}({10}^{10})–\mathcal{O}({10}^{11})\mathrm{g}/{\mathrm{cm}}^3$] is ejected from the tidal tail and as the unbound matter expands its density decreases. Bottom left (yellow): Post-merger phase when the cores of the NSs have merged and formed a hypermassive NS.

Figure 4

Binding energy vs specific angular momentum curves $E_b(\ell )$ for SLy EOS for nonspinning (top panel) and aligned spin configurations (bottom panel). The moments of merger are marked by circles.

Figure 5

Binding energy vs specific angular momentum curves $E_b(\ell )$ for MS1b EOS for nonspinning (top panel) and aligned spin configurations (bottom panel). The moments of merger are marked by circles.

Figure 6

Bolometric luminosity for the SLy setups with zero spin configurations and varying initial eccentricity. The luminosities are computed using the publicly available BNS Kilonova Lightcurve Calculator [128]. The initial times for the bolometric lightcurves differ due to the different ejecta masses; see discussion in [115].

Figure 7

Absolute magnitudes in the ugriz-bands and JHK-bands for the SLy setups with zero spin configurations and varying initial eccentricity. The magnitudes are computed with the BNS Kilonova Lightcurve Calculator available at [128]; see Ref. [115] for details.

Figure 8

Real part of the (2,2) mode of the gravitational wave strain $rh$ vs the retarded time $u$. The color code refers to Fig. 1 and corresponds to different eccentricities. See, e.g., [26] for eccentric BBH waveforms, for comparison (albeit for smaller eccentricities).

Figure 9

$r{\mathrm{\Psi }}_4^{22}$ extracted at $r=900M_{\odot }$ and instantaneous GW frequency $M{\omega }_{22}$ computed from the phase of $r{\mathrm{\Psi }}_4^{22}$ (note that the curves at different extraction radii are coincident) plotted alongside $2\times M\mathrm{\Omega }$ as computed from the coordinate tracks of the NSs, both for the inspiral part of the $\mathrm{SLy}-e{0.50}^{00}-\mathrm{R}4$ case. We find a good agreement for the relation $M{\omega }_{22}=2\times M\mathrm{\Omega }$ and similar results hold for all the other cases. The higher frequency in-between the periastron encounters correspond to the quasinormal mode excitations of the NSs. We also plot the perturbation theory (PT) estimate for the $f$-mode frequency as a black dashed line.

Figure 10

Spectrogram (left) and power spectral density (right) of the curvature scalar $r{\mathrm{\Psi }}_4$ for the $\mathrm{SLy}-e{0.50}^{00}-\mathrm{R}2$ configuration. The labels in the right panel plot refer to the (22), (21), and (20) modes of the curvature scalar. The black dashed line marks the PT estimate of the $f$-mode frequency and the gray line marks the moment of merger. The PSD is split into a part before the moment of merger (black) and a part after the merger (crimson).

Figure 11

Spectrogram (left) and power spectral density (right) of the curvature scalar $r{\mathrm{\Psi }}_4$ for the $\mathrm{SLy}-e{0.50}^{\uparrow \uparrow }-\text{R}2$ configuration. The labels in the right panel plot refer to the (22), (21), and (20) modes of the curvature scalar. The black dashed line marks the PT estimate of the $f$-mode frequency for nonspinning stars and the gray line marks the moment of merger. The PSD is split into a part before the moment of merger (black) and a part after the merger (crimson).

Figure 12

Instantaneous frequency of ${\mathrm{\Psi }}_{2,2}^{4,f\text{−}\mathrm{mode}}$ for $\mathrm{SLy}-e{0.50}^{00}-\text{R}4$, compared to the prediction of the $f$-mode frequency of an isolated star, as well as the $f$-mode frequency including the estimated gravitational redshift due to the star’s companion. These curves should only be compared away from the periastron bursts.

Figure 13

Illustration of the effects of displacement-induced mode mixing on the waveforms extracted from $\mathrm{SLy}-e{0.50}^{00}-\text{R}4$: Coordinate distance from star 1 to the origin (top); mode mixing coefficients (second-from-top) ($|{\mu }_{2,\pm 2;2,\mp 2}|$ is almost indistinguishable from 0 at this scale, so we do not show it); illustrating the averaging used to separate the $f$-mode and orbital signals (third-from-top); comparing the full mode mixing calculation of the intrinsic $f$-mode contribution to the real part of the (2,2) mode from a single star to the naive calculation of taking half of the extracted mode (second-from-bottom) (the result for the imaginary part is exactly analogous); comparing the contribution to the real part of the (2,0) mode from mode mixing to that extracted from the evolution of the binary (bottom) (the imaginary part vanishes). In all cases, the approximations used to compute the mode mixing are only valid during the inter-burst $f$-mode oscillations, and the predictions of the mode mixing calculation for the periastron bursts should not be considered.

Figure 14

Estimate of the energy stored in $f$-mode oscillations as a function of the retarded time. (Left panel): Comparison of the results after accounting for displacement-induced mode mixing to those obtained using half the extracted signal, both plain (i.e., the pure extracted signal), and with only the $f$-mode part for $\mathrm{SLy}-e{0.50}^{00}-\text{R}4$. The energy estimate is only valid in the times between the periastron bursts. (Middle and right panels): The energy estimate computed accounting for displacement-induced mode mixing for the nonspinning SLy and MS1b cases that have multiple encounters before merger, considering the R2 resolutions for which we have data for the most eccentricities. In order to clarify these plots, we have excluded the portions of the traces that are close to the periastron bursts.

Figure 15

Rest mass conservation of $\mathrm{SLy}-e{0.50}^{00}$ case on refinement level $l=1$. The panel shows the absolute error of the rescaled mass evolution of the baryonic mass. We mark the moment of merger as circles and the moment when a BH forms with a star. The error in the total mass stays below 0.6% for the lowest resolution and below 0.01% for the highest resolution over the entire evolution time. Mass conservation is not violated once a BH forms.

Figure 16

ADM constraints for the $\mathrm{SLy}-e{0.50}^{00}$ case. The upper panel shows the $L^2$ volume norm of the Hamiltonian constraint, $||\mathcal{H}|{|}_2$. The lower panel shows the Euclidean norm of the $L^2$ volume norms of the Cartesian components of the momentum constraint, $||\overrightarrow{\mathcal{M}}|{|}_2=\sqrt{||{\mathcal{M}}^x|{|}_2^2+||{\mathcal{M}}^y|{|}_2^2+||{\mathcal{M}}^z|{|}_2^2}$. The constraints are evaluated on refinement level 4 and are decreasing for increasing resolution during the inspiral of the neutron stars. Note that here we restrict our convergence analysis to the refinement level 4, cf. Table 1. This is needed since for the coarser refinement levels even the higher resolution setups are only barely able to resolve the star whereas for the coarser resolutions not a single point covers the NS leading to a nonconvergent behavior. We have filtered our data with an average filter to give a better visualization and reduce high frequency noise and mark the times of merger by circles.

Figure 17

(Top panel): Real part of the curvature multipole $r{\mathrm{\Psi }}_4^{22}$ for the $\mathrm{SLy}-e{0.50}^{00}$ configuration for four different resolution, see Table 2. The inset figure shows two periastron encounters and the region in between where the $f$-mode oscillations are present. (Bottom panel): The instantaneous dimensionless frequency $M{\omega }_{22}$ computed from $r{\mathrm{\Psi }}_4^{22}$. The plots show the robustness of the NS oscillations and the different evolution of the postmerger part. The merger times also vary depending on the resolution due to varying numerical dissipation, which is an inevitable numerical phenomenon for all simulations.