Accurate models for recoil velocity distribution in black hole mergers with comparable to extreme mass-ratios and their astrophysical implications

TL;DR

This work addresses the challenge of predicting black-hole remnant recoil kicks across the full mass-ratio spectrum by marrying analytic insight with data-driven methods. It introduces two aligned-spin kick models valid up to $q\in[1,200]$—an analytic model gwModel_kick_q200 and a GPR-based gwModel_kick_q200_GPR—and a probabilistic precessing-spin model gwModel_kick_prec_flow built with a normalizing flow, trained on NR and BHPT data. The aligned-spin models achieve high accuracy and robust extrapolation, with the analytic version favored for its smooth behavior, while the precessing-spin flow model reproduces the full distribution of kicks and remains efficient for population synthesis, even when extrapolated to extreme mass ratios. Together, these models improve upon existing surrogate and HLZ approaches, enable reliable hierarchical-merger and retention studies in varied environments, and are publicly available in the gwModels package, facilitating broader astrophysical applications across $q$ and spin configurations.

Abstract

Modeling the remnant recoil velocity (kick) distribution from binary black hole mergers is crucial for understanding hierarchical mergers in active galactic nuclei or globular clusters. Existing analytic models often show large discrepancies with numerical relativity (NR) data, while data-driven models are limited to mass ratios of q<=8 (aligned spins) and q<=4 (precessing spins) and break down when extrapolated outside their training ranges. Using ~5000 of NR simulations from the SXS and RIT catalogs up to q=128 and ~100 black hole perturbation theory simulations up to q=200, we present two classes of models: (i) gwModel_kick_q200 (gwModel_kick_q200_GPR), an analytic (Gaussian process regression) model for aligned-spin binaries. (ii) gwModel_kick_prec_flow, a normalizing-flow model for kick distribution from precessing binaries with isotropic spins. Our approach combines analytic insights from post-Newtonian theory with data-driven techniques to ensure correct limiting behavior and high accuracy across parameter space. Both gwModel_kick_q200 and gwModel_kick_prec_flow are valid from comparable to extreme mass ratios. Extensive validation shows all three models outperform existing ones within their respective domains. Finally, using both back-of-the-envelope estimates and 1404 detailed star cluster simulations incorporating our kick models, we find that the black hole retention probability in low mass globular clusters can vary noticeably when the gwModel_kick_prec_flow model is employed. The models are publicly available through the gwModels package.

Figures (12)

• Figure 1: We show the kick-velocity predictions from the analytic gwModel_kick_q200 model (blue solid line) along with the corresponding $95\%$ confidence interval (blue shaded region) in the nonspinning limit, as a function of mass ratio from $q = 1$ to $q = 256$. For comparison, we include a combination of NR data from the SXS (circles) and RIT (squares) catalogs, as well as BHPT results (triangles). The red, orange, purple, and green lines correspond to predictions from the NRSur7dq4Remnant, NRSur3dq8Remnant, BHPTNRSurRemnant, and HLZ models, respectively, all shown in their nonspinning limits. The red and orange dashed lines indicate extrapolations of the NRSur7dq4Remnant and NRSur3dq8Remnant models beyond their region of validity. The inset shows the residuals of the gwModel_kick_q200 predictions. More details are in Section \ref{['sec:aligned_spin_accuracy']}.
• Figure 2: (Upper panel) We show the histogram of errors in kick-velocity predictions from gwModel_kick_q200 (blue) and gwModel_kick_q200_GPR (black), computed with respect to NR and BHPT results for aligned-spin binaries across five different models in their aligned-spin limit. The NR error (gray) is estimated by comparing simulations at two different numerical resolutions when available. We find that gwModel_kick_q200 performs as well as, or in some cases better than, the NR surrogate models NRSur7dq4Remnant (orange) and NRSur3dq8Remnant (red) within their respective training domains, and shows significantly improved agreement with the data compared to the HLZ model (green). The errors from our model are comparable to the intrinsic NR uncertainties. (Lower panel) We show the kick-velocity distribution as a function of the mass ratio $q$ and the effective inspiral spin $\chi_{\mathrm{eff}}$. Dashed lines indicate regions where additional data are available from SXS NR, RIT NR, and BHPT simulations. More details are in Section \ref{['sec:aligned_spin_accuracy']}.
• Figure 3: We show the predictions from gwModel_kick_q200 (solid lines) and gwModel_kick_q200_GPR (dashed lines) for five different aligned-spin configurations as a function of mass ratio, ranging from $q = 1$ to $q = 256$. We find that the analytic model gwModel_kick_q200 exhibits smooth and physically consistent behavior across the full mass-ratio range, with the kick velocity decreasing monotonically with increasing $q$. In contrast, the data-driven gwModel_kick_q200_GPR model produces several unphysical features, even though it achieves a slightly lower validation error. Similar unphysical trends are also observed in one-dimensional parameter-space slices for other GPR-based models, such as NRSur7dq4Remnant and NRSur3dq8Remnant. An example of such behavior for NRSur3dq8Remnant is shown in Fig. \ref{['fig:nospin']}. More details are in Section \ref{['sec:aligned_spin_accuracy']}.
• Figure 4: We show the kick-velocity distribution from precessing-spin SXS (blue stars) and RIT (orange circles) NR simulations up to $q = 15$ for different spin magnitudes and orientation configurations. Our training dataset also includes BHPT simulations spanning mass ratios from $q = 40$ to $q = 100$ (not shown). For reference, we plot the empirical upper envelope of the kick velocities as a dashed gray line, given by $V_{\mathrm{upper}} = 7.6 \times 10^{4} \eta^{2}$. The inset shows the distribution of these data points as a function of the spin magnitudes $|\chi_1|$ and $|\chi_2|$. More details are in Section \ref{['sec:flow_model_accuracy']}.
• Figure 5: We show the training ($p_{\rm loss}^{\rm train}$; crimson solid line) and validation ($p_{\rm loss}^{\rm val}$; green dashed line) losses for the precessing-spin normalizing-flow model gwModel_kick_prec_flow as a function of training steps. Additionally, we plot the Wasserstein distance ($W_p$; black solid line) between the precessing-spin training data and the corresponding gwModel_kick_prec_flow outputs at various training stages. The inset shows the distribution of these data points as a histogram (gray), along with the corresponding distribution predicted by the trained normalizing-flow model gwModel_kick_prec_flow (blue). More details are in Section \ref{['sec:flow_model_accuracy']}.