Close encounters of black hole - star binaries with stellar-mass black holes

Dynamical interactions involving binaries play a crucial role in the evolution of star clusters and galaxies. We continue our investigation of the hydrodynamics of three-body encounters, focusing on binary black hole (BBH) formation, stellar disruption, and electromagnetic (EM) emission in dynamical interactions between a BH-star binary and a stellar-mass BH, using the moving-mesh hydrodynamics code {\small AREPO}. This type of encounters can be divided into two classes depending on whether the final outcome includes BBHs. This outcome is primarily determined by which two objects meet at the first closest approach. BBHs are more likely to form when the star and the incoming BH encounter first with an impact parameter smaller than the binary's semimajor axis. In this case, the star is frequently disrupted. On the other hand, when the two BHs encounter first, frequent consequences are an orbit perturbation of the original binary or a binary member exchange. For the parameters chosen in this study, BBH formation, accompanied by stellar disruption, happens in roughly 1 out of 4 encounters. The close correlation between BBH formation and stellar disruption has possible implications for EM counterparts at the binary's merger. The BH that disrupts the star is promptly surrounded by an optically and geometrically thick disk with accretion rates exceeding the Eddington limit. If the debris disk cools fast enough to become long-lived, EM counterparts can be produced at the time of the BBH merger.


INTRODUCTION
Dynamical interactions between stars and the compact objects they leave behind play an important role in dense environments, such as globular and nuclear star clusters and disks of Active Galactic Nuclei (AGNs). On global, large scales they can influence cluster thermodynamics (Hut et al. 1992), while on local, small-scales close interactions can alter the original birth composition of isolated stars and binaries.
Despite the relatively high fractions of stars and compact objects in binaries, hydrodynamic simulations of close encounters involving binaries have begun only recently. Lopez et al. (2019) and Ryu et al. (2022, Paper 1 in the following) studied close encounters between BBHs and single stars. They found that, in addition to altering the spin of the accreting BHs, tidal disruption events (TDEs) can have a significant impact on the binary BH's orbit, in ways which can be quantitatively different than the case of pure scattering. The EM signatures produced by these close encounters can also differ significantly from those of TDEs by isolated BHs: depending on the geometry of the encounter, the accretion rate can display periodic modulations with the orbital period. Detections of such events can provide constraints on the formation of BBH mergers .
More recently, Ryu et al. (2023, Paper 2 in the following) performed the first investigation of close encounters between binary stars and single BHs. Their hydrodynamic simulations showed a variety of possible outcomes, from full disruptions of both stars, to a full disruption of one star and a partial disruption of the other, to dissociation into bound and unbound single stars. Among these cases of dissociation, interesting outcomes include the formation of a runaway star, and of a fast-moving BH that accretes the tidally disrupted debris of the other star. In other outcomes, the binary stars are dissociated, and one of the stars is exchanged with the intruding BH, resulting in the formation of an X-ray binary.
Here we extend the line of investigation begun in Paper 1 and Paper 2 by performing a suite of hydrodynamic simulations of nearly parabolic close encounters between BH-star binaries and single BHs. Similarly to Paper 2, we use the moving-mesh code AREPO (Springel 2010;Pakmor et al. 2016;Weinberger et al. 2020), whose quasi-Lagrangian approach to hydrodynamics is well adapted to the problem. Our study aims to elucidate how such encounters can lead to a variety of outcomes, including EM transients due to the disruption of the star and the formation (via member exchange) of tight BBHs surrounded by debris material, potentially leading to a situation in which the BBH merger could be accompanied by an EM counterpart.
Our paper is organized as follows: § 2 presents the estimate for the rate of this type of encounters in globular clusters. § 3 describes the details of the numerical simulations and the initial conditions. We present our simulation results in § 4, with particular emphasis on a classification of the outcomes and its dependence on encounter parameters. We discuss these results in the context of their possible EM counterparts in § 5, and we finally summarize our work and conclude in § 6.

ENCOUNTER RATE IN GLOBULAR CLUSTERS
We begin by estimating the encounter rate of three-body interactions between BH-star binaries and single stars in globular clusters. Following Paper 2, we first calculate the differential rate of a single BH encountering a BH-star binary as dR/d • Σ rel . Here, is the binary number density in the vicinity of the BH, b s , where b is the non-interacting star -BH binary fraction, b 10 −4 − 10 −5 Kremer et al. 2018), and s is the number density of stellar-mass objects near the center of the clusters. The variable rel represents the relative velocity between the binary and the BH while Σ is the encounter cross-section. For p = √︁ 2 ( • + ★ )/ p , we can write Σ π ( • + ★ ) p / 2 where is the velocity dispersion. Adopting our results that strong encounters occur when p < , and assuming that this relation applies to binaries of any size and mass ratio, we can approximate Σ π ( • + ★ ) −2 . Then, we find that dR/d • can be expressed as Assuming more than a few tens of single stellar-mass black holes exist in dense clusters at the present day Askar et al. 2018;Kremer et al. 2018), and 150 globular clusters in the Milky Way (Harris 2010), the rate of strong three-body encounters per Milky Way-like galaxy is, Two outcomes that can be produced in this type of encounters are EM transients due to disruption of the star and formation of BBHs. In particular, because BBHs would likely form in 25% of all encounters (see § 4.4), the rate R of BBH-forming events would be (10 −9 ) yr −1 . The rate for encounters involving massive stars would be relatively high, compared to that for low-mass stars before all the massive stars turn into compact objects. However, as discussed in §5.3, the total number of this type of encounters over the full cluster lifetime would be higher for less massive stars because of their longer lifetime and higher abundance. Note that b depends on cluster parameters such as the initial binary fraction and the cluster age , and calculating R requires a detailed modeling of cluster evolution, as well as the star formation history. Thus, for a more precise estimate of R a more careful consideration of cluster evolution history is required.

Numerical methods
Our numerical methods and setup are essentially the same as in Paper 2. We perform a suite of 3D hydrodynamic simulations of close encounters using the massively parallel gravity and magnetohydrodynamics moving-mesh code AREPO (Springel 2010;Pakmor et al. 2016;Weinberger et al. 2020), which combines advantages of the two conventional hydrodynamical schemes, the Eulerian finite-volume method and the Lagrangian smoothed particle method, such as shock capturing without introducing an artificial viscosity, low advection errors, an efficient tracking of supersonic flow, and an automatically adaptive adjustment of spatial resolution. We use the HELMHOLTZ equation of state (Timmes & Swesty 2000) which accounts for radiation pressure, assuming local thermodynamic equilibrium. We include 8 isotopes (n,p,4 He,12 C,14 N,16 O,20 Ne, 24 Mg, Pakmor et al. 2012). We follow the advection of the elements which are then used for the update of the thermodynamics quantities (e.g., pressure). We do not follow the nuclear reactions, which should be fine given the short duration of the simulations and the reaction rates expected for the temperatures and densities that occur in our simulations.

Stellar model
The initial state of the star was taken as an evolved main-sequence (MS) star computed using the stellar evolution code MESA (version r22.05.1) (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019Jermyn et al. 2023). The star has an initial mass of 10 M and a metallicity = 0.006, which is lower than solar and consistent with what is found for globular clusters (e.g., VandenBerg et al. 2013), whose high stellar density facilitates dynamical interactions. 1 . Convection is modeled according to the mixing length theory with a mixing length parameter of 1.5. We use the Ledoux (1947) criterion to determine the boundary of the convective regions and include an exponential overshoot prescription (Herwig 2000) with parameters = 0.014 and 0 = 0.004. We treat semiconvection as in Langer et al. (1983Langer et al. ( , 1985 assuming it is fully efficient. Wind mass loss is modeled with the prescription from Vink et al. (2001). We evolve the star until about halfway through the main sequence, which we choose as the time when the central hydrogen mass fraction drops to 0.3, at which point the star has developed a 2.6 M convective helium core with a radius of 0.8 R . The stellar radius is ★ = 5.4 R , and the central density is 11 g cm −3 . Since we evolve the stellar model as a single star, we neglect possible past interactions that could have affected the structure. For example, if the black hole binary system formed through binary evolution, the star may have accreted mass (e.g., Renzo & Götberg 2021). If the system formed through dynamical capture, the star may have lost mass. We expect such effects to be small and not affect our results significantly.
We use the MESA stellar model as initial state for the AREPO simulation. After mapping the 1D MESA model into a 3D AREPO grid with 5 × 10 5 cells, the 3D single star is first relaxed. It usually takes up to five stellar dynamical time until it is fully relaxed.
The stellar dynamical time is defined as √︃ 3 ★ / ★ where ★ and ★ are the radius and mass of the star, respectively. Note that we increase the resolution for each single star by almost a factor of 2 compared to that in Paper 2, showing the results were converged with 2.5 × 10 5 . This is to more conservatively guarantee the convergence of our results, and to better resolve stars that may be partially disrupted during encounters. In addition to that, the resolution becomes finer over time by adopting refinement near the BHs ( § 3.4): some simulations with violent interactions have 10 7 cells at the end of the simulations. The density profiles of the relaxed stars considered in our simulations are depicted in Figure 1. As shown in the figure, the relative difference of our 3D star with respect to the MESA model is 1% for the inner region up to 2 M . The match is better than 10% throughout most of the star, except for the surface. We expect that this is sufficient for the aims of this study.

Black holes
As in Paper 2, we model the BH using a sink particle assuming it is not rotating initially. It only interacts gravitationally with gas and grows in mass via accretion of gas. We set the gravitational softening length of the BH ( 0.01 R ) to be ten times the minimum softening length of the cells of the stars.
We follow the same procedure for accretion described in Paper 2.
1 Note that the exact value of the metallicity does not affect our main results because the two most often cases in our simulations are fly-bys of stars around the BHs, or near-collisional disruptions of stars, and the outcomes of these two cases are not significantly affected by a slight change in the stellar internal structure due to a different metallicity.
However, we significantly improve the resolution near the BH using refinement (see § 3.4), which leads to more accurate estimates of the accretion rate with stricter conditions for accretion than in Paper 2. We search for cells bound to the BH (i.e. negative orbital energy relative to the BH) within 10 3 g (c.f., 1.5 × 10 4 g in Paper 2) where g = • / 2 is the gravitational radius of the BH and • denotes the mass of the BH. We still apply the same inverse-distance kernel (Monaghan & Lattanzio 1985) to put more weight onto closer cells. Although the change in the momentum and the mass of the BH due to accretion is taken into account, our simulations do not include potential radiative feedback produced by accretion.

Mesh refinement
The simulation code can adjust the local mesh resolution by adaptively splitting or merging cells if certain prescribed refinement criteria are satisfied (for more details, see § 6 in Springel 2010). We apply the refinement technique to cells in the vicinity of each BH to better resolve the stream structure there. At every time step, the code refines cells near a BH if all of the following conditions are met: (i) the distance from the BH fulfils < 5000 g , (ii) the cell density is > 2 × 10 −4 g cm −3 , (iii) the cell mass exceeds > 6 × 10 22 g, (iv) and Δ / > 0.26 for 500 g < < 5000 g and Δ > 500 g for < 500 g , where Δ is the cell size.
The refinement radius in condition (i) is chosen to be larger than the accretion radius (1000 g in this work) to ensure that gas streams inside the accretion radius are well resolved. Condition (ii) is designed to apply the refinement only to the cells that represent "real" gas, not vacuum regions. Criterion (iii) avoids a runaway creation of low-mass cells. Finally, the resolution limit imposed through condition (iv) guarantees that there are at least (10 2 ) cells within the accretion radius. On the other hand, at every time step, the code can also derefine cells within < 5000 g around each BH if the cell mass is < 1.5 × 10 22 g, meaning the mass of cells within < 5000 g never becomes smaller than this mass resolution limit.
We ran a few simulations with five different resolution limits within the range 0.05 ≤ Δ ≤ 0.4. This confirmed that the global evolution of the systems, such as their final interaction outcomes, is not affected by the refinement. The accretion rate has converged when the cell size fulfils Δ / < 0.3. Note that the number of cells within a volume at distance increases approximately by a factor of 8 when Δ decreases by a factor of 2.

Star -black hole binary
Before we carry out our encounter experiments, we relax binaries consisting of a fully relaxed star and a BH for 10 stellar dynamical times. We parameterize the binary's semimajor axis using an approximate analytic estimate of the Roche lobe radius (Eggleton 1983), where RL is the volume averaged Roche lobe radius of the star, = ★ / • is the mass ratio, and is the orbital separation. We define RL ≡ ( RL = ★ ) as the separation at which the star fills its Roche lobe. For = 0.5 and RL = ★ , RL 3.12 ★ 16.9 R .
We have performed this binary relaxation process for every binary with different orbital parameters (3 different binaries in total). The semimajor axis and the eccentricity of the relaxed binaries differ by less than 1% from their initial values.

Initial conditions
Following the same terminology as in Paper 2, we refer to quantities with a subscript containing − • as those relating to the orbit between a binary and a single BH. We assume a parabolic encounter with eccentricity 1− −• = 10 −5 between a single 10 M BH with a binary consisting of a 20 M BH and a 10 M star. The exact choices of the system parameters are somewhat arbitrary, but BHs with such masses have been found in X-ray binaries (Binder et al. 2021, e.g.). Encounters between objects of comparable masses are expected in the dense centers of young mass-segregated star clusters. We later discuss potential effects of different masses and orbits of encountering objects in § 5.3, based on our simulation results. We consider three semi-major axes for the initial binary systems: / RL = 2, 4 and 6, corresponding to an orbital period of 4, 12, and 22 days, respectively. We assume the binaries are circular at the start of our simulations. This is primarily to simplify the initial conditions, but this may not be unreasonable given that close binaries are often found to be circular (Almeida et al. 2017).
The distance between the binary's center of mass and the BH at the first closest approach p,b−• is parameterized using the impact parameter , i.e., p,b−• = 0.5 where is the binary semimajor axis. We consider = 1/4, 1/2, 1, and 2 for / RL = 4, and 1/2 for / RL = 2 and 6. The binary's angular momentum direction is always along the -axis in our simulations. We illustrate the initial configuration of the stellar binary and the BH in Figure 2.
We investigate the dependence of encounter outcomes on key encounter parameters, that is inclination angle = 0, 30 • , 60 • , 120 • and 180 • , = 1/4, 1/2, 1 and 2, and the phase angle = 0 • − 180 • with Δ = 45 • . We define as the initial angle between the line connecting the two members in the binary and the coordinate −axis (see Figure 2). We start by studying the dependence on the two phase angles of the binary ( = 0 • and 180 • ) by fixing all the other parameters. To achieve this, we initially rotate the binary while the initial separation between the center of mass of the binary and the BH is fixed at = 5 . This allows us to examine the outcomes from the first contact of the single BH with a different member of the binary. However, given the relatively high computational costs, instead of simulating encounters with every combination of and , we perform simulations for the encounters of the intermediate-size binaries ( / RL = 4) with different combinations of = 1/4, 1/2, 1 and 2, and = 30 • , 150 • , and = 0 • and 180 • . For the smallest and largest binaries ( / RL = 2 and 6), we only consider = 30 • and 150 • while = 1/2. In addition, we further examine the dependence of on the outcome properties by considering = 0, 60 • , 120 • and 180 • (for = 1/2). Last, we also study the impact of the phase angle on the encounter outcomes by simulating encounters with six additional phase angles ( = 45 • , 90 • , 135 • , 225 • , 270 • , and 315 • ).
In Table 1, we summarize the initial parameters considered in our simulations. Each of the models is integrated in time up to a few 100 p as needed to identify the final outcomes. Here, p = √︃ 3 p / is the dynamical time at = p , where is the total mass of the three objects (40 M ). The value of p for each model is given in Table 1.
The total computational cost for each run varies, mainly depending on how long the interactions last until the final outcomes are produced. Using 200 -300 CPU-cores of the Intel Xeon CascadeLake-AP processor (Xeon Platinum 9242), the total compute time per run has been around 70000 -100000 core hours.  (4), the numerical values encode (1) the initial semimajor axis of the binary / RL , (2) the impact parameter , (3) the initial phase angle in degrees, and (4) the initial inclination angle in degrees. Here, RL 3.12 ★ 16.9 R is the separation when the star in the binary fills its Roche lobe. The last three columns show the dynamical time p at pericenter (see its definition in the text), in units of hours, the orbital period of the binary, and the relative velocity orb of the binary members.

Classification of outcomes
The outcomes of three-body encounters between BH-star binaries and single BHs can be divided into three classes, depending on the final products.
(i) BBH-forming encounters: This class refers to encounters in which a BBH emerges. In this case, the impact parameter is mostly 1/2 − 1. The incoming single BH frequently interacts first with the star by the time the binary's center of mass and the single BH arrive at pericenter (models with "Yes" in the fifth column in Table 2). In this situation, the star in the binary nearly collides with the incoming single BH. We show one example for this type of encounter in Figure 3. The incoming single BH loses a significant amount of its kinetic energy and is gravitationally captured by the other BH initially in the binary. Because of the member exchange due to a violent star-removing encounter, the size of the final binary is not necessarily correlated with the size of the initial binary. To illustrate this, we compare in the top-left panel of Figure 4 the final of the BBHs with the semimajor axis of the initial BH-star binaries. The final covers over a wide range of values and is not necessarily comparable to of the initial binary. These violent interactions can lead to the formation of merging BBHs, as illustrated in the top-right panels: the GW-driven merger time scale of 5 out of 14 final BBHs is less than a Hubble time. Note that the absolute magnitude of the binding energy of the merging BBHs is much larger than the typical kinetic energy of stars in both globular and nuclear clusters ( 2 where is the velocity dispersion). This suggests that subsequent interactions with other stars would not dissociate these "hard" binaries, but rather make them more compact (Heggie 1975) and more eccentric (Valtonen & Karttunen 2006), which would facilitate their mergers. In addition, the disruption of the star prior to the BBH formation means that at least one member of the BBH is frequently surrounded by gas upon binary formation. When the BBH is compact, both BHs accrete gas.
(ii) Non-BBH-forming encounters: In this class, the outcomes are member exchanges between the two BHs or perturbations of the initial binary's orbit (models with "No" in the fifth column of Table 2). This mostly occurs when the two BHs interact at the first contact The color bar gives the logarithmic density in g cm −3 . The time is measured since the expected pericenter passage between the binary's center of mass and the single BH. At / p −16 (top-1 st ), the binary (star -green dot) and the single BH (yellow dot) approach each other. At / p = −3.21 (top-3 rd ), the incoming BH strongly encounters with the star in the binary, followed by the collision of the star (top-4 th ). The BH that disrupted the star is gravitationally captured by the other BH, forming a merging binary with 6.70 R and 0.943, corresponding to GW 10 4 yr (bottom panels).
between the binary and the single BH. We show one example for this type of encounter in Figure 5, resulting in an orbit perturbation. The impact of the encounters is relatively weak compared to the BBHforming encounters. As shown in the bottom-left panel of Figure 4, the final value of scatters within less than a factor of 2 around the initial . The eccentricity of the final BH-star binary is widely distributed between 0.1 and 0.9 (bottom-right panel), similarly to those of final BBHs (top-right panel). In this type of encounters, EM transient phenomena, such as tidal disruption events, collisions, or interacting binaries, can be created (e.g., Model 17. 2 1/2 0 30). In addition, the single BHs are ejected at 60 km s −1 , comparable to or higher than the escape velocity of globular clusters (i.e., tens of km s −1 ; Gnedin et al. 2002;. (iii) Undetermined: This class refers to cases where final outcomes are not determined (12 models in total, Models with "-" in the fifth column in Table 2 and with superscript ★ or ★★). Among these 12 models, there are eight encounters (models designated with superscript ★) in which the three objects form an unstable hierarchical triple, which we define as a triple where the outer binary is on a very large eccentric orbit so that the pericenter distance of the outer binary is smaller than the semimajor axis of the inner binary. In the table, we provide the orbital parameters of the inner binary. In the rest (models with superscript ★★), interactions become extremely prolonged so that a final outcome has not (yet) emerged.
From now on, we will focus on the first two classes, i.e., BBHforming and Non-BBH-forming encounters. These types of final outcomes and their properties are summarized in Table 2.

Dynamical processes
The two most crucial factors to determine the outcomes for the parameter space considered in our simulations are 1) the types of objects that meet at the closest encounter (BH -BH or BH -star), and 2) the net direction of the momentum kick relative to the bystander object (i.e., the object in the binary that does not interact with the incoming BH at the first closest approach) imparted by the interaction between the two meeting objects. As explained in the previous section, the first aspect substantially affects the chances of the survival of the star. The second aspect determines which objects end up in the final binary (i.e., member exchange, binary perturbation) and how the final binary's orbit looks like.
For the non-BBH-forming encounters with 1/2 − 1, the most frequent outcomes are either a member exchange or a perturbation of the original binary. The latter happens in retrograde encounters. This case can be categorized into three configurations, which are illustrated in Figure 6. In the first configuration (top panels), the incoming BH strongly interacts with the other BH and turns around at a small pericenter distance compared to the binary semimajor axis. The initially single BH is rapidly ejected from the system in the direction roughly opposite to the incoming direction. This interaction imparts a momentum kick to the BH that perturbs the binary orbit. In the other two configurations, the incoming BH either moves around or passes through the binary, without significant interactions with any of the binary members (bottom panels).
The dominant channel for member exchange is depicted in Figure 7. For the non-BBH-forming encounters in a prograde orbit, the two BHs meet first and pass through their points of closest approach. Like the first configuration for orbit perturbation, their relative motion gives a momentum kick to the motion of the BH originally in the binary, relative to the star. The momentum kick gives an additional acceleration in the BH's receding motion from the star. The initially single BH, after turning around the other BH, moves in a similar direction with the star and gravitationally captures it.
For the BBH-forming encounters, the star and the initially single  BH undergo close encounters, naturally resulting in a tidal disruption event or stellar collision. Both events can also impart a momentum kick to the disrupting BH. In our simulations, the momentum kick is not large enough to prevent the two BHs from forming a bound pair. For example, if the star and the incoming BH undergo a head-on collision, the incoming BH dramatically slows down and forms a merging BBH with the other BH (e.g., Model 20. 2 1/2 180 150). We caution that the head-on collision between the two equal-mass objects in our simulations is an extreme case yielding a dramatic drop in the BH's kinetic energy. The net effect of such star-removing events on the motion of the disrupting BH and the subsequent formation of a BBH depends on the mass ratio, relative velocity, and the direction of the momentum kick.

Binary black hole formation
Typical semimajor axes of BBHs formed in the BBH-forming encounters range within 10 − 400 R while eccentricities vary within 0.1 − 0.97. Correspondingly, the GW-driven merger time scales of those merging binaries are in the range 10 4 −10 13 yr. Five of our models among these encounters are merging BBHs with 7 − 150 R and 0.6 − 0.97. As explained in § 4.2, the dominant formation channel is the gravitational capture of the incoming BH by the BH originally in the binary after strong interactions between the incoming BH and the star. Naturally, a disruption event or a collision precedes the BBH formation. As a result, either the disrupting BH or both BHs accrete matter by the time they form a stable binary.

Dependence of outcomes on parameters
We examine the dependence of outcomes on a few key encounter parameters, phase angle , impact parameter , inclination angle , and semimajor axis , by varying one parameter at a time, keeping the rest of them fixed. Our simulations suggest that the two most important parameters that affect the formation of BBHs in this scenario of three-body encounters are the impact parameter and the phase angle.
(i) Phase angle : this is found to be one of the key parameters that separates BBH-forming encounters from non-BBH-forming encounters. For the former, very likely outcomes are BBHs, frequently accompanied by a disruption of the star. On the other hand, for the latter, frequent outcomes are eccentric BH-star binaries produced via member exchange or weak tidal perturbations of the initial stellar orbit. In addition, even for the BBH-forming encounters, the direction of the encounter between the initially isolated BH and the star at the first closest approach relative to the other BH determines the size of the semimajor axis of the BBH: if the momentum kick imparted on the encountering BH is given such that it adds to the encountering BH's momentum, a large binary forms (e.g., Model 6. 4 1 180 30 and Model 22. 6 1/2 180 30). Although our study is not appropriate for rate estimates, the dependence of outcomes on the phase angle may indicate that roughly 25% of these three-body encounters between objects of similar mass with 1 may possibly lead to BBH formation with a high chance of creating EM transients.
(ii) Impact parameter : in general, the initial binary and the single BH can interact significantly (member exchange or stellar collisions) at the first closest approach when p , which is also found in Ryu et al. (2022Ryu et al. ( , 2023. Fly-by only occurs at p > (Models 1. 4 2 0 30, and Model 9. 4 2 0 150). For this case, the initial binary orbit is weakly perturbed, resulting in a 10 -20% change in the semimajor axis. Relatively weak interactions also take place when the impact parameter is too small compared to the size of the binary, i.e., p < /8 (e.g., Model 4. 4 1/4 0 30, and Model 16. 4 1/4 180 150), as the single BH penetrates through the binary without interacting strongly with any of the binary members (see the bottom-right panel of Figure 6).
(iii) Inclination angle : prograde encounters tend to result in strong interactions between the first two encounter objects, frequently Figure 5. An example of a non-BBH-forming encounter, Model 17. 2 1/2 0 30. We depict the density distribution in the binary orbital plane at a few different times in units of p . At / p −16 (top-1 st ), the binary (star -green dot) and the single BH (yellow dot) approach each other. At / p = −3.21 (top-3 rd ), the two BHs encounter, resulting in the ejection of the initially single BH, while the initial binary orbit is significantly perturbed to become an interacting binary (bottom panels). Because of periodic interactions at pericenter, the binary orbit continues to evolve until the end of the simulation. The semimajor axis and eccentricity measured at the end of the simulation are 18 R and 0.4, respectively. leading to outcomes that involve member exchange (e.g., models with = 0 • and = 30 • ) or stellar collisions (e.g., models with = 180 • and = 150 • ). This is because the relative velocity between the two encountering objects is smaller, implying a larger gravitational focusing cross section (∝ 2 / 2 where is a typical relative velocity at infinity). A typical configuration for member exchange in prograde encounters is drawn in Figure 7. On the other hand, the first interactions in retrograde orbits are relatively weak due to the large relative speed between the two encountering objects. As a result, frequent outcomes are perturbations of the initial binary orbit, as depicted in Figure 6.
(iv) Semimajor axis : given the same pericenter distance relative to ( p 0.25 ) for the simulations with varying , the type of the final outcomes does not show a strong dependence on . However, the size of the final binary is closely correlated with that of the initial binary, e.g., 76 R of final binaries in models with / RL = 6 (or = 101 R ) and 61 R in models with / RL = 2 (or = 32 R ).

Accretion
Our simulations show that stars can be disrupted in three-body interactions between BH-star binaries and single BHs via strong interactions with very small impact parameters, i.e., collisions. In such events, a merging BBH can subsequently form, and at least one of the BHs is surrounded by an accretion disk which can create EM transient phenomena. To zero-th order, the disk structure and features of the accretion rate can be imprinted onto light curves of such events.
The refinement scheme adopted for the simulations allows us to resolve the gas structure down to 0.01 R 10 3 g from the BH. Although the regions that we can resolve are still too far from the BH to be directly related to the accretion process, we can provide an accurately resolved large scale structure of the disks formed in stardestroying events, which can be used as initial conditions for detailed disk simulations. Here, we define a disk as a group of gas cells tightly bound to the BH and coherently orbiting in the azimuthal direction. The outer edge of the disk is defined as the radius containing 99% of the total bound mass orbiting at a velocity exceeding 1% of the local Keplerian speed kep ( ) = √︁ [ (< ) + • ]/ , where (< ) is the mass enclosed within .
We show in Figure 8 both the face-on (left panels) and edge-on (right panels) density distributions of the disks around the BH that destroys the star at the first encounter in four example models, and in Figure 9 the radial profiles of the aspect ratio, the density, the temperature, and the rotational velocity for all models where an accretion disk forms. The aspect ratio ℎ/ is defined as the ratio of the first-moment density scale height, averaged over a given cylindrical radius, to the cylindrical radius. Here, we excluded Model 20.
2 1/2 180 150 in this analysis because the BH in that model is surrounded by a nearly spherical gas cloud, not by a disk. But we provide the accretion rate for that model also, shown in Figure 10.
We find that the disks are thick and pressure-supported, and mostly confined within 30 R . In general, the aspect ratio ℎ/ (top-left panel of Figure 9) is comparable to or greater than order unity up to the outer edge of the disks. ℎ/ declines from ℎ/ 3 − 5 to ℎ/ 1 outwards. The rotational velocity near the mid-plane is sub-Keplerian ( / kep 0.1 − 0.6), indicating the disk is not rotationally supported. The velocity ratio remains the same out to the outer disk edge. The density of the inner region stays flat at (0.1 − 5) g cm −3 up to 0.1 -0.2 of the disk size, then declines steeply following a −4 power-law. On the other hand, the temperature does not show such flatness at R , but continuously decreases following a −1 power-law.
Finally, we present in Figure 10 the accretion rate of the initially Figure 6. Schematic diagram showing three dominant configurations resulting in the perturbation of the original binary's orbit in our simulations. The black arrows indicate the direction of motion of objects, and long grey arrows the trajectory of the incoming BH (black circle). In the first configuration (top panels), for the prograde encounter with < 90 • and 1, the incoming BH undergoes a strong encounter with the BH (blue circle) initially in the binary (small closest approach distance compared to the binary semimajor axis), quickly turns around and advances to the left (top-right panel). This quick turn-around motion gives a momentum kick (green arrow) to the blue circle to the right with respect to the star (red circle). The orbit of the initial binary is perturbed. The second configuration (bottom-left panel) is a distant fly-by where the incoming BH does not significantly interact with any of the binary members, and this happens when 1. The last configuration (bottom-right panel) shows the case where the incoming BH passes through the binary without strong interactions with any of the binary members (e.g., 1/4 and / RL = 4).
single BHs that fully destroy the star at the first closest encounter. The general trend is that, upon disruption or collision, the accretion rate dramatically increases up to (10 −6 − 10 −5 ) M s −1 and it takes around 80-100 hours until declines by a factor of 100 from its peak. When the binary is eccentric and the pericenter distance is sufficiently close, a periodic perturbation from the other BH at periastron results in periodic bursts on a time scale the orbital period (e.g., Model 15. 4 1/2 180 150, and Model 20. 2 1/2 180 150). Although the accretion rate is super-Eddington, the total accreted mass is at most 0.1 M ( 0.4%) until the end of the simulation, and the magnitude of the BH spin driven by accretion can be as large as 0.01.
We have to caution that such extremely high accretion rates for stellar-mass BHs would result in strong outflows (e.g., Sądowski et al. 2014), which would regulate the accretion rate. Although we have  (Figure 6), the incoming BH strongly interacts with the BH (blue circle) initially in the binary, quickly turns around and advances to the right. This motion results in a momentum kick (green arrow) to the blue circle to the left with respect to the star (red circle). The initially single BH gravitationally captures the star and forms a binary. realized a significant improvement in resolving gas motions near the BHs compared to Paper 2 thanks to using refinement, since feedback from the BHs is not included in our simulations it is likely that our accretion rates are overestimated. Nonetheless, if the luminosity is mostly driven by accretion, the features revealed in the accretion rate (e.g., periodic bursts) could possibly be imprinted in the light curves.

Formation of merging binary black holes
Our simulations show that close three-body encounters between a BH-star binary and a single BH can create a merging BBH (see the top-left panel of Figure 4). One possibly dominant formation process we identified is the close interaction between the star and the incoming BH at the first closest approach, resulting in a stellar disruption, followed by the formation of a BBH. 5 out of 11 BBHs formed in our simulations would merge in a Hubble time via GW emission. The semimajor axes of the merging BBHs are 114 R and their eccentricities are quite high, 0.66 0.97. If the required conditions are met ( p 0.5 , encounters between the star and the incoming BH at the first closest approach), this type of encounters can form, albeit likely rarely, a very compact eccentric BBH: GW 10 4 yr in Model 20. 2 1/2 180 150.
To see whether the merging BBHs can have residual eccentricities when they enter the frequency band of LIGO (10 Hz to 10kHz), we evolve the five binaries assuming their orbits evolve purely via GW emission until GW = , where is the binary orbital period. We solve Equations 5.6 and 5.7 in Peters (1964) simultaneously using a 4th-order Runge-Kutta method with an adoptive step size of 10 −3 GW . As a sanity check, we confirmed that our numerical solutions are consistent with the analytic solution (Equation 5.11 in Peters 1964) within fractional errors of 10 −8 . Figure 11 shows the evolution of and of the five merging BBHs, starting from those found in our simulations (marked as scatters near the top-left corner). As shown in the figure, by the time the BBHs enter the LIGO frequency band, their residual eccentricities would be very small ( < 10 −5 ).
Nonetheless, the circumbinary gas produced by the disruption of the star may affect the (at least early) orbital evolution, which may hence deviate from the purely GW-driven evolution considered above. The gas-binary interaction and resulting binary evolution remains an active topic of study. A growing number of numerical works have suggested that a binary surrounded by a circumbinary disk can expand (e.g., Miranda et al. 2017;Muñoz et al. 2019;Duffell et al. 2020) and can be driven into an eccentric orbit (e.g., Zrake et al. 2021;D'Orazio & Duffell 2021), depending on the disk and binary parameters, as opposed to the predictions from the commonly held picture of surrounding gas driving binaries into shrinking circular orbits (e.g., Armitage & Natarajan 2002). However, given the limited parameter space explored in previous work, it is not straightforward to predict the evolution of our unequal-mass, very eccentric BBHs surrounded by a possibly misaligned disk, based on the results from the previous work.
The remaining 6 BBHs with GW-driven merger timescales longer than a Hubble time are hard binaries in typical stellar cluster environments. This means that those binaries could become potential GW event candidates via weak interactions with other objects and a few strong interactions like the ones considered in this study.

Electromagnetic counterparts of binary black hole merger
The close association of BBHs and stellar disruptions can have important implications for EM counterparts of BBH mergers. At the   time the BBH forms, there would be a prompt EM transient phenomenon due to the stellar disruption. The very high accretion rate (Figure 10), along with the accretion-driven BH spin and magnetic field of debris inherited from the star, suggests that a jet can be launched. For such a case, the luminosity powered by the jet would track the accretion rate as ∝ 2 with an uncertain efficiency factor. We also found that both BHs can be surrounded by the stellar debris and accrete, possibly suggesting that both BHs may be able to launch jets simultaneously, potentially leading to a unique observational signature.
In addition to the prompt EM emission, the existence of the surrounding gas when the BBH forms may result in a possible EM counterpart at the time of merger. This is a quite similar situation as found in Ryu et al. (2022) where an initially hard BBH encounters with a single star and becomes surrounded by gas debris after disrupting the star. Perna et al. (2016) studied the evolution of an initially hyper-Eddington accretion disk which cools and shuts down the magnetorotational instability before the disk material is fully accreted. Under these conditions, the "dead disk" is expected to survive until the BBHs merge, and to heat up and re-ignite during the merger process, hence yielding a possible EM counterpart to the GW event. Figure 11. The evolution of and of the five merging BBHs formed in three-body interactions due to GW emission. The markers depict and of the final merging binary black holes. The four grey horizontal lines indicate the semimajor axes at which the rest-frame GW frequency (twice the orbital frequency) is GW = 10 −4 , 10 −2 , 1, and 10 2 Hz, respectively.

Varieties of encounters
Although we consider three-body encounters between a circular binary and a single object with similar masses (the largest mass ratio is 0.5), there could be a variety of these types of events involving, e.g., initially eccentric binaries and a wide range of masses of encountering objects.
Encounters involving massive stars (i.e., 10 M ) are likely to occur during the early evolutionary stages of star clusters unless there is another episode of star formation, since stars with mass > 10 M would collapse to compact objects in tens of Myrs. Therefore, over the full cluster lifetime, the overall rate would indeed be higher for encounters involving less massive MS stars because such binaries would survive longer. Using Monte Carlo simulations of globular clusters, Kremer et al. (2018) showed that up to 10 detached BH-MS binaries can exist in clusters at an age of 10 − 12 Gyr, and the typical mass of the companion MS stars is 1 − 2 M , depending on the cluster properties. Even for this case, strong interactions between a low-mass MS star and the incoming BH at the first closest approach would have higher chances of forming BBHs than for the other cases where two BHs meet first.
At later times when all massive stars collapse to compact objects, interactions between stars and BHs with significantly different masses would be more probable. If the star is significantly less massive than both BHs, the interactions would be effectively two-body with small perturbations of the BH orbits by the star. However, if the star was disrupted by the incoming BH as in the BBH-forming encounters, the stellar disruption would generate bright EM flares. Furthermore, resulting momentum kicks and gas dynamical friction would facilitate the formation of BBHs, unless the momentum kick is given to increase the relative kinetic energy of the BHs. This process would be most efficient when the star and the incoming BH have comparable masses.
If the encountering binary is eccentric, the binary members would spend most of their orbital time near apocenter, indicating that the cross-section would be enhanced by a factor of 1 + . Unlike the increase in in our models when initially circular binaries are considered, the final binaries can be circularized depending on the direction of the momentum kick associated with close interactions between the two objects at the first closest approach. We already demonstrated in Figures 3 and 5 that the momentum kick acts to add to or remove the momentum of the BH in the binary, depending on whether the orbit is initially in a prograde or retrograde direction.

SUMMARY AND CONCLUSIONS
In this work we have investigated the outcomes of three-body encounters between a 20 M BH -10 M star circular binary and a 10 M stellar-mass BH, using a suite of hydrodynamical simulations with the moving-mesh code AREPO. We have focused on the formation of BBHs, the conditions required for their formation, and the EM emission from those systems. We have considered a wide range of encounter parameters, i.e., varying the binary size ( 34, 68, 101 R ), the impact parameter ( /4 − ), the inclination angle ( = 0 • , 30 • , 120 • , 150 • , and 180 • ), and the phase angle ( = 0 • , 45 • , 90 • and 135 • , 180 • , 225 • , 270 • , 315 • ), while we have kept fixed the masses of the star and of the BHs.
We have categorized the encounters into two classes depending on their outcomes. This classification is primarily determined by which types of objects meet at the first closest approach. When the star and the incoming single BH encounter first, their close interaction imparts a momentum kick to the BH, resulting in a dramatic decrease in the BH's speed. The BH is subsequently captured by the other bystander BH, forming a BBH. In this case, the star is frequently destroyed due to its close encounter with the BH. On the other hand, when two BHs encounter first, either the original binary's orbit is simply perturbed (prograde encounters), or the originally single BH captures the star, forming a new binary (member exchange, retrograde encounters). Although the most frequent outcomes are BH-star binaries, a disruption of the star and BBH formation are still possible.
The most important factors that determine the outcomes are the phase angle and the impact parameter. As explained above, the phase angle primarily demarcates the boundary between "BBH-forming" encounters and "non-BBH-forming" encounters. The impact parameter on the other hand affects the strength of interactions: for p > , the incoming BH interacts weakly with the binary. As a result, the binary orbit is perturbed, or the binary members are exchanged. For p , interactions can become significant, possibly resulting in a disruption of the star when the star and the BH meet at the first closest approach. Although our simulations do not cover the entire parameter space for this type of encounters, the key dynamical processes can be extrapolated within this class of encounters to other initial parameters, and possibly also to other astrophysical systems (e.g., three-body encounters involving a massive black hole having a stellar companion and an isolated BH, forming extreme mass ratio inspirals).
The close correlation between BBH formation and stellar disruption in our systems has interesting implications for the formation channel of BBHs and EM counterparts of their merger. We confirm that three-body encounters between a BH-star binary and a BH can produce merging BBHs. In addition, we find that the BH that disrupts the star in the BBH-forming encounters is promptly surrounded by an optically and geometrically thick disk with accretion flows towards the BH exceeding the Eddington limit. If a jet is launched from the system, the jet luminosity would likely track the accretion rate. If the disk remains long-lived and revives at merger, EM counterparts can be produced at the time of the BBH merger.
Our order-of-magnitude estimate for the encounter rate suggests that this type of encounters may be rarer than other types of three-body encounters considered in Paper 1 (i.e., between binary BHs and single stars) and Paper 2 (i.e., between stellar-binaries and single BHs). However, given the simplified assumptions made here, more detailed estimates should be made for these encounters, taking their specific astrophysical environments accurately into account.