Long-term Evolution of Tightly-Packed Stellar Black Holes in AGN Disks: Formation of Merging Black-Hole Binaries via Close Encounters

We study the long-term evolution of two or more stellar black holes (BHs) on initially separated but unstable circular orbits around a supermassive BH (SMBH). Such a close-packed orbital configuration can naturally arise from BH migrations in the AGN disk. Dynamical instability of the orbits leads to recurring close encounters between two BHs, during which the BH separation $r_{\rm p}$ becomes less than the Hill radius $R_{\rm H}$. In the rare very close encounters (with $r_{\rm p}$ several orders of magnitude less than $R_{\rm H}$), a tight merging BH binary can form with the help of gravitational wave emission. We use $N$-body simulations to study the time evolution of close encounters of various degrees of"closeness"and the property of the resulting binary BH mergers. For a typical"SMBH + 2 BHs"system, the averaged cumulative number of close encounters (with $r_{\rm p} \lesssim R_{\rm H}$) scales approximately as $\propto t^{0.5}$. The minimum encounter separation $r_{\rm p}$ follows a linear cumulative distribution $P(<r_{\rm p}) \propto r_{\rm p}$ for $r_{\rm p} \ll R_{\rm H}$. From these, we obtain a semi-analytical expression for the averaged rate of binary captures that lead to BH mergers. Our results suggest that close-packed BHs in AGN disks may take a long time ($\gtrsim 10^7$ orbits around the SMBH) to experience a sufficiently close encounter and form a bound binary, although this time can be shorter if the initial BH orbits are highly aligned. The BH binary mergers produced in this scenario always have high eccentricities when entering the LIGO band, and have a broad distribution of orbital inclinations relative to the original AGN disk. We also explore the effects of the gas disk and find that simple gas drags on the BHs do not necessarily lead to an enhanced BH binary capture rate.

The possibility of BBH mergers in AGN accretion disks around supermassive black holes (SMBHs) has received much attention in recent years (e.g., McKernan et al. 2012McKernan et al. , 2014Bartos et al. 2017;Stone et al. 2017;Secunda et al. 2019Secunda et al. , 2020Yang et al. 2019a,b;Gröbner et al. 2020;Tagawa et al. 2020a,b;Ford & McKernan 2021). BBHs in flat disks may be hardened, or even driven to merger, by a series of nearly co-planar binary-single scatterings (e.g., Stone et al. 2017;Leigh et al. 2018;Samsing et al. 2020). Hydrodynamical interaction between the gaseous AGN disk and a pre-existing BH binary may also influence the orbital evolution of the binary (Baruteau et al. 2011;Li et al. 2021a,b;Li & Lai 2022). BH Mergers that happen inside AGN disks could have several observable properties, such as the possible associations with electromagnetic counterparts (Stone et al. 2017;McKernan et al. 2019;Graham et al. 2020) and distinct mass and spin distributions (e.g., McKernan et al. 2018;Yang et al. 2019a).
The AGN disk channel of BBH mergers typically relies on the disks being an ideal environment for dynamical formation of BH binaries. BBHs could be "preexisting" in nuclear star clusters and get captured in the inner AGN disks (Bartos et al. 2017) or form in situ in the extended region of the disks (Stone et al. 2017). It was suggested that AGN disks (Sirko & Goodman 2003;Thompson et al. 2005) may contain migration traps where stellar-mass BHs (sBHs) can accumulate (Bellovary et al. 2016). Orbiters inside such traps can have close encounters with each other due to their mutual gravity and potentially form binaries (Secunda et al. 2019(Secunda et al. , 2020. A direct gas-capture channel has also been proposed to form binaries from single sBHs (Tagawa et al. 2020b).
An important and unsolved issue of the AGN disk BBH merger scenario is the evolution of sBHs during close encounters and how bound BBHs actually form. Previous studies tend to adopt very simple prescriptions for BBH formation and follow-up evolution. For example, Secunda et al. (2019) performed N -body simulations of multiple sBHs and included disk force prescriptions to mimic the effects of eccentricity damping and migration traps, but they assumed that all close encounters between two BHs within the mutual Hill radius and with negative relative binary energy lead to the formation of bound binaries that quickly merge. This assumption is problematic as the vast majority of such bound binaries are short-lived and are quickly destroyed by the tidal force of the SMBH. Tagawa et al. (2020b) considered gas-assisted binary formation in their population synthesis study, including the time delay between BBH formation and merger, and allowing newly formed BBHs to be disrupted before mergers. However, their studies are one-dimensional and they assumed that the relative orbits of all BBHs are circular for simplicity.
In this paper, we study how often two or more sBHs in closely-packed, initially circular and nearly co-planar orbits around a SMBH can be captured into a "perma-nent" binary and merge with the aid of the gravitational radiation. Such close-packed orbits (with the difference in orbital radii less than a few times the Hill radius R H ) can be naturally produced by the differential migrations of sBHs in the AGN disk and/or the migration traps (Bellovary et al. 2016;Secunda et al. 2019). Since the sBHs in close orbits are dynamically unstable, they exhibit chaotic orbital motion around the SBMH and undergo repeated close encounters with each other (with their separation less than R H ). Binary capture occurs only in the very rare occasion when two sBHs experience an extreme close encounter (with their separation several orders of magnitude less than R H ), during which energy dissipation through GW emission is effective. We perform a large number of N -body simulations to study the occurrence rate of such close encounters and the properties of the captured BBHs. Since the dynamical influence of the disk gas on sBHs during close encounters is highly uncertain, the major part of our paper focuses on the clean problem where the only dissipative effect is GW emission. Nevertheless, we also carry out an exploratory study on the gas effects by adding "frictional" forces on the sBHs to mimic the BH-gas interactions.
The rest of this paper is structured as follows. In Section 2, we introduce our scenario for BBH formation in AGN disks via close encounters between sBHs. In Section 3, we describe our fiducial "SMBH + 2 BHs" simulations with no gas effects included. We use these simulations to obtain the time evolution of the close encounter rate, the distribution of the minimum sBH separations during encounters, and the timescale (and the probability) to form long-lived or merging BH binaries. Section 4 discusses how our results depend on the initial inclinations of the BH orbits around the SMBH. In Section 5 we apply simple models of the disk forces on the sBHs to assess the effects of gas disks on the evolution of the embedded sBHs and the formation rate of merging BBHs. Section 6 explores the rate and properties of the close encounters in systems with more than two sBHs around the SMBH. In Section 7, we summarize our findings.

SCENARIO
AGN disks can help bringing stellar-mass black holes (sBHs) circulating around a supermassive BH (SMBH) into close orbits due to the differential migrations of the BHs and migration traps (Bellovary et al. 2016), therefore promoting close encounters between the sBHs. While the encountering sBHs typically have too much relative kinetic energy to become bound in the presence of the tidal field of the SMBH, they may occasionally have a very close encounter, during which gravitational wave (GW) emission can take away the excessive energy. Dynamical friction from the disk gas may also play a role, but its effect is more difficult to quantify. These very close encounters may turn the two BHs into a bound binary and lead to BBH merger. We consider a system with a central SMBH of mass M , around which orbit two or more sBHs on nearly circular and nearly co-planar trajectories. For simplicity, henceforth 'BHs' always refer to stellar-mass black holes that are orbiting around the SMBH. Due to their migrations in the AGN disk, the BHs may have orbits very close to each other. For the most part of this paper, we set up two BHs with masses m 1 , m 2 and initial semi-major axes a 1 , a 2 around the SMBH. If the dimensionless orbital separation (in units of the mutual Hill radius R H12 ) is less than a critial value, i.e. where the orbits are dynamically unstable, such that the two BHs will soon develop orbital crossing and start chaotic evolution. The boundary between "stable" and "unstable" can be fuzzy but the critical value K crit is of order unity and depends on the "frictional" force acting on the BH from the disk gas (Li, Rodet and Lai 2022).
In the absence of BH-disk interaction, the Hill stability criterion gives K crit = 2 √ 3 (Gladman 1993). There are three possible outcomes for the two BHs in unstable orbits: (i) The lighter BH (m 2 ) is ejected from the system; (ii) The two BHs experience a sufficiently close encounter such that GW emission and/or gas drag makes them into a bound binary and eventually merge; (iii) One of the BHs moves very close to the SMBH and gets "swallowed" by it. Outcome (iii) has a negligible probability when a 1 , a 2 GM/c 2 , i.e. when the horizon radius of the SMBH is much less than the BH orbital distances. In our "SMBH + 2 BHs" systems with M m 1 , m 2 , outcome (i) will take many orbital periods to happen. This can be understood as follows. Close encounters between m 1 and m 2 (< m 1 ) cause m 2 to experience energy diffusion, with the change (loss or gain) of energy during each encounter given by ∆E ∼ α(Gm 1 m 2 /a 1 ), where α 1. Thus the average number of close encounters required for m 2 to be ejected is N ej ∼ (GM m 2 /2a 1 ) 2 (∆E) −2 ∼ (4α 2 ) −1 (M/m 1 ) 2 . Indeed, extensive numerical experiments carried out by Pu & Lai (2021) show that N ej has a broad distribution, with the mean value given by (see their Eq 24) and the distribution has a long tail at the larger values (the 68% quantile of N ej ranges from 0.25 N ej to 13 N ej ). The ejection time is usually much longer then N ej P 2 (where P 2 is the initial orbital period of m 2 ) since the semi-major axis of m 2 increases as it approaches ejection. Thus, for M/m 1 10 6 , the ejection timescale t ej 10 10 P 2 , which means ejection almost never happens.
We are thus left with outcome (ii), i.e. BH binary formation due to dissipative processes. If we neglect the possible effect of gas drag, the only dissipation is GW emission ("gravitational bremsstrahlung"). During a very close encounter (i.e. the separation between m 1 and m 2 is much less than R H ), the two BHs lose their relative energy by the amount (Peters 1964;Quinlan & Shapiro 1989) where m 12 = m 1 + m 2 , µ = m 1 m 2 /m 12 , and r p is the periastron separation of the relative trajectory of the BHs. To form a binary that is stable in the presence of the SMBH tidal field, we require with η of order unity. This implies that the pericenter distance of the relative orbit of m 1 and m 2 must be less than the critical value for capture, given by , where a 12 = (a 1 + a 2 )/2 a 1 . Thus, for the two unstable BHs to be captured into a bound binary requires an extremely close encounter with r p 10 −4 R H12 . A major goal of our paper is to evaluate the rate of such extremely close encounters and the properties of the resulted merging BH binaries.

FORMATION OF BOUND BINARIES IN "SMBH + 2 BHS" SYSTEMS
In this section, we study the formation rate of bound binaries for two BHs in unstable orbits around a SMBH. We neglect the effects of gas disk (if any) here -these effects will be studied in Section 5.

Setup of simulations
The fiducial system of our simulations consists of a central SMBH with mass M and two smaller BHs with mass m 1 = 2 × 10 −5 M and m 2 = 10 −5 M . The initial orbital separation is set as so that their orbits are unstable. For convenience, henceforth we define to be the Hill radius of the m 1 at the beginning of the simulation and use it as the unit of distance. We let P 1 be the initial orbital period of m 1 and use it as the unit of time. We note that P 1 at a 1 = 100GM/c 2 for M = 10 6 M is 10 −3 yr. The BHs are given initial eccentricities e 1 = 0, e 2 = 10 −5 , and inclinations i 1 = i 2 = R H /a 1 . We carry out 2000 runs and sample the initial values of the argument of the peripasis, the longitude of the ascending node, and the mean anomaly randomly in the range [0, 2π] for each BH, assuming they all have uniform distributions. We simulate the evolution of the system using the Nbody code REBOUND (Rein & Liu 2012) with the IAS15 integrator (Rein & Spiegel 2015). Each simulation runs for 10 5 P 1 .

Close encounters
To characterize encounters with various degrees of "closeness", we designate a close encounter event to be "CEα" when the separation between the BHs, r rel = |r 1 − r 2 |, becomes less than 10 −α R H and numerically keep track of CE0, CE1 and CE2 in each simulation. A CEα event starts when r rel changes from > 10 −α R H to < 10 −α R H , and ends when r rel becomes greater than 10 −α R H again and the relative energy The whole process is regarded as a single CEα event no matter how long it elapses. In the simulation, a new CEα is logged only if the previous one has ended.
The left panels of Figure  other. The right panel of Figure 1 shows the probability distribution of N α (t = 10 5 P 1 ) from all 2000 simulations. The total numbers of CE0, CE1, CE2 events in those simulations are 248790, 43153, 5238, respectively. Thus, at t = 10 5 P 1 , the average number per simulation of CE0, CE1, and CE2 are N α (10 5 P 1 ) = 124, 22, and 2.6, respectively. Note that N α has a wide distribution: For example, while N 2 (10 5 P 1 ) = 2.6, 5% of the runs have N 2 (10 5 P 1 ) > 10. Despite the difference in the numbers, the three kinds of CEs all have higher occurrence rate at the early times than at the later times. The time evolution of the average can be described by a power law N α (t) ∝ t nα . We perform least-square fits using such power-law to the results at t 2 × 10 4 P 1 to exclude the initial transient stage. We find

Very close encounters and BBH formation rate
During each CEα event, we take a simulation snapshot every 10 −3 P 1 and monitor the separation r rel = |r 1 −r 2 | between the two BHs. The exact periapsis passage moment may lie between two of the snapshots. We use the snapshot right after the true pericenter passage (when r rel first increases) to calculate the relative energy and angular momentum of the encountering BHs: We then compute the semi-major axis a rel of the orbit via E rel = −Gm 1 m 2 /(2a rel ), and the pericenter distance r p using rel = Gm 12 a rel (1 − e 2 rel ) 2Gm 12 r p .  Figure 2 shows the distribution of the semi-major axis a rel of the BH binary undergoing close encounters. The most likely a rel is about R H , and nearly all close encounters have R H /a rel ≤ 2.3. Given that the BH pairs during the encounters are either unbound or have a rel R H , they quickly depart from each other after reaching their minimum separation and get disrupted by the SMBH tide. This implies that the lifetime of most binaries formed by CEs is less than their orbital period, which is approximately equal to their period around the SMBH when a rel ∼ R H . Indeed, we find that, among the 248790 CE0 events in our 2000 simulations, the vast majority (96.35%) disentangle within one orbital period (P 1 ) around the SMBH, and only 0.02% survive for more than 10P 1 .  (14). Figure 3 shows the distribution of the pericenter ("closest") distance of the two BHs during close encounters. We find that, for a CEα event, the cumulative distribution of r p (i.e., the probability for the pericenter distance to be less than r p ) is approximately given by (14) Equation (14) becomes increasingly accurate for closer encounters. The difference in N α (t) for CE0, CE1 and CE2 discussed in Section 3.2 is a direct consequence of this cumulative distribution of r p . Equation (14) is equivalent to a constant probability distribution dP α /dr p = const. Since the relative angular momentum is rel 2Gm 12 r p for nearly parabolic encounters (for r p a rel ∼ R H ; see Figure 2), this is equivalent to the probability distribution in rel given by This relation has been previously obtained both analytically and numerically in the context of planetary collisions (Li & Lai 2020). It can be understood as follows: When the two BHs reach a close separation r rel R H , the angular momentum is rel = v rel r ⊥ , where v rel 2Gm 12 /r rel and r ⊥ is the projection of r rel perpendicular to v rel ; if the initial mutual inclination of the two BHs is much larger than r rel /a 1 , the vector r ⊥ is sampled uniformly in the plane perpendicular to v rel (Li & Lai 2020), i.e. dP/(r ⊥ dr ⊥ ) = const, which implies Equation (15).
While Figure 3 depicts the r p -distribution during the whole simulation time span, Figure 4 shows that the CEs from different time intervals also follow to the same distribution.  Thus, combining Equations (9) and (14), the averaged cumulative number of very close encounters with r p < r cap given by Evaluating Equation (16) using N 2 (t) and P 2 (< r cap ), we find where we have scaled r cap according to the estimate given by Equation (6). Thus, for r cap 10 −4 R H , it would take about 10 8 P 1 on average for the two BHs to capture into a merging binary by GW emission. Note that, because of the wide spread of N (t) (see Figure 1), it is possible for an unstable pair of BHs to reach N (t) 4 N (t) (about 5% of the systems) and be captured in ∼ 10 7 P 1 .

Orbital obliquities of BH binaries formed in close encounters
The orbital obliquities (i bin ) of the BH binaries formed by CEs can be inferred from the inclination of the angular momentum vector rel of the CE orbits with respect to the initial orbital axis (ẑ) around the SMBH, . Probability density function of the binary obliquity i bin of BHs in close encounters. Note that i bin measures the angle between the close encounter plane and the orbital plane around the SMBH, i.e. cos i bin ≡ˆ rel ·ẑ, whereˆ rel is the relative angular momentum axis of the two BHs, andẑ is the initial orbital axis around the SMBH.
i.e. cos i bin =ˆ rel ·ẑ. Figure 5 shows the distribution of cos i bin for the close encounters in our fiducial simulations. Encounters within Hill radius (CE0, black) are predominantly prograde with the most probable i bin being zero. But closer encounters (CE1 and CE2) have nearly uniform distribution of cos i bin . This indicates that merging BH binaries formed in very close encounters (r p R H ) have a broad range of inclinations with respect to the AGN disk plane.
The broad distribution of cos i bin is consistent with the findings of Li & Lai (2020) in the context of planet collisions. Such a broad distribution arises when the initial mutual inclination of the BH orbits ∆i is much larger than r p /a 1 . Since our simulations have ∆i ∼ R H /a 1 , all CEs with r p R H will have a broad cos i bin distribution. 1

Mass dependence of the close encounter rate
The mass ratio between the SMBH and the BHs in AGN disks can span several orders of magnitude. We repeat our fiducial experiment (which has m 1 = 2m 2 = 2 × 10 −5 M ) with four different sets of BH masses: (m 1 , m 2 ) = (8, 4) × 10 −5 M , (m 1 , m 2 ) = (2, 1) × 10 −6 M , (m 1 , m 2 ) = (2, 1.5) × 10 −5 M , and (m 1 , m 2 ) = (2, 0.5) × 10 −5 M . We use them to test the effect of having larger 1 Li & Lai (2020) showed that when ∆i rp/a 1 and ∆i R H /a 1 , the distribution of cos i bin is f (cos i bin ) ∝ (1−cos 2 i bin ) −1/2 , and that when ∆i rp/a 1 and ∆i ∼ R H /a 1 , f (cos i bin ) becomes more uniform. . The texts on the right panels show the average number of N2 at t = 10 5 P1 and the percentage of the runs that experience more than 10 CE2 events.
BH mass, smaller BH mass, smaller BH mass ratio (m 1 /m 2 = 4/3), and larger BH mass ratio (m 1 /m 2 = 4), respectively. Figure 6 compares the CE2 event rates in the fiducial system and in the four systems with different BH masses. In all cases, the N 2 -distributions at t = 10 5 P 1 show large spreads, with the mean N 2 (t = 10 5 P 1 ) in the range of 1.6-2.9 (see Figure 6). The probabilities for high N 2 s (beyond the mean) are more different: for example, the probability for each simulation to get more than 10 CE2 events is 5%, 9%, 1%, 4%, and 6% for the five cases displayed in Figure 6.
We find that, for four systems with different BH masses, the probability distributions of R H /a rel , r p /R H and i bin are all similar to the fiducial results as depicted in Figures 2, 3 and 5.
Overall, our results suggest that the BBH formation rate obtained for our fiducial system in Section 3.3 can be applied to other systems with different BH masses (as long as m 1 , m 2 M ) to within a factor of 2. The biggest effect that we observe is a drop in the CE2 rate by a factor of 2 when the BH masses are lowered by a factor of 10 compared to the fiducial system.

Results with different initial inclinations
Our simulations in Section 3 assume that the initial orbits of the two BHs have inclinations i 1 = i 2 = R H /a 1 and random longitudes of ascending nodes. We expect the CE rates to increase when i 1 and i 2 are smaller. Figure 7 shows the cumulative number of CE2 counts (averaged over 2000 runs) in simulations with different initial values of i 1 and i 2 . We see that, as expected, N 2 (t) increase as i 1 and i 2 decrease. This increase is highly substantial in the exact coplanar systems (i 1 = i 2 =0), for which we find an average of 56 CE2s at t = 10 5 P 1 and the CE2 rate evolve in time as This should be compared to Equation (11) for our fiducial runs (which assume i 1 = i 2 = R H /a 1 ). We notice that a smaller (but non-zero) initial inclination leads to more CE2 at the beginning, but the various N 2 (t) curves become roughly parallel after a few ten-thousand orbits, indicating that the rate d N 2 /dt "saturates" to the fiducial value (see below). ).
The distributions of the pericenter ("closest") distance r p of CE2 for the simulations with different initial inclinations are shown in Figure 8. For the exact co-planar systems, we find that the cumulative distribution of r p is P α (< r p ) r p 10 −α R H 1/2 (for r p ≤ 10 −α R H ). (19) This should be contrasted with Equation (14) for our fiducial simulations. Equation (19) is equivalent to a uniform distribution of angular momentum, dP/d rel = const. Such a uniform distribution is expected because for i 1 = i 2 = 0, the dynamics is confined to the original orbital plane, and at r rel R H the projected separation r ⊥ is restricted to a line with dP/dr ⊥ = const.
Combining Equations (18) and (19), we find that for the exact co-planar systems, the cumulative BH binary formation rate in very close encounters is It would only take about 6 × 10 3 P 1 on average for the two BHs to form a merging binary. Compared to Equation (17), we see that for such perfectly coplanar systems, the BH binary formation rate is much higher than our fiducial systems (with i 1 = i 2 = R H /a 1 ).
All of the other small initial inclination simulations yield between 0.1 and 1.0 critical CEs with r p < 10 −4 R H on average in their first 10 5 P 1 . We expect their rates to be similar to the fiducial rates over longer timescales (see below). Figure 9. Time evolution of the mutual inclination of the two BH orbits in the fiducial system (green) and in the nearly coplanar system with initial i1 = i2 = 10 −5 RH/a1 (blue). The colored lines represent 200 individual simulations for each system. The black curves represent their averages.

Evolution of nearly coplanar systems
Figures 7 and 8 show that for nearly coplanar systems (with initial inclinations 0 < i 1 , i 2 R H /a 1 ), both the cumulative CE2 rate and the r p distribution P α (< r p ) lie between those of the exact coplanar system and our fiducial system. However, Figure 7 also indicates that while the cumulative number of CE2 events in the first ∼ 10 4 orbits is much larger than the fiducial case, the rate d N 2 /dt seems to settle down to the fiducial rate at later times. This suggests that the mutual inclinations of the BH orbits grow in time in systems with nearly coplanar initial orbits. Figure 9 shows the evolution of mutual inclination θ 12 of the BH orbits for our fiducial systems and for a system with initial i 1 = i 2 = 10 −5 R H /a 1 . The mutual inclination is computed from cos θ 12 =ˆ 1 ·ˆ 2 , whereˆ 1 andˆ 2 are the unit angular momentum vectors of m 1 and m 2 around the SMBH. We see that in the nearly coplanar system, the mutual inclination θ 12 gradually increases and then saturates at θ 12 ∼ R H /a 1 . For our fiducial system, the average mutual inclination remains at ∼ R H /a throughout the simulation.
Due to the time evolution of the mutual inclination θ 12 , we expect that the CE statistics, such as the r pdistribution, also evolve with time in the nearly coplanar systems. Figure 10 shows the time dependence of the r p -distribution for the i 1 = i 2 = 10 −5 system. It is clear that the initial small inclinations (i 1 , i 2 R H /a 1 ) only affect the results at earlier times. The long-term statistics of CEs for nearly coplanar systems is better 10 4 10 2 10 0 r p /R H 10 4 10 3 10 2 10 1 10 0 CDF coplanar t = 10 2.0 3.5 P 1 t = 10 3.5 4.0 P 1 t = 10 4.0 4.5 P 1 t = 10 4.5 5.0 P 1 fiducial Figure 10. Same as Figure 4, showing the time evolution of the cumulative distribution of rp. The colored lines show the results for a system with initial inclination i1 = i2 = 10 −5 RH/a1, and different colors indicate different time intervals of the simulations. The black lines show the total distribution for the exact coplanar system (dash-dot) and the fiducial system (solid) for comparison.
represented by our fiducial simulations with i 1 = i 2 = R H /a 1 (see Section 3).

EFFECTS OF FRICTIONAL DISK FORCES
The AGN disk can affect the dynamical evolution of the embedded BHs. For example, a BH (with mass m 1 and semi-major axis a 1 ) experiences eccentricity damping on the timescale (Tanaka & Ward 2004) τ e M 2 h 4 2πm 1 Σa 2 1 P 1 1.2 × 10 6 a 1 10 2 GM/c 2 −2 Σ 10 5 g/cm where Σ is the disk surface density and h is the disk aspect ratio, and we have adopted some representative parameters for AGN disks (e.g., Sirko & Goodman 2003; see Figure 1 of Secunda et al. 2019). In the previous sections, we have studied BH binary captures via very close encounters in the absence of any disk force on the BHs. A full exploration of the effects of disk forces on BH binary formation would require long-term hydrodynamical simulations and is beyond the scope of this paper. Here, to qualitatively assess the effect of the disk, we apply simple prescriptions of disk forces on the BHs in our N -body simulations. We consider two simple models for the disk forces: 1. The first model includes the frictional force (per unit mass): where v K = GM/r 3θ is the Keplerian velocity and r is the instantaneous distance of the BH to the SMBH. This force tends to damp the BH velocity v to v K and damp its eccentricity around the SMBH at the rateė −e/τ drag .
2. The second model includes a force that mimics a migration trap at radius r 0 : where Ω K,0 = GM/r 3 0 is the Keplerian frequency at r 0 . Equation (23) assumes that the torque on the BH at r is approximately linear in (r − r 0 ) near the trap. In the following, we set r 0 to a 1 , the initial semi-major axis of m 1 .
The constants τ drag and τ trap in Equations (22) to (23) characterize the strengths of the disk forces. We apply these disk forces to our fiducial systems (ses Section 3), considering different values of τ drag and τ trap . For each value case, we perform 200 simulations up to 10 5 orbit.  Figure 6, but for the fiducial systems (see Figure 1) and including the effects of disk forces. The type and the strength (timescale) of the adopted disk force for each panel is as labeled. Figure 11 presents three example simulations: F drag only with τ drag = 10 6 P 1 , F trap only with τ trap = 10 6 P 1 , and both F drag and F trap with τ drag = τ trap = 10 6 P 1 . The left panel shows the time evolution of N 2 and the right panel shows distribution of N 2 at t = 10 5 P 1 . In all three cases, we have the average N 2 (t = 10 5 P 1 ) 2. Due to the stochastic nature of the evolution, all three simulations exhibit large variation in the individual N 2 , with about 2% of the runs having N 2 (t = 10 5 P 1 ) 10. Average number of CE2 encounters fiducial: no disk force drag = 10 5 P 1 drag = 10 6 P 1 trap = 10 5 P 1 trap = 10 6 P 1 drag = trap = 10 5 P 1 drag = trap = 10 6 P 1 Figure 12. Average number of CE2 events as a function of the time in the fiducial systems (see Fig. 1), including different types and strengths of disk forces. For example, the red solid line refers to the case that includes only F drag (with τ drag = 10 6 P1), the green solid line includes only F trap (with τtrap = 10 6 P1), and the blue solid line includes both F drag and F trap (with τ drag = τtrap = 10 6 P1). Figure 12 shows the effects of different disk forces on the CE2 rates by comparing the time evolution N 2 (t) from our simulations with different force types and strengths. The drag force (Equation 22) tends to stabilize the system by preventing crossing orbits between the BHs. When F drag with τ drag = 10 5 P 1 (dashed red) or 10 6 P 1 (solid red) is applied, the system is still unstable initially and the two BHs experience CE2 events (also see Li, Rodet and Lai, in prep). However, no more CE2s are found after about t = 10 4 P 1 for τ drag = 10 5 P 1 and t = 8 × 10 4 P 1 for τ drag = 10 6 P 1 . To check why CEs cease, Figure 13 compares the orbital separation of the two BHs in the no-drag simulations (fiducial, left panel) and the with-drag simulations (τ drag = 10 5 P 1 , right panel) at t = 4 × 10 3 P 1 , 1.2 × 10 4 P 1 , and 2 × 10 4 P 1 . In the fiducial simulations, r in = a in (1 + e in ) (the apocenter distance of the inner BH to the SMBH) and r out = a out (1−e out ) (the pericenter distance of the outer BH) spread to the region with r out − r in 2.5R H12 during the first 4000 orbits and remain in the same region at later time. This allows the two BHs to continue to "engage" with each other. However, in the τ drag = 10 5 P 1 simulations, the BH orbits evolve in time toward smaller r in and larger r out . Eventually, the difference between r out and r in becomes too large and CE becomes very rare. As a result, no CE2 happens between t = 2 × 10 4 P 1 and 10 5 P 1 in the with-drag simulations. with drag drag = 10 5 P 1 Figure 13. Apocenter of the inner BH rin = ain(1 + ein) and the pericenter of the outer BH rout = aout(1 − eout) from the fiducial (left) and the τ drag = 10 5 P1 simulations (right) at three different times. The region between the dashed lines corresponds to where (rout − rin)/RH12 ∈ (2, 3), with RH12 given by Equation (2) and assuming a1,2 rin,out. It is marked in both panels to help comparing the distributions.
One may expect that the trapping force (Equation 23) can accelerate the CEs by keeping the BH orbits close to the trapping radius. Our simulations show that when F trap is applied, our systems indeed have more CE0 events. However, Figure 12 shows that the CE2 rate is actually smaller than the fiducial simulations without F trap : The average number of CE2 events becomes 2.1 for τ trap = 10 6 P 1 (solid green) and 0.9 for τ trap = 10 5 P 1 (dashed green) after 10 5 orbits. Unlike in the simulations with F drag = 0, the CE2 rate in our F trap = 0 simulations do not drop to zero at later times.
The stabilizing effects of the drag force and the trapping force can balance each other. The blue curves in Figure 12 show that applying τ drag = τ trap = 10 5 P 1 produces more CE2 events than applying either of the two forces. With τ drag = τ trap = 10 6 P 1 , the two BHs experience slightly more CE2 events than in the fiducial simulations (black curve). Figure 14 shows that the r p -distributions of our simulations with drag forces are similar to the fiducial result. In particular, for r p /R H 0.01, the linear scaling of the distribution with r p remains accurate in all cases. In the simulations with the trap force, the mild encounters with r p /R H 0.01 are relatively more frequent than in the fiducial system with no disk forces (as the CDFs 10 4 10 3 10 2 10 1 10 0 r p /R H 10 4 10 3 10 2 10 1 10 0 CDF fiducial: no disk force drag = 10 5 P 1 drag = 10 6 P 1 trap = 10 5 P 1 trap = 10 6 P 1 drag = trap = 10 5 P 1 drag = trap = 10 6 P 1 Figure 14. Similar to Figure 3, but showing the cumulative distribution of rp in CE0 events for systems with different types and strengths of disk forces. with the trap force have a larger slope at r p /R H 0.01).
Thus the simple power-low distribution of r p given by Equation (14) is still valid for CE2 but becomes less accurate for CE0 when F trap is applied. In a real AGN disk, multiple effects can take place at the same time. The outcome of their competition depends on the detailed disk properties. Our results in this section suggest that disk forces have little effect on the r p -distribution for very close encounters (Figure 14), but can influence the CE2 event rate (Figure 12), and therefore affecting the BBH formation rate. Our simulations also suggest that in order for the two unstable BHs in a AGN disk to capture into a merging binary, a combination of the drag force and trapping force from the disk is needed, or we would rely on the chance that an individual system is in the high-N 2 tail of the distribution. We emphasize that our results in this section are based on simple prescriptions of disk forces. Hydrodynamics simulations will be needed to fully capture the effects of disk forces.

THREE AND MORE BHs AROUND SMBH
While a "SMBH + 2BHs" system is unstable only if a 2 − a 1 2 √ 3R H (see Equation 1), no precise stability criterion exists for systems with more than two BHs. In fact, such systems always exhibit instability eventually, except that the instability time grows with increasing orbital spacings. Numerical integrations suggest that a system of 3 or more bodies on nearly circular, coplanar orbits around a central massive object is stable for at least N orbital periods if the separation between adjacent bodies satisfies |a j+1 − a j | K(N )R H , where the constant K(N ) increases with increasing N (e.g., K 9 − 12 for N = 10 10 ; Smith & Lissauer 2009).
To explore close encounters in systems with more than two BHs, we consider a "SMBH + 3BHs" system with (m 1 , m 2 , m 3 ) = (2, 1, 0.5) × 10 −5 M and a "SMBH + 5 BHs" system with (m 1 , m 2 , m 3 , m 4 , m 5 ) = (2, 1, 0.5, 0.25, 0.125) × 10 −5 M , where M is mass of the central SMBH. The initial orbital separation is set as where R H,mut is mutual Hill radius of m j and m j+1 (see Equation 2) and we set K = 2 or 4. The BHs are given initial eccentricities e 1 = 0, e j = 10 −5 for j > 1 and inclinations i = R H /a 1 for all BHs. Similar to our fiducial simulations (Section 3), we carry out 2000 runs and sample the initial values of the argument of the peripasis, the longitude of the ascending node, and the mean anomaly randomly in the range [0, 2π] for each BH, assuming they all have uniform distributions.  Figure 15 shows the averaged cumulative CE2 rates for systems with two, three and five BHs. The two-BH (fiducial) systems experience nearly 3 to 5 times more CE2 events than the other systems. For the K = 2 cases, the average number of CE2 events at t = 10 5 P 1 is 0.9 with 3 BHs and 1.1 with 5 BHs. The fact that the CE2 rate is highest for the two-BH systems is somewhat surprising. It may be because the two-BH systems are subjected to conservation constraints, which cause the BH orbits to repeatedly overlap. A three-or-more-BH architecture has more degrees of freedom and there is no conservation law between any two BHs.
Adopting K = 4 in the three-and five-BH simulations leads to no CE2 events before t 2 × 10 3 P 1 . This is be-cause the initial conditions with larger K require longer times for the BHs to develop their first Hill sphere crossing. After the instability develops, the K = 4 systems catch up in the growth rate. At t = 10 5 P 1 , the three-BH and five-BH simulations with K = 4 yield N 2 = 0.5 and 0.62, respectively. For systems with larger K, we expect a longer time before the CEs occur, but the CE rate will eventually converge.
We fit the time evolution of the average number of CE2 events for t > 5 × 10 4 P 1 , by Least-square fit gives n 2 = 0.26, 0.50, 0.30, 0.50 and T = 1.2 × 10 5 , 4.0 × 10 5 , 8.5 × 10 4 , 2.6 × 10 5 for 3 BHs with K = 2, 3 BHs with K = 4, 5 BHs with K = 2 and 5 BHs with K = 4, respectively. Using N (t; r p < r cap ) = N 2 (t) P 2 (< r cap ) (see Equation 16) with P 2 (< r cap ) r p /(10 −2 R H ) (see Equation 14), we find that it would take about 5 × 10 12 P 1 , 4 × 10 9 P 1 , 4 × 10 11 P 1 , 3 × 10 9 P 1 on average to form a merging BH binary in each of the four cases.  Figure 16 shows that the r p -distributions of CE2 for all of the simulations in this section are similar to the fiducial result (see Figure 3). Figure 17 shows that CE2 always have nearly uniform distributions of cos i bin and thus the merging BH binaries formed in our scenario have a wide range of inclination angles with respect to the AGN disk plane. This is expected: Regardless of the number and the initial spacings of the BHs, the relative orbit of the captured binary is dominated by the gravity of the BHs during the CEs. Having more BHs and different initial spacings only affects the CE rate.

SUMMARY AND DISCUSSION
In this paper, we have studied the long-term evolution of two or more stellar-mass BHs on closely-packed, dynamically unstable circular orbits around a SMBH, with the initial semi-major axes separated by a few times the Hill radius R H . Such an orbital configuration may naturally arise from BH migrations in the AGN disk and leads to recurring close encounters between the BHs. We use N -body simulations to study the statistics and rate of close encounters (of various degrees of "closeness" compared to R H ), the properties of the relative orbits of the encountering BHs, and the probability for two BHs to be captured into a merging binary with the help of gravitational wave emission during very close encounters. Our fiducial simulations focus on "SMBH + 2 BHs" systems with m 1 = 2m 2 = 2 × 10 −5 M and initial inclinations i 1 = i 2 = R H /a 1 (where a 1 and a 2 are the BHs' semi-major axes around the SMBH). Additions simulations with different BH masses, initial orbital inclinations, prescriptions for disk forces, and different number of BHs are also performed.
Our simulations show that close encounters (CEs) between the BHs in such systems exhibit three general characteristics: 1. Close encounters (with the separation between two BHs less than R H ) occur stocastically and the rate of CEs generally declines in time (see Fig. 1). The average cumulative number of CEs, N (t) , is approximately a power-law function of time, although there are wide spreads in N (t) for different systems. For our fiducial "SMBH + 2 BHs" setup, we find N (t) ∝ t 0.5 (see Eqs. 9-11).
2. The vast majority of the CEs result in the formation of short-lived binaries with a small binding energy (of order Gm 1 m 2 /R H ; see Fig. 2). Such binaries are quickly destroyed by the SMBH tide and do not lead to binary mergers.
3. The closest separation (r p ) during a close encounter follows a cumulative distribution P (< r p ) ∝ r p for r p R H (see Figs. 3, 4, and Eq. 14). This distribution is robust, regardless of the BH masses (Section 3.5) and the number of BHs in the system (Section 6). For systems with very small but non-zero initial mutual inclinations, the same r p -distribution applies at later times as the mutual inclinations grow. In a system with two exact co-planar BHs, the distribution becomes P (< r p ) ∝ r 1/2 p ( Fig. 8 and Eq. 19).
These results imply that, to capture an encountering BH pairs into a "permenant" binary, a fast dissipative mechanism is required. Given the high close encounter rate, a promising mechanism is the GW emission when r p is very small (less than the critical "capture" radius r cap ∼ 10 −4 R H , depending on the system parameters; see Eq. 6): 4. We provide a semi-analytical formula for the averaged cumulative rate of binary captures via GW emission in a "SMBH + 2 BHs" fiducial system, Equation (16) or (17), which is the product of the close encounter rate N (t) and the capture probability P (< r cap ) during each encounter.
5. Our formula suggests that, in the fiducial systems, the timescale for two BHs to be captured is 10 8 P 1 on average (assuming r cap = 10 −4 R H ) and 10 7 P 1 for the 5% of the systems with N (t) 4 N (t) . In the exact co-planar systems, the capture rate is much higher (see Eq. 20), and we find that the average binary capture time is 6 × 10 3 P 1 .
6. After the two BHs are captured, we expect these BH binaries to merge in a few binary orbits. Their mergers will exhibit high eccentricities (see below) when entering the LIGO band ( 10 Hz) and will have a board distribution of orbital inclinations relative to the original AGN disk (see Figure 5).
We have carried out additional simulations to assess how the above results may be influenced by various system parameters and by the gas disk effects: 7. The masses of the BHs (relative to the mass of the SMBH) affect the close encounter rate only in a modest way (Section 3.5). For example, the rate of CEs with r p ≤ 10 −2 R H decreases by a factor of ∼ 2 when the BH masses are lowered by a factor of ten (see Figure 6). The r p -distribution (Eq. 14) is robust. Thus, we expect that our fiducial binary capture rate, Equation (17), remain valid for systems with different BH masses.
8. The most optimal setup to get more binary captures in our scenario is with two exact co-planar BHs (see points 3 and 5 above). Such an exact co-planar system would have a higher rate of close encounters (see Fig. 7 and Eq. 18) and a flatter r p -distribution ( Fig. 8 and Eq. 19), leading to a greatly enhanced binary capture rate (Eq. 20). However, if the mutual inclination between the BH orbits is initially small but non-zero, it will inevitably grow and evolve to an "equilibrium" value ∼ R H /a 1 (see Fig. 9), causing P (< r p ) to converge to the fiducial result (see Fig. 10).
9. We have explored the effects of gas disks by applying simple prescriptions of disk forces (see Eqs. 22 and 23) on the BHs in our N-body simulations. We find that such prescribed disk forces do not necessarily lead to an enhanced binary capture rate. In fact, simple gas drags on the BHs may stop close encounters at late times. A "migration trap" force can sometimes balance the drag force and maintain the close encounter rate (see Fig. 12).
10. The number of BHs and their initial spacings in a closely-packed system do not affect the r pdistribution ( Fig. 16) and the orbital obliquity distribution ( Fig. 17) during the very close encounters. We find that systems with more than 2 BHs have a lower close encounter rate than systems with two BHs.
Regarding point 6 above: Using Equation (6) for r cap , we find that the GW emitted by the BH binary at capture has a frequency For the adopted parameter values, this frequency lies slightly below the LIGO band. Such a newly captured binary will merge within a few binary orbits. Using the standard gravitational radiation formulae (Peters 1964), it is easy to see that when the frequency of the GWs from the binary enters the LIGO band (¿10 Hz), the binary can retain a very significant ( 0.5) eccentricity. The event rate of binary BH mergers produced in our scenario depends on the population of BHs in AGN disks, which is very uncertain. Nevertheless, given the binary capture timescale obtained in this paper, a nonnegligible event rate for such BH mergers may be expected. Perhaps the recently claimed eccentric merger event GW190521 (Gayathri et al. 2022) is an example.
The most uncertain aspects of the present study concern the effects of gas disks. Our conclusion on the gas effects (point 9 above) should be considered tentative. For future work, hydrodynamics simulations should be used for a more in-depth study of the disk effects. Although we have only considered GW emission in this paper, the broad range of r p in close encounters may allow the physical processes at different distance scales to operate. For example, interaction between the circum-BH disks may generate dissipation for the BH pairs during close encounters. Our results of the close encounter rate and the r p -distribution can be applied to these alternative dissipation mechanisms to explore different possibilities of binary BH formation in AGN disks.