Spin–Orbit Alignment in Merging Binary Black Holes Following Collisions with Massive Stars

Observational breakthroughs over the past decade have provided compelling evidence for the presence of black holes (BHs) in globular clusters (GCs). X-ray and radio detections (J. Strader et al. 2012), along with radial velocity measurements (B. Giesers et al. 2018), have confirmed that BHs can remain bound to GCs, contrary to previous theoretical predictions that dynamical interactions would eject them early in the cluster's evolution. Complementing these observations, advanced N-body simulations (e.g., M. Morscher et al. 2015; L. Wang et al. 2016; A. Askar et al. 2018; K. Kremer et al. 2020c) have demonstrated that BHs influence the long-term dynamics and structure of GCs, particularly in their dense cores. It is now widely accepted that most old GCs harbor dozens to hundreds of dynamically active BHs today, playing a pivotal role in delaying gravothermal collapse and shaping the cluster's overall evolution (D. Merritt et al. 2004; A. D. Mackey et al. 2007; P. G. Breen & D. C. Heggie 2013; K. Kremer et al. 2018, 2019a, 2020a).

Additionally, recent gravitational-wave (GW) detections by LIGO/Virgo/KAGRA (LVK) have provided groundbreaking evidence for the formation and merger of stellar-mass binary BHs (BBHs; R. Abbott et al. 2023). Although the formation pathways of these mergers remain an active area of research (I. Mandel & F. S. Broekgaarden 2022), dynamical formation in dense stellar systems, such as GCs, has emerged as a likely key contributor (e.g., C. L. Rodriguez et al. 2016a), competing with traditional isolated binary evolution scenarios (e.g., K. Belczynski et al. 2016). These dense environments provide the necessary conditions to drive frequent close encounters, enabling the formation, hardening, and eventual merger of BBHs.

Various recent studies have highlighted several key properties of the BBH population expected to form in clusters, including their masses (e.g., C. L. Rodriguez et al. 2019; K. Kremer et al. 2020a), eccentricities (e.g., M. Zevin et al. 2021), and spins (e.g., E. Payne et al. 2024). If observed, these features may provide crucial tests of the role of clusters in forming GW sources. Regarding spins in particular, the spin magnitude and orientation of BBH components are determined by the angular momentum accreted from their progenitors and their interactions with nearby astrophysical objects. Canonically, random interactions in dense environments have been linked to uncorrelated spin–orbit orientations and an isotropic distribution in effective spin (e.g., C. L. Rodriguez et al. 2016b), which is defined by

$$\begin{eqnarray}&&{\chi }_{{\rm{eff}}}=\frac{{M}_{1}{\chi }_{1}\cos {\theta }_{1}+{M}_{2}{\chi }_{2}\cos {\theta }_{2}}{{M}_{1}+{M}_{2}}\in (-1,1),\end{eqnarray} \tag{ 1 }$$

where M_i is the BH mass; ${\chi }_{i}\equiv {S}_{i}/{M}_{i}^{2}\in (0,1)$ denotes the dimensionless spin of each BH, with S_i representing the spin angular momentum of the BH; and θ_i is the angle between the vector S_i and the orbital angular momentum. However, evidence suggesting a preference against a purely isotropic spin distribution (R. Abbott et al. 2023) poses a challenge to the hypothesis that BBHs form predominantly in dense star clusters through random dynamical encounters (e.g., C. L. Rodriguez et al. 2021).

Recent studies reveal that BHs in clusters often undergo close encounters with stars, leading to unique electromagnetic (EM) transients and, in some cases, significantly altering BH properties. These encounters may cause a star to pass by a BH within its tidal disruption radius or to collide physically with the BH, resulting in the complete or partial disruption of the star (e.g., H. B. Perets et al. 2016; K. Kremer et al. 2019b, 2022, 2023; M. Lopez et al. 2019; T. Ryu et al. 2020; F. Kıroğlu et al. 2023). Our recent work demonstrates that collisions with massive stars, comparable in mass to the BHs, can result in substantial accretion and spin-up, with important implications for GW detections by LVK (F. Kıroğlu et al. 2025).

N-body simulations show that the majority of BH–star collisions occur during multibody encounters (J. M. Fregeau & F. A. Rasio 2007; K. Kremer et al. 2021; F. Kıroğlu et al. 2025). This trend is particularly pronounced in clusters with a high initial binary fraction among massive stars (e.g., E. González et al. 2021), consistent with observational data indicating that nearly all O- and B-type stars in the Galactic field are born in binaries (H. Sana et al. 2012). Motivated by these trends, in this Letter, we conduct hydrodynamic simulations to study encounters between BBHs and massive stars. We explore a portion of the parameter space for BBH properties (e.g., semimajor axis, eccentricity, and component masses) guided by N-body simulations of dense star clusters.

Our focus is on BBHs that are initially wide enough to interact with stars before merging via GW emission but also compact enough that both BHs can dynamically interact with the disrupted stellar debris. Figure 1 sets the stage for our Letter, where we demonstrate the parameter space for the BBH semimajor axis and total mass of interest. In gray, we highlight the "excluded" region where the BBHs are likely to merge before undergoing a stellar encounter that may lead to a collision. We calculate the GW inspiral time t_GW of the binaries following P. C. Peters (1964) and the encounter timescale using ${t}_{{\rm{enc}}}={(n{{\rm{\Sigma }}}_{{\rm{bs}}}{\sigma }_{v})}^{-1}$, where n is the cluster's central number density, σ_v is the central velocity dispersion, and Σ_bs is the cross section for binary–single encounters. The scatter points represent all BBH–star collisions from our recent N-body simulations of dense star clusters (F. Kıroğlu et al. 2025) using the Cluster Monte Carlo code, a Hénon-style N-body code for stellar dynamics (see C. L. Rodriguez et al. 2022 for a detailed review). These simulations are performed with an initial total number of stars N = 8 × 10⁵, metallicity Z = [0.1, 1] Z_⊙, an initial virial radius r_v = [0.5, 1] pc, and a primordial binary fraction for massive stars f_b (>15 M_⊙) = 100%. In these models, up to 20% of all stellar disruptions by BHs occur during BBH+star three-body encounters, potentially enhancing the merger rates of these BBHs through orbital parameter modifications. We find that approximately 60% of these BBHs are dynamically assembled, while the remaining 40% originate from dynamically shaped primordial binaries.

Zoom In Zoom Out Reset image size

Our Letter is organized as follows. In Section 2, we describe the computational methods used for hydrodynamic simulations and the calculation of BH spins resulting from stellar collisions and subsequent accretion. In Section 3, we present the outcomes of the hydrodynamic simulations and discuss the varying results across the parameter space. Finally, we conclude and summarize our findings in Section 4.

Hydrodynamic simulations of the BBH+star collision parameter space mapped out in Figure 1 are the focus of this study. For this task, we use the StarSmasher code (F. A. Rasio & S. L. Shapiro 1991; E. Gaburov et al. 2018), which employs smoothed particle hydrodynamics (SPH) to model stellar gas as particles with density profiles defined by a smoothing kernel. We adopt a Wendland C4 kernel (H. Wendland 1995) with compact support for smoothing and gravitational softening. To prevent unphysical particle interpenetration, we implement artificial viscosity with a Balsara switch (D. S. Balsara 1995) as described in J. C. L. Hwang et al. (2015). Gravitational forces are computed using direct summation on NVIDIA GPUs, which provides improved accuracy compared to tree-based methods (E. Gaburov et al. 2010). The simulations use an equation of state that incorporates contributions from both ideal gas and radiation pressure (J. C. Lombardi et al. 2006). The SPH particles evolve using variational equations of motion (J. J. Monaghan 2002; V. Springel & L. Hernquist 2002; J. C. Lombardi et al. 2006).

2.1. Stellar Profiles

2.2. Initial Conditions and Parameters

2.3. Analysis of Hydrodynamics

In postprocessing, we compute at each time snapshot the mass bound to each of four components: the two individual BHs, the star (if it survives the collision), and the common envelope that forms around the BBH. We employ an iterative procedure similar to that of J. C. Lombardi et al. (2006) and K. Kremer et al. (2022). A necessary condition for a particle to be bound to one of the components is that it has negative specific mechanical energy relative to the center of mass of that component. In particular, the specific mechanical energy of particle i with respect to the center of mass of component j is calculated as

$$\begin{eqnarray}&&\frac{{v}_{ij}^{2}}{2}-\frac{G{M}_{j}}{{d}_{ij}},\end{eqnarray} \tag{ 2 }$$

where v_ij is the speed of particle i relative to the center of component j, G is Newton's gravitational constant, M_j is the mass of component j, and d_ij is the distance between the particle and the center of component j. We choose not to include the specific internal energy u_i in the boundness criterion, the preferred approach of J. L. A. Nandez et al. (2014). If the specific mechanical energy is negative for either the BH or the star, the particle is bound to the corresponding component (or possibly to the common envelope instead, as discussed below); if it is negative for more than one component, it is bound to the one for which it is most negative.

The iterative nature of this procedure arises from the need to refine the gravitational potential used in Equation (2). Initially, at each snapshot, we assume the masses bound to the components are the same as those from the previous snapshot. We then calculate the specific mechanical energy of each particle and update the set of bound particles for each component. These updated bound masses are then used to refine the gravitational potential and the specific mechanical energy calculations. The process repeats until the set of bound masses converges, ensuring self-consistency between the calculated bound particles and the gravitational potential associated with them.

Unlike in K. Kremer et al. (2022), where each simulation contained only a single BH, the presence of a BBH requires that we consider a common envelope. There are two ways by which a particle could be counted as part of the common envelope. First, particles that are not bound to either BH individually but are bound to the center of mass of the binary system are considered part of the common envelope. This binding criterion uses the center-of-mass position of the binary and the total mass of the two BHs to determine whether a particle has negative specific mechanical energy relative to the center of mass of the BBH. Second, particles may be considered bound to the common envelope if they lie outside the Roche lobe of either BH but still have negative specific energy with respect to one or both BHs individually. The Roche lobe radius for each BH is calculated using the standard Eggleton approximation (P. P. Eggleton 1983), with the separation between the two BHs taken as the periapsis separation of the BHs; this criterion ensures that the mass bound to each BH remains stable over an orbit. Particles that are not bound to either of the BHs or the common envelope are considered part of the ejecta.

Asymmetric mass loss during the collision imparts a net momentum to the BBH, causing it to drift within the center-of-mass frame in which the calculation is performed. To quantify this kick, we calculate the center-of-mass velocity of the BBH, including the mass bound to it, at the final simulation snapshot, defining this velocity as v_kick.

We note that the calculation of the orbital parameters by the postprocessing code treats the BHs, including the mass bound to them, in the two-body approximation. This simplification occasionally leads to an assignment of a BBH eccentricity e > 1 during the disruption of the star (see Figure 2). This temporary artifact arises when the relative velocity of the BHs exceeds the escape velocity calculated in a two-body framework, that is, calculated without regard for the existence of the star (or common envelope). In three-body interactions, the star can boost the BHs' relative speed, resulting in values of eccentricity e > 1 being calculated even though the BBHs will ultimately remain as a binary. The eccentricity can then drop below 1 as the star or stellar debris repositions itself and decreases the relative speed of the two BHs. The eccentricity values being shown are most reliable once the star has been fully disrupted and the system is well modeled as two bodies.

Zoom In Zoom Out Reset image size

2.4. Mass Accretion and Spin-up

In order to compute χ_eff (Equation (1)) for each BBH, we must define how the spin vector of each BH evolves following this accretion. We calculate the spin attained by an initially nonrotating BH with mass M due to accretion of the disrupted material following J. M. Bardeen et al. (1972) and K. S. Thorne (1974),

$$\begin{eqnarray}&&\chi (M,{M}_{{\rm{f}}})={\left(\frac{2}{3}\right)}^{1/2}\frac{M}{{M}_{{\rm{f}}}}\,\left[4-{\left(\frac{18{M}^{2}}{{M}_{{\rm{f}}}^{2}}-2\right)}^{1/2}\right],\end{eqnarray} \tag{ 3 }$$

where M_f is the final BH mass after accreting. This prescription assumes that the specific angular momentum accreted is at most the angular momentum of the last stable orbit. The final BH mass is derived from our hydrodynamic simulations, under the assumption that all material determined to be bound to each BH is eventually accreted. Our estimates therefore represent an upper limit to the BH accretion and spin values. We assume that the direction of each BH spin is in the same direction as the final angular momentum of the disk.

Although in reality, accretion and spin-up occur on longer timescales than our SPH simulations, it is also useful to define an "instantaneous" χ_eff value, calculated from the angular momentum of stellar debris bound to each BH at any given time. This parameter evolves as both the amount of material and its orientation change over time and helps inform how the angular momentum is exchanged between the BBH and the gas over the course of the encounter. This is distinct from the final effective spin parameter χ_eff,f (reported in column (18) of Table 2) measured at the end of our simulations, once the system's evolution has stabilized.

In this section, we describe the results of hydrodynamic simulations of BBH–star collisions and explore their effects on the BBH orbital parameters and the BH spins, both in direction and in magnitude. We provide a list of all hydrodynamic simulations performed in this study in Table 2, including the initial conditions and outcomes.

In the region of the parameter space we consider, where a_i ≲ 0.2 au, each BH in a binary acquires a rotating envelope following the disruption of the star. This outcome can occur under two distinct scenarios. In the first case, the star is fully disrupted by one BH, and a subfraction of the stellar debris is captured and bound to the nondisrupting BH. This outcome is equivalent to the "overflow scenario" described in M. Lopez et al. (2019). The second case involves multiple pericenter passages, where the star survives the first encounter with one BH and subsequently undergoes additional pericenter passage(s) around both BHs, ultimately leading to its complete disruption. For very wide BBH orbits, i.e., a_i ≳ 1 au, the collision is expected to resemble a single BH–star interaction with only the disrupting BH accreting material (M. Lopez et al. 2019). In this regime (where a_i ≫ R_⋆), binary torques cannot influence the angular momentum of the disks, which is determined by the angular momentum of the star at the time of disruption. Consequently, these cases are not of interest for this study.

In Figure 2, we present the evolution of the eccentricity, instantaneous χ_eff, and separation between the BHs following their collision with a 10 M_⊙ star. The four colors show separate simulations with different initial semimajor axes. For the binaries with a_i ≥ 0.1 au, the instantaneous χ_eff is negative after the initial disruption of the star. We quantify this phase with the parameter t_χ<0 (shown in column (17) in Table 2), which is the longest consecutive time interval that the instantaneous χ_eff is less than 0. As the BBH undergoes subsequent pericenter passages, the binary exerts a torque on the misaligned stellar debris and ultimately tilts its orientation into alignment with the BBH orbit. In column (16) of Table 2, we list the time (typically around 1−10 days) between the disruption of the star and the next pericenter passage of the BHs. As tidal torques that lead to angular momentum exchange are most prominent at the pericenter, this timescale is typically comparable to the timescale for the angular momentum of the stellar debris to align with the orbit.

In Figure 3, we show a series of SPH time snapshots for the a_i = 0.1 au model of Figure 2 (model 13 in Table 2). Here, color denotes the mass density of stellar material, and arrows denote velocity vectors of the fluid elements in the plane shown. As notated in the various panels, the instantaneous χ_eff is initially negative but ultimately flips to a positive value as angular momentum exchanges from the BBH orbit to the gas. In Figure 4, we show trajectories for the two BHs in this same simulation, with colors representing the instantaneous χ_eff. Again, we see that after multiple pericenter passages, each providing a significant tidal torque, the instantaneous χ_eff ultimately flips from negative to positive.

Zoom In Zoom Out Reset image size

Conversely, for the binary in Figure 2 with an initial semimajor axis of a_i = 0.02 au, the stellar debris is aligned with the orbit after the initial disruption of the star. In this case, the instantaneous χ_eff is positive at all times (t_χ<0 = 0 in Table 2). These collective findings show that spin–orbit alignment is ultimately achieved in all BBH–star collisions involving stars and BBHs with the mass ratios and orbital properties explored in this work, where the BBH orbit is comparable to the characteristic size of the envelope, roughly R_⋆.

In most of the simulations performed (nearly head-on encounters), the initial angular momentum of the star's orbit is less than that of the BBH orbit. Consequently, the stellar debris field will tend to align with the orbital direction of the BBH. However, in cases where the angular momentum in the star's orbit is large (grazing encounters; models 28 and 31 in Table 2), angular momentum exchange occurs in the opposite direction: the orbital direction of the BBH will tend to align with the debris field. This follows a general principle that the angular momentum of the orbit and the debris field ultimately tend toward alignment, regardless of which component contributes more significantly.

In Figure 5, we present the final inspiral timescale versus the final effective spin χ_eff,f of the postcollision BBHs for all SPH simulations. In response to the stellar collision, we find that BBH orbits can decrease by up to 10 times (the ratio of final to initial semimajor axis, a_f/a_i, is shown in column (14) of Table 2) and attain moderate eccentricity e ≈ 0.5. Both of these factors reduce their GW inspiral timescales significantly. As a result, many of the postcollision BBHs shown in Figure 5 are likely to merge with positive χ_eff before subsequent dynamical encounters can randomize their spin–orbit orientation.

Zoom In Zoom Out Reset image size

A fraction of the material from the disrupted star becomes bound to the BHs, while a fraction, with positive total energy, becomes unbound from the system entirely. This unbound ejecta imparts a recoil kick to the center of mass of the BBH. As shown in column (12) of Table 2, the amount of unbound material can be comparable to the BH mass.

The cases with high kick velocities are associated with encounters that have a large impact velocity and a direct, close to head-on collision between the star and one of the BHs. Such cases result in asymmetric mass ejection that gives a significant kick to the BBH. In the context of a star colliding with a single BH, this kick is studied in some detail in K. Kremer et al. (2022, see Section 3.2). The main difference here is that the kick is ultimately given to the BBH as a whole instead of a single, isolated BH. In cases when the impact velocity is small or the encounter with the BH is less direct, much of the mass ejection occurs during the subsequent common-envelope-like phase of the BBH, leading to more axisymmetric ejection and a smaller final kick.

The resulting recoil velocities can reach up to 80 km s⁻¹, potentially sufficient to eject the BBH from its host cluster core or, in extreme cases, from its host cluster entirely. In these cases, even the BBHs with relatively long postcollision inspiral times (t_GW ≳ 100 Myr) may still retain positive χ_eff values at merger, since the chance of a subsequent encounter is reduced in the lower-density cluster halo (or eliminated fully in the case of ejection from the cluster). The kick velocities of BBHs after each collision are illustrated in the color bar of Figure 5.

BBHs formed either dynamically or through binary evolution in dense star clusters are expected to undergo a series of interactions, including collisions and mergers with stars. These interactions can significantly alter their mass, spin, and orbital properties, thereby shaping their GW signatures during inspiral and merger. To investigate these processes, we performed a suite of SPH simulations of close encounters between massive stars and BBHs typical of those expected to occur in young and dense cluster environments.

Previous studies have explored the impact of stellar tidal disruption events on BBH spin–orbit orientations (M. Lopez et al. 2019; T. Ryu et al. 2022). However, these studies have focused mainly on the disruption of 1 M_⊙ MS stars, motivated by both the higher abundance of low-mass stars in clusters and the shorter lifetimes of more massive stars. In this work, we explore interactions with massive (10 M_⊙) stars that lead to direct collisions with the BBH components. We obtain qualitatively different results, as the impact on the binary orbit and mass accretion is greater in nearly equal-mass encounters. We find that collisions between BBHs and massive stars result in the preferential alignment of the BH component spins with the orbital angular momentum, provided the binaries are initially sufficiently compact (≲1 au) to form individual accretion disks that can interact with and be torqued by the orbit during the close pericenter passage. On the other hand, T. Ryu et al. (2022) found that initially circular BH binaries experience only a ≈10% increase in their eccentricity following the disruption of a 1 M_⊙ star. In that case, the periapsis passage may never get close enough to induce significant torques for the disk to realign with the orbital angular momentum.

Our N-body cluster simulations demonstrate that a typical cluster features approximately 10 (3) BBH+star collisions with mass ratios q > 0.1 (q > 0.5) at early times (see Figure 1). Consequently, for a typical cluster with roughly 100 BBH mergers over its full ∼12 Gyr lifetime (e.g., K. Kremer et al. 2020c), we expect roughly 10% of the BBH mergers to be affected by these collisions. Following angular momentum exchange and ultimately accretion, we find that these collisions tend to result in spin–orbit alignment of the BBH, with important implications for the BBH's effective spin parameter χ_eff at the time of merger. While canonical gas-free cluster dynamics provides the best route for forming binaries with a significant component of the spins antialigned with the orbital angular momentum (e.g., E. Payne et al. 2024), BBHs undergoing an accretion phase after colliding with massive stars in young clusters could contribute toward the observed bias toward small positive χ_eff values seen in current LVK data. The fraction of BH mergers with aligned spins due to interactions with massive stars, as predicted in this study (≈10%), is roughly comparable to the fraction of mergers that are asymmetric around χ_eff = 0, as predicted by S. Banagiri et al. (2025).

A key limitation of the calculations presented here is the exclusion of accretion feedback effects on the hydrodynamics and ultimate accretion. Accretion releases energy over various timescales and may drive outflows, resulting in an uncertain fraction of mass loss. While our hydrodynamic models provide insight into the amount of stellar material bound to each BH, we assume 100% accretion efficiency of this material in our spin calculations, ignoring the potential impact of these outflows. Consequently, our estimates represent an upper limit to the BH spin values. Additionally, although our simulations are hydrodynamic, in reality, magnetic fields are present, which may drive accretion either via turbulence that acts like an effective viscosity or via winds that extract mass and angular momentum via torque.

For cases where the initial disks are misaligned with the orbit, we assume that the disks persist for at least 10 days, long enough for the next BH periapsis passage to occur. However, it is possible that the material in the misaligned disk might have already been accreted by the BH. In this case, the disk could remain misaligned long enough for the BH to accrete the antialigned angular momentum, possibly leading to a final negative effective spin parameter. The viscous accretion timescale remains highly uncertain, as it depends on the poorly constrained alpha viscosity and disk scale height. Further simulations may provide further insight into these quantities.

Last, if the common envelope around the BBH formed after the collision survives until the time of the BBH merger, it could produce an EM counterpart, resembling a weak, short gamma-ray burst (R. Perna et al. 2016). A detailed exploration of such counterparts associated with accreting BBH mergers in dense star clusters will be presented in a forthcoming paper (F. Kıroğlu et al. 2025, in preparation).

This work was supported by NSF grant AST-2108624 and NASA grant 80NSSC22K0722, as well as the computational resources and staff contributions provided for the Quest high-performance computing facility at Northwestern University. Quest is jointly supported by the Office of the Provost, the Office for Research, and Northwestern University Information Technology. F.K. acknowledges support from a CIERA Board of Visitors Graduate Fellowship. We thank Emma Chambers and Cade Snyder for their contributions to earlier related work that aided the development of analysis routines used in this research. This work made use of the SPLASH visualization software (D. J. Price 2007).

A.1. Initial Conditions for BH and Star Positions and Velocities

To facilitate comparisons with future works, we describe here the initial conditions for the positions and velocities of the BHs and the star. The overall strategy uses the analytic solution to the Kepler two-body problem (e.g., S. Thornton & J. Marion 2004) twice. The first usage of the Keplerian solution is for the BBH: we treat the BHs as a two-body system orbiting their common center of mass in an XY-plane. The orbit of the BBH is then rotated out of that plane at an inclination angle i, as described below. The second usage of the Keplerian solution is for the parabolic orbits of the star and the BBH about their common center of mass. To keep this center of mass of the three-body system at the origin, we shift the BBH and star in opposite directions away from the origin of the xyz-frame used for the SPH calculation.

To initiate the BBH, we set values for the following parameters: the masses of the two BHs (M_BH1 = M_BH2 = M_BBH/2), the semimajor axis a, the orbital eccentricity e, the inclination i relative to the orbital plane of the star, and the true anomaly f. In all of our simulations, the initial argument of periapsis and longitude of the ascending node are zero. An XY-coordinate system is used to define, temporarily, the BBH orbital dynamics. In our convention, a true anomaly f = 0 corresponds to the second BH being on the +X-axis. Because we consider equal-mass BHs, the position coordinates of the first BH in the XY-frame are simply ${X}_{1}=-\frac{R}{2}\cos f$ and ${Y}_{1}=-\frac{R}{2}\sin f$, where the separation of the BH components is

$$\begin{eqnarray}&&R=\displaystyle \frac{a(1-{e}^{2})}{1+e\cos f}.\end{eqnarray} \tag{ A1 }$$

The position of the second BH is then chosen so that the center of mass of the binary is at the origin in the XY-frame: X₂ = −X₁, Y₂ = −Y₁. The velocity components of the first BH in this orbital plane are obtained by differentiating X₁ and Y₁ with respect to time, treating the semimajor axis a and eccentricity e as constants, with

$$\begin{eqnarray}&&\dot{f}={(1+e\cos f)}^{2}\sqrt{\displaystyle \frac{G{M}_{\mathrm{BBH}}}{{a}^{3}{(1-{e}^{2})}^{3}}}\end{eqnarray} \tag{ A2 }$$

as the initial rate of change of the true anomaly and

$$\begin{eqnarray}&&\dot{R}=e\sin f\sqrt{\displaystyle \frac{G{M}_{\mathrm{BBH}}}{a(1-{e}^{2})}}\end{eqnarray} \tag{ A3 }$$

as the initial rate of change of the separation. The velocity components of the second BH are similarly set but with the opposite sign so that the total momentum of the BBH in the XY-frame is zero.

For the initial conditions of the star, the BBH is treated as a single point mass at its center of mass, and both the BBH and the star follow parabolic orbits about their common center of mass. Their separation is given by $r=2{r}_{p}/(1-\cos \theta )$, where the angular position θ of the star is measured counterclockwise with respect to the +x-axis in the xy-plane and the periapsis separation r_p occurs at θ = π. The parabolic trajectory of the star is aligned such that it approaches the BBH from infinity in the first quadrant and would recede toward infinity in the fourth quadrant. (When r_p = 0, the star approaches from along the +x-axis.) To set the initial position of the star, we determine its coordinates for the initial separation r = r₀ = 50 R_⊙ and the initial angular position $\theta ={\theta }_{0}={\cos }^{-1}(1-2{r}_{p}/{r}_{0})$. In particular, so that the center of mass of the entire three-body problem is at the origin in the xyz-frame, the star is placed with position coordinates ${x}_{3}={(1+Q)}^{-1}r\cos \theta$, ${y}_{3}\,={(1+Q)}^{-1}r\sin \theta$, and z₃ = 0, while the initial center-of-mass position of the BBH is x_BBH = −Qx₃, y_BBH = −Qy₃, z_BBH = 0. Here the mass ratio Q = M₃/M_BBH, and M₃ is the star mass. We note that ${x}_{3}-{x}_{{\rm{BBH}}}=r\cos \theta$ and ${y}_{3}-{y}_{{\rm{BBH}}}=r\sin \theta$.

The initial velocity components of the star are obtained by differentiating x₃ and y₃ with respect to time and making use of $\dot{\theta }=\sqrt{2G{r}_{p}{M}_{{\rm{tot}}}}/{r}^{2}$ and $\dot{r}=-\sqrt{2G{M}_{{\rm{tot}}}(1-{r}_{p}/r)/r}$ (from angular momentum and energy conservation), where M_tot = M_BBH + M₃. The center-of-mass velocity of the BBH is then chosen so that the center of mass of the entire three-body system will remain at the origin in the xyz-frame.

The 3D positions of the binary components, accounting for the center-of-mass shift and the inclination i, are

$$\begin{eqnarray}&&\left[\begin{array}{l}{x}_{j}\\ {y}_{j}\\ {z}_{j}\end{array}\right]=\left[\begin{array}{l}{x}_{\mathrm{BBH}}\\ {y}_{\mathrm{BBH}}\\ 0\end{array}\right]+\left[\begin{array}{rcl}\cos i & 0 & -\sin i\\ 0 & 1 & 0\\ \sin i & 0 & \cos i\end{array}\right]\left[\begin{array}{l}{X}_{j}\\ {Y}_{j}\\ 0\end{array}\right],\end{eqnarray} \tag{ A4 }$$

where j = 1, 2 corresponds to the two BHs and X_j, Y_j are the coordinates of the BHs in the orbital plane, with Z_j = 0. The transformation of the velocity components is of the same form as the transformation of the position coordinates. We note that i = 0 corresponds to a prograde approach, while i = 180° corresponds to a retrograde approach.

Putting everything together, the initial position of the star is

$$\begin{eqnarray}&&\left[\begin{array}{l}{x}_{3}\\ {y}_{3}\\ {z}_{3}\end{array}\right]={(1+Q)}^{-1}\left[\begin{array}{l}{r}_{0}-2{r}_{p}\\ 2{({r}_{0}{r}_{p}-{r}_{p}^{2})}^{1/2}\\ 0\end{array}\right],\end{eqnarray} \tag{ A5 }$$

while the initial positions of the BHs are

$$\begin{eqnarray}&&\left[\begin{array}{l}{x}_{j}\\ {y}_{j}\\ {z}_{j}\end{array}\right]=\left[\begin{array}{l}-Q{x}_{3}\mp \displaystyle \frac{R}{2}\cos i\cos f\\ -Q{y}_{3}\mp \displaystyle \frac{R}{2}\sin f\\ \mp \displaystyle \frac{R}{2}\sin i\cos f\end{array}\right],\end{eqnarray} \tag{ A6 }$$

where R is given by Equation (A1), the minus sign is chosen from the ∓ for BH j = 1, and the plus sign is chosen for BH j = 2. The initial velocity of the star is

$$\begin{eqnarray}&&\left[\begin{array}{l}{v}_{x3}\\ {v}_{y3}\\ {v}_{z3}\end{array}\right]=-{(1+Q)}^{-1}\sqrt{2G{M}_{\mathrm{tot}}}\left[\begin{array}{l}\sqrt{\left(1-{r}_{p}/{r}_{0}\right)/{r}_{0}}\\ \sqrt{{r}_{p}}/{r}_{0}\\ 0\end{array}\right],\end{eqnarray} \tag{ A7 }$$

while the initial velocities of the BHs are

$$\begin{eqnarray}&&\left[\begin{array}{l}{v}_{xj}\\ {v}_{yj}\\ {v}_{zj}\end{array}\right]=\left[\begin{array}{l}-Q{v}_{x3}\pm \displaystyle \frac{1}{2}(R\,\dot{f}\sin f-\dot{R}\cos f)\cos i\\ -Q{v}_{y3}\mp \displaystyle \frac{1}{2}(R\,\dot{f}\cos f+\dot{R}\sin f)\\ \pm \displaystyle \frac{1}{2}(R\,\dot{f}\sin f-\dot{R}\cos f)\sin i\end{array}\right],\end{eqnarray} \tag{ A8 }$$

where R, $\dot{f}$, and $\dot{R}$ come from Equations (A1)–(A3), the top sign is chosen from the ∓ or the ± for BH j = 1, and the bottom sign is chosen for BH j = 2. For example, in case 3, we have in code units (G = M_⊙ = R_⊙ = 1) Q = 0.5, a = 3.95, e = 0, f = 30°, r_p = 0, i = 0, M_BBH = 20, which yields initial positions for the star and BHs of (33.3, 0.0) and (−16.7 ∓ 1.7, ∓ 1.0), respectively, as well as initial velocities of (−0.73, 0) and (0.37 ± 0.56, ∓ 0.97). In case 25, we have in code units Q = 1/3, a = 21.5, e = 0.5, f = 0°, r_p = 0.5, i = 180°, M_BBH = 30, which yields initial positions for the star and BHs of (36.75, 7.46) and (−12.250 ∓ 5.375, −2.488), respectively, as well as initial velocities of (−0.944, −0.095) and (0.315, 0.032 ∓ 1.023). In Appendix Table 2, we present detailed information for each simulation. Figure 6 shows the density and mass profiles of the MESA stellar model used in our simulations.

Zoom In Zoom Out Reset image size

Table 2. Collision Outcomes and Postcollision Properties

Case	MBBH	Star	rp	ai	eBBH	f	i	tGW,i	MBH 1,f	MBH 2,f	mej	mCE	af/ai	ef	tdtp	tχ<0	χeff,f	tGW,f	vkick
	(M⊙)		(au)	(au)		(deg)	(deg)	(Myr)	(M⊙)	(M⊙)	(M⊙)	(M⊙)			(days)	(days)		(Myr)	(km s−1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)	(20)
1	20	10 MS	0.012	0.018	0.0	0	0	18	11.25	10.21	7.8	0.77	0.78	0.12	11	0	0.23	5.1	17
2	20	10 MS	0.000	0.018	0.0	0	0	18	10.20	10.07	8.9	0.84	0.33	0.48	0.47	0.64	0.05	0.081	13
3	20	10 MS	0.000	0.018	0.0	30	0	18	10.18	10.08	8.4	1.29	0.27	0.48	0.24	0.35	0.05	0.035	20
4	20	10 MS	0.000	0.018	0.0	0	30	18	10.01	10.00	9.6	0.43	0.21	0.73	0.3	0.33	0.00	0.0022	30
5	20	10 MS	0.000	0.018	0.0	0	180	18	10.07	10.18	9.0	0.75	0.32	0.49	0.46	0.62	0.04	0.074	14
6	20	10 MS	0.007	0.018	0.0	0	0	18	10.95	10.28	8.1	0.66	0.45	0.01	1.6	0.0092	0.20	0.64	41
7	20	10 MS	0.000	0.037	0.0	0	0	2.9e+02	10.19	10.59	8.5	0.72	0.19	0.05	0.36	0	0.13	0.36	54
8	20	10 MS	0.000	0.037	0.0	0	180	2.9e+02	10.95	10.19	5.6	3.24	0.26	0.31	0.083	0	0.18	0.83	55
9	20	10 MS	0.002	0.100	0.0	0	0	1.6e+04	10.07	11.40	6.6	1.97	0.66	0.50	21	27	0.22	8.9e+02	13
10	20	10 MS	0.000	0.100	0.0	30	0	1.6e+04	10.86	10.20	8.4	0.52	0.12	0.17	8.2	8.3	0.17	2.2	45
11	20	10 MS	0.000	0.100	0.0	60	0	1.6e+04	10.05	11.22	7.1	1.66	0.67	0.63	19	25	0.20	4.7e+02	5.1
12	20	10 MS	0.005	0.100	0.0	0	0	1.6e+04	10.05	12.44	5.0	2.52	1.5	0.65	88	89	0.33	8.9e+03	14
13	20	10 MS	0.009	0.100	0.0	0	0	1.6e+04	10.01	12.79	4.8	2.38	0.75	0.37	23	31	0.36	2e+03	24
14	20	10 MS	0.000	0.200	0.0	0	0	2.6e+05	10.01	11.46	4.4	4.13	0.36	0.46	37	39	0.22	1.4e+03	4.1
15	30	10 MS	0.000	0.018	0.0	0	0	5.4	15.48	15.38	6.8	2.35	0.47	0.28	0.17	0.16	0.10	0.17	34
16	30	10 MS	0.000	0.037	0.0	0	0	91	15.21	16.09	8.0	0.71	0.34	0.18	0.38	0.67	0.14	0.92	82
17	30	10 MS	0.000	0.037	0.0	0	30	91	15.15	15.25	9.1	0.53	0.34	0.46	0.35	0.055	0.05	0.52	23
18	30	10 MS	0.000	0.037	0.0	0	60	91	15.26	15.10	9.0	0.61	0.47	0.58	0.26	0.32	0.04	1	12
19	30	10 MS	0.000	0.037	0.0	0	90	91	15.14	15.15	9.0	0.75	0.36	0.41	0.3	0.14	0.03	0.78	16
20	30	10 MS	0.000	0.037	0.0	0	120	91	15.10	15.25	9.0	0.65	0.47	0.59	0.25	0.32	0.04	0.95	14
21	30	10 MS	0.000	0.037	0.0	0	150	91	15.26	15.15	9.0	0.55	0.35	0.47	0.36	0.065	0.05	0.55	26
22	30	10 MS	0.000	0.037	0.0	0	180	91	16.26	15.21	7.7	0.86	0.36	0.26	0.39	0.68	0.16	1	82
23	30	10 MS	0.000	0.074	0.0	0	0	1.5e+03	15.10	15.84	7.6	1.45	0.83	0.52	3.2	6.4	0.10	2.1e+02	13
24	30	10 MS	0.002	0.100	0.0	0	180	4.8e+03	15.00	15.45	8.9	0.70	0.15	0.38	7.5	0.68	0.05	1.3	17
25	30	10 MS	0.002	0.100	0.5	0	180	1.7e+03	15.13	15.27	9.0	0.58	0.2	0.67	0.26	0.014	0.05	0.94	6.5
26	30	10 MS	0.002	0.100	0.9	0	180	16	15.00	15.00	10.0	0.04	0.042	0.99	0.59	0.54	0.00	1.1e-08	65
27	30	10 MS	0.000	0.100	0.0	0	180	4.8e+03	15.08	15.61	8.5	0.81	0.76	0.82	8.3	3.8	0.08	31	8.7
28	30	10 MS	0.047	0.100	0.0	0	180	4.8e+03	15.00	15.00	10.0	0.03	0.037	0.85	0.6	0.6	0.00	0.00012	27
29	30	10 MS	0.012	0.100	0.0	0	180	4.8e+03	15.04	15.00	9.9	0.07	0.075	0.79	0.16	0.0092	0.00	0.005	17
30	30	10 MS	0.012	0.100	0.5	0	180	1.7e+03	15.00	15.00	9.9	0.05	0.057	0.87	0.36	0.35	0.00	0.00038	30
31	30	10 MS	0.012	0.100	0.9	0	180	16	15.00	15.06	9.9	0.02	0.06	0.81	15	0.18	0.01	0.0016	51
32	20	10 post	0.000	0.040	0.0	0	0	4.1e+02	10.92	10.08	8.3	0.30	0.41	0.05	0.99	1.6	0.16	9.8	22
33	20	10 post	0.000	0.100	0.0	0	0	1.6e+04	10.06	11.31	7.2	1.00	0.41	0.28	1.7	3.9	0.21	2.7e+02	25
34	20	10 post	0.000	0.150	0.0	0	0	8.1e+04	10.09	10.90	7.5	1.11	0.23	0.54	10	0.77	0.15	55	18
35	30	10 post	0.000	0.040	0.0	0	0	1.2e+02	15.66	15.18	8.1	0.65	0.55	0.30	0.52	0.65	0.10	7.1	11
36	30	10 post	0.000	0.100	0.0	0	0	4.8e+03	15.60	15.05	8.5	0.40	0.36	0.69	1.7	1.1	0.07	7.8	35
37	30	10 post	0.000	0.150	0.0	0	0	2.4e+04	15.84	15.07	6.6	2.05	0.54	0.66	4.6	0.9	0.10	2.6e+02	6.5
38	30	10 post	0.000	0.200	0.0	0	0	7.6e+04	15.05	16.07	6.5	1.94	0.76	0.56	8.7	15	0.12	6.2e+03	9
39	40	10 post	0.000	0.040	0.0	0	0	51	20.51	20.51	8.0	0.56	0.62	0.27	0.4	0.092	0.08	5.4	16
40	40	10 post	0.000	0.100	0.0	0	0	2e+03	20.48	20.10	7.6	1.46	0.77	0.68	0.36	0.75	0.03	79	7.5
41	40	10 post	0.000	0.200	0.0	0	0	3.2e+04	20.11	20.68	7.9	0.90	0.95	0.54	9.1	9.3	0.06	7.3e+03	9.9

Note. Columns (1)–(8) list the initial conditions for each simulation, including the total BBH mass, stellar model, pericenter distance, BBH semimajor axis, BBH eccentricity, BBH true anomaly, and BBH inclination. Columns (9) and (19) report the inspiral timescale of the BBH before and after the collision, respectively. Columns (10)–(15) describe the final outcomes, including the final masses of the BHs, the mass of ejected material, the envelope mass surrounding the BBH, the ratio of the final to initial BBH semimajor axis, and the final BBH eccentricity. Column (16) lists the amount of time between the disruption of the star and pericenter passage of the BHs. We define "disruption" to be when the mass bound to the more disruptive BH reaches 10% of the peak bound mass to that BH. Column (17) lists the amount of time the stellar debris is misaligned with the orbit (instantaneous χ_eff > 0) before becoming aligned (instantaneous χ_eff > 0). Column (18) lists the final effective spin parameter, calculated using Equation (1). Column (20) lists the kick velocity the center of the BBH receives due to mass ejection from the system.

Download table as:  ASCIITypeset image