Recoil Velocity of Binary Neutron Star Merger Remnants

Abstract The LIGO-Virgo gravitational wave detectors have confidently observed 4 events involving neutron stars: two binary neutron star (BNS) mergers (GW170817 and GW190425), and two neutron star-black hole mergers (GW200105 and GW200115). However, our theoretical understanding of the remnant properties of such systems is incomplete due to the complexities related to the modeling of matter effects and the very high computational cost of corresponding numerical relativity simulations. An important such property is the recoil velocity, which is imparted onto the remnant due to the anisotropic emission of gravitational radiation and the dynamical ejection of matter in the kilonova. In this work, we combine gravitational radiation as well as dynamical ejecta distributions, computed by the Computational Relativity numerical simulations, to get accurate estimates for BNS remnant recoil velocities. We find that recoils due to ejection of matter dominate those caused by gravitational wave emission. Knowledge of BNS remnant recoil velocities is important in determining if the remnant is retained by its environment for future hierarchical mergers which, in turn, can form binaries with black holes in the so-called lower mass gap of ∼ 3 – 5M ⊙.


Introduction
Binaries of compact objects containing neutron stars are among the most intriguing events known in the Universe, and understanding these mergers provide insights into astrophysical phenomena such as gamma-ray bursts, the origin of heavy elements in the local Universe, and formation channels of compact object binaries.
Knowing the properties of binary neutron star (BNS) merger remnants has applications in constraining the astrophysical environments in which these mergers take place and their remnants participate in hierarchical mergers. The remnant properties of binary black hole (BBH) mergers are very well studied, thanks to many numerical relativity (NR) simulations performed over the years for BBHs. Consequently, we have reliable and accurate predictions for properties of the BBH merger remnants such as their final mass, spin and recoil velocity [Berti (2007), Buonanno (2008), Gerosa (2016)]. On the other hand, one has only a limited understanding of the remnants of BNS mergers due to the complexities related to the modeling of matter effects and very high computational cost of corresponding NR simulations. However, there are some predictions for the final mass and spin of BNS merger remnants [Coughlin (2019)] but none for the remnant's recoil velocity, mainly because of the additional uncertainty in the computation of the ejecta properties [Endrizzi (2016)]. Therefore, we are estimating recoil velocities of BNS merger remnants using Computational Relativity (CoRe) database of numerical simulations [Deitrich (2018)] and developing a fit for them covering an extended BNS parameter space.

Overview
The loss of energy and angular momentum through gravitational waves (GWs) causes the orbit of a compact binary to shrink. In addition, mass and spin asymmetries in the binary system also lead to the anisotropic emission of linear momentum through GW radiation. This causes the center-of-mass of the system to recoil in the opposite direction. As a result, the remnant object at the end of the merger is imparted with a recoil velocity, or a kick. These kicks have been shown to be particularly large (∼ 10 2 − 10 3 km/s) for BBHs with unequal masses and high spins [see, e.g., Gerosa (2018)].

Kicks due to gravitational wave emission
Compared to BBH systems, BNS systems have low spins and comparable masses. These characteristics are expected to impart low GW recoils due to anisotropic linear momentum emission. Given a strain waveform model h(t), the linear momentum radiated throughout the binary evolution can be calculated using its time derivative.
The emission of momentum can be obtained by summing over the different (l, m) modes ofḣ(t), the time derivative of the GW strain expanded in the basis of spin-weighted spherical harmonics, as given in Ruiz (2008), Gerosa (2018) where P + = P x + iP y , and r is the extraction radius for the NR simulation. The total radiated momentum can be obtained by integrating Eqs. (2.1) and (2.2) over the waveform duration. The final radiated energy is subtracted from the binary's total mass, to obtain the mass of the remnant. The recoil velocity is then the final linear momentum divided by the remnant mass.

Kicks due to ejecta emission
A BNS merger is also accompanied by the asymmetric ejection of matter. The amount of matter ejected depends upon binary characteristics such as the mass ratio (q) and the mass weighted tidal deformability (Λ) of the system as well as the equation of state [Nedora (2021)]. This asymmetric ejecta also imparts a kick onto the central remnant, which we can estimate by doing a simple calculation using the conservation of linear momentum. The net momentum of the ejecta is given by integrating the flux density over the ejecta mass distribution: where M NS is the mass of the remnant (with the net ejecta mass and radiated GW energy subtracted), and v rec is the kick imparted due to ejecta emission.

Results
Figure 1 displays recoil estimates for a set of 50 NR simulations of BNS systems, plotted against the binary mass ratio q, with different equations of states. The mean GW emission recoil is 15.67 km/s, and that for the ejecta emission is 28.9 km/s. This implies that the recoils due to ejecta emission dominate the GW recoils, and the net recoil velocity (assuming same direction for both kicks) can be as high as ∼ 100 km/s. This is important for two reasons: (i) It tells us that we cannot rely upon earlier recoil estimates obtained using only GW emission and BBH waveforms, and (ii) it emphasizes the need to have richer dynamical post-merger data from NR simulations.
The above results come with important caveats. These results do not yet include systematic errors that can be estimated using the simulations at different numerical resolutions. In addition, we hope to study the variation in kicks of BNS systems due to choice of the NS equation of state more systematically.

Recoil velocity fit using Gaussian Process Regression
Having obtained NR-based estimates of BNS remnant recoil velocities for systems from certain areas of the parameters space, the next step is to generate a fit for recoil as a function of BNS parameters such as the mass ratio (q), and the mass-weighted tidal deformability (Λ). To achieve this, we are applying Gaussian Process Regression (GPR), a machine-learning based Bayesian, non-parametric regression tool. GPR provides a probability distribution over possible functions that fit a set of data points. In our case, the data are NR kick estimates for BNS parameters used in the NR simulations. The mean calculated by the posterior distribution of the possible functions is the function used for regression. Figure 2 shows an example 1-parameter fit for recoils across BNS mass ratios using GPR.

Astrophysical Implications
The two main channels for the formation of a compact binary system are the isolated channel wherein components of a massive main-sequence stellar binary individually collapse into compact objects, and the dynamical channel where compact objects present in a dense stellar environment such as a Globular Cluster, Nuclear Star Cluster or an Active Galactic Nucleus (AGN) disk can interact dynamically to come within mutual gravitational potential and form a binary. The latter scenario presents the exciting possibility of hierarchical binary mergers, where at least one of its two components is itself the remnant of a past binary merger. These successive generations of compact objects have interesting distinguishing properties such as higher masses and higher spins. Remnants of BNS mergers can occupy the so-called lower mass gap of compact objects between 3 − 5M that divides masses of neutron stars and black holes. A crucial requirement of hierarchical mergers is the retention of merger remnants in their environment to be available for successive mergers.
For an environment with a given escape velocity (v esc ), the remnant recoil velocity (v rec ) determines whether it would be retained (v rec < v esc ) or whether it would escape the environment (v rec > v esc ). Remnants which escape are not available for future mergers in the hierarchical scheme. For the three stellar cluster environments discussed above, typical values of v esc are: (i) Open Clusters < 10 km/s, (ii) Globular Clusters ∼ 50 km/s, and (iii) Nuclear Star Clusters ∼ 200 km/s Mapelli (2020). We find that recoil speeds of the order of 100 km/s are possible due to BNS ejecta emission, and hence, such remnants can be retained in environments like nuclear star clusters and AGN disks. Some population models show that such remnants can be detected as mass gap objects using future third-generation ground-based GW detectors [Gupta (2019)].

Summary
The CoRe database provides an expanded set of NR simulations to extract remnant properties of a large number of BNS systems. A highlight of this dataset is the rich dynamical ejecta information, which lets us estimate the post-merger recoil due to ejecta -thus filling a gap in the literature. We find that these ejecta kicks are greater by an order of magnitude than recoils imparted due to GW emission. With these more accurate kick calculations, we aim to obtain a fit for a recoil velocity as a function of BNS parameters. This fit will be useful for population modeling of compact binaries, and in particular for studying hierarchical mergers in dynamical environments.