\section{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 [\cite[Berti (2007)]{Berti_07},  \cite[Buonanno (2008)]{Buonanno_etal08}, \cite[Gerosa (2016)]{Gerosa_etal16}]. 
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 [\cite[Coughlin (2019)]{Coughlin_19}]
but none for the remnant's recoil velocity, mainly because of the additional uncertainty in the computation of the ejecta properties [\cite[Endrizzi (2016)]{Endrizzi_16}].
Therefore, we are estimating recoil velocities of BNS merger remnants using Computational Relativity (CoRe) database of numerical simulations [\cite[Deitrich (2018)]{Dietrich_etal2018}] and developing a fit for them covering an extended BNS parameter space. 






\section{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 \textit{kick}. 
These kicks have been shown to be particularly large ($\sim10^{2}-10^{3}~\text{km/s}$) for BBHs with unequal masses and high spins 
[see, e.g., \cite[Gerosa (2018)]{Gerosa2018Surrkick}]. 



\subsection{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 $\dot{h}(t)$, 
the time derivative of the GW strain expanded in the basis of spin-weighted spherical harmonics,  
as given in~\cite[ Ruiz (2008), Gerosa (2018)]{Ruiz_etal08,Gerosa2018Surrkick}

\begin{align}
\label{eq:dpdt1}
\frac{dP_{+}}{dt} &= \lim_{r\to\infty}\frac{r^{2}}{8\pi} \sum_{l,m}\dot{h}^{l,m} \Big[a_{l,m}\dot{h}^{l,m+1} + b_{l,-m}\dot{h}^{l-1,m+1} - b_{l+1,m+1}\dot{h}^{l+1,m+1} \Big], \\
\label{eq:dpdt2}
\frac{dP_{z}}{dt} &= \lim_{r\to\infty}\frac{r^{2}}{16\pi} \sum_{l,m}\dot{h}^{l,m} \Big[c_{l,m}\dot{h}^{l,m} + d_{l,m}\dot{h}^{l-1,m} - b_{l+1,m}\dot{h}^{l+1,m} \Big],
\end{align}
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.~(\ref{eq:dpdt1}) and (\ref{eq:dpdt2}) 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.

 
\subsection{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 ($\tilde \Lambda$) of the system as well as the equation of state [\cite[Nedora (2021)]{Nedora_etal21}]. 
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:

\begin{equation}
P_{\rm ej} = \int_{M_{\rm ej}} \rho |v| dV,
\end{equation}
where the volume integration is carried over $M_{\rm ej}$, the relevant ejecta mass, while $\rho$ and $v$ are the mass density and velocity respectively. Conservation of linear momentum then implies~\cite[Janka (2017)]{Janka_17}
\begin{equation}
v_{\rm rec} = \frac{P_{\rm ej}}{M_{\rm NS}},
\end{equation}
where $M_{\rm NS}$ is the mass of the remnant (with the net ejecta mass and radiated GW energy subtracted), and $v_{\rm rec}$ is the kick imparted due to ejecta emission.

\noindent\begin{minipage}{.48\textwidth}
        \centering
        \includegraphics[width=0.99\linewidth]{GWkicks.pdf}
    \end{minipage}\hfill
    \begin{minipage}{0.48\textwidth}
        \centering
        \includegraphics[width=0.99\linewidth]{Ejectakicks.pdf}
    \end{minipage}
\captionof{figure}[BNS remnant recoil velocity estimates]{\label{fig:kickestimates} Estimates of the kick magnitude due to GW emission (left) and Ejecta emission (right) plotted against the mass ratio of available NR simulations. Colors represent the different equations of state of neutron stars.}

\section{Results}
Figure \ref{fig:kickestimates} 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~  \text{km/s}$, and that for the ejecta emission is $28.9~  \text{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 $\sim100~\text{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. 


\subsection{Recoil velocity fit using Gaussian Process Regression}

\begin{wrapfigure}[17]{r}{0.5\textwidth}
  \begin{center}
    \includegraphics[width=0.49\textwidth]{GPR.pdf}
    \caption[GPR]{Example of a GPR fit for remnant kicks across a limited range of BNS mass ratios. \label{fig:gpr}}
  \end{center}
\end{wrapfigure}

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 ($\tilde \Lambda$).
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 \ref{fig:gpr} shows an example 1-parameter fit for recoils across BNS mass ratios using GPR.

 

\section{Astrophysical Implications}
The two main channels for the formation of a compact binary system are the \textit{isolated} channel wherein components of a massive main-sequence stellar binary individually collapse into compact objects, and the \textit{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 \textit{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-5 M_{\odot}$ 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_{\rm esc})$, the remnant recoil velocity $(v_{\rm rec})$ determines whether it would be retained ($v_{\rm rec} < v_{\rm esc}$) or  whether it would escape the environment ($v_{\rm rec} > v_{\rm 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_{\rm esc}$ are: (i) Open Clusters $< 10 ~\text{km/s}$, (ii) Globular Clusters $\sim 50~\text{km/s}$, and (iii) Nuclear Star Clusters $\sim 200~\text{km/s}$~\cite[Mapelli (2020)]{Mapelli_20}. 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 \textit{mass gap} objects using future third-generation ground-based GW detectors [\cite[Gupta (2018)]{Gupta_etal18}].


\section{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.

\medskip
\noindent
\paragraph{Acknowledgments: We thank David Radice and Rahul Kashyap 
for fruitful discussions and technical inputs for this project.}

