Dynamical Friction Models for Black Hole Binary Formation in Active Galactic Nucleus Disks

Stellar-mass black holes (sBHs) embedded in gaseous disks of active galactic nuclei (AGN) can be important sources of detectable gravitational radiation for LIGO/Virgo when they form binaries and coalesce due to orbital decay. In this paper, we study the effect of dynamical friction (DF) on the formation of BH binaries in AGN disks using N-body simulations. We employ two simplified models of DF, with the force on the BH depending on Δ v , the velocity of the sBH relative to the background Keplerian gas. We integrate the motion of two sBHs initially on circular orbits around the central supermassive black hole (SMBH) and evaluate the probability of binary formation under various conditions. We find that both models of DF (with different dependence of the frictional coefficient on ∣Δ v ∣) can foster the formation of binaries when the effective friction timescale τ satisfies ΩK τ ≲ 20–30 (where ΩK is the Keplerian frequency around the SMBH): prograde binaries are formed when the DF is stronger (smaller τ), while retrograde binaries dominate when the DF is weaker (larger τ). We determine the distribution of both prograde and retrograde binaries as a function of initial orbital separation and the DF strength. Using our models of DF, we show that for a given sBH number density in the AGN disk, the formation rate of sBH binaries increases with decreasing τ and can reach a moderate value with a sufficiently strong DF.

Several mechanisms have been proposed for facilitating the mergers of BBHs in AGN disks.For example, the dynamical evolution of BBHs can be affected by the gravity from the highly perturbed gas around them (e.g., Baruteau et al. 2011;Stone et al. 2017); high-resolution hydrodynamics simulations have shown that, for some AGN disk parameters, such BBH-disk interaction may indeed help the BBHs contract in orbit (and possibly merge; Li et al. 2021Li et al. , 2022b;;Dempsey et al. 2022;Li & Lai 2022, 2023); and recent three-dimensional hydrodynamics simulations suggest that retrograde BH binaries may undergo rapid gas-driven mergers and experience runaway eccentricity growth (Calcino et al. 2023;Dittmann et al. 2023).BBHs embedded in an AGN disk may also be hardened and driven to mergers by a series of scatterings with other single BHs (e.g., Stone et al. 2017;Leigh et al. 2018;Samsing et al. 2022), especially when the mutual orbital inclinations between the BHs are small.These mechanisms, however, all assume preexisting binaries.The origin of sBH binaries in AGN disks remains an open question.While some binaries may be formed in nuclear star clusters and then captured into the AGN disk (Bartos et al. 2017), it has been suggested that many binaries may form within the AGN disk through close encounters between embedded single BHs (e.g., Stone et al. 2017;Fabj et al. 2020).Since AGN disks may have nonmonotonic radial density profiles (Sirko & Goodman 2003;Thompson et al. 2005;Dittmann & Miller 2020;Gilbaum & Stone 2022), single BHs embedded in a disk may naturally approach each other due to differential migrations or migration traps (Bellovary et al. 2016;Secunda et al. 2019Secunda et al. , 2020)).The resulting tightly packed BH orbits can be dynamically unstable, leading to frequent close encounters and opportunities for binary formation.
A few recent studies have investigated the dynamics and the probability of forming tightly bound BBHs through close encounters.Assuming negligible interactions between the BHs and the gaseous disk, Li et al. (2022a) showed by a suite of longterm N-body simulations that BH-BH encounters almost always lead to weakly bound and short-lived binaries; only in the rare event of very close encounters, which may occur over longer timescales, can the two BHs be captured into tightly bound Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
binaries via GW emission.Such binaries are formed with very high eccentricities and remain eccentric when entering the LIGO band.Combining the simulation results with a theoretical scaling analysis, they obtained the probability of binary formation via this "GW bremsstrahlung" mechanism (see Rom et al. 2023 for an analytical study).Binary formation without gas was also studied by Boekholt et al. (2023) using a phase space fractal structure analysis.Tagawa et al. (2020b) estimated that most merging BBHs in AGN disks are formed in close encounters through the dissipation of kinetic energy in the disk gas, although their results are based on one-dimensional N-body simulations and rough semianalytical prescriptions (i.e., with simplified binary evolutions and BH-disk interactions).Li et al. (2023) used high-resolution global hydrodynamics simulations to study the close encounters between two single sBHs in AGN disks.They showed that when the disk density is sufficiently high, bound BH binaries can be formed by the collision of circumsingle disks.They also found that the resulting binaries tend to have compact semimajor axes, large eccentricities, and retrograde rotations.A similar investigation using smoothed particle hydrodynamics simulations was done by Rowan et al. (2023), who found that prograde binaries can also be formed.Both Li et al. (2023) and Rowan et al. (2023) have only surveyed limited parameter spaces because their simulations are very expensive.DeLaurentiis et al. (2023) studied the BH-BH close encounters in AGN disks using a large number of N-body simulations that model the gas effect with simple gas dynamical friction (GDF; Ostriker 1999).They characterized the properties of the resulting binaries as functions of the impact parameter of the close encounters and explored the dependence of BBH formation on various parameters such as BH mass, disk temperature, and disk density.Recently, Whitehead et al. (2023) carried out extensive hydrodynamical simulations in a shearing box and determined the condition of binary formation as a function of the impact parameter of the two sBHs and the gas density of the background AGN disk.
Although previous works have found that BH-BH encounters can form both prograde and retrograde binaries, a complete picture of the conditions and the rates of their formation is still missing.Prograde and retrograde binaries (relative to the original AGN disk) can affect the spin-orbit misalignment angles of the merging BBHs.It has also been recognized that, compared to the prograde ones, retrograde binaries may inspiral at faster rates and can grow their masses more rapidly through gas accretions (Li et al. 2021;Li & Lai 2022).Hence, understanding the conditions and rates for the formation of prograde and retrograde binaries is important for assessing the AGN channel of BBH mergers.
In this paper, we study the gas-assisted formation of prograde and retrograde BH binaries in an AGN disk through N-body simulations.Using two different gas force models, we systematically investigate the binary formation conditions as functions of the BHs' initial orbital separation and masses and the strength of the gas forces.Our simulations employ a similar setup as in some previous works (e.g., Li et al. 2022a;DeLaurentiis et al. 2023).However, we focus on the results that can be scaled to various parameters/conditions, and our calculations reveal the different pathways and probabilities to form binaries with different rotational orientations.
The rest of this paper is organized as follows.In Section 2, we describe the setup of our N-body simulations and the two different gas force models, including the linear frictional force and the GDF.Sections 3 and 4 show the results of the N-body simulations with the linear frictional force and the GDF, respectively.We use a local shearing box model to analyze our simulation results in Section 5. We include GW emission in Section 6 and examine the relative roles of gas friction and GW bremsstrahlung.Then, based on these results, Section 7 estimates the formation rates of prograde and retrograde binaries in AGN disks (given the number density of sBHs).Section 8 summarizes our findings.

Simulation Setup and Initial Conditions
We consider a system consisting of a central SMBH (mass M) and two sBHs (m 1 and m 2 ) that orbit around M on initially separate, nearly circular orbits.The sBHs are embedded in a gaseous accretion disk around the SMBH and subjected to gas drags.We investigate how the two sBHs dynamically form a bound binary with the assistance of the gas.
The two sBHs are given masses of m 1 = 10 −6 , m 2 = 0.5 × 10 −6 M for most of this paper.They are initially placed on coplanar circular Keplerian orbits around the SMBH with orbital radii a 1 and a 2 .We set the initial orbital separation of m 1 and m 2 using the dimensionless parameter K, defined as is the mutual Hill radius of the two sBHs with m 12 = m 1 + m 2 .We also specify the initial azimuthal separation of the sBHs by the angle Δf 0 = f 2 − f 1 , where f 1 , f 2 are the initial orbital phases.
We choose the values of K from the range [1.0, 2.5] and Δf 0 from [10R H /a 1 , 2π − 10R H /a 1 ].These values allow m 1 and m 2 to be well separated (?R H ) initially and become much closer (R H ) at a later time.The sBHs are then expected to experience a close encounter and can possibly form a binary.Note that for K outside this range, the orbits of the sBHs are either weakly perturbed (for K  2.5) or experience strong deflection on a horseshoe trajectory (for K  1; see, e.g., Petit & Henon 1986), and binary formation is nearly impossible, except for a few cases enabled by multiple close encounters (see Section 3).

Running the Simulation
We simulate the orbital evolution of the "SMBH+2sBHs" systems using the N-body software REBOUND (Rein & Liu 2012) and the IAS15 integrator (Rein & Spiegel 2015).The calculations are repeated for many values of K and Δf 0 .Each initial condition is run for 150P K , where experience about two close encounters for most cases (see Section 3.3; also see Figures 12 and 13 of Li et al. 2022a).
The gas effects are implemented in our N-body simulations as additional forces on the sBHs.We apply two models of gas forces, as described in the next subsection.

Models for Gas-BH Interaction
In this work, we consider two models for the gas drags experienced by the sBHs in the Keplerian AGN disk.

Model 1
The first model applies a linear frictional force.The force per unit mass is given by where v is the velocity of the sBH and v K is the local Keplerian velocity of the disk around the SMBH.The strength of this friction is determined by the timescale τ.

Model 2
The second model is the canonical GDF (Ostriker 1999).We express the GDF as which is similar to model 1 except that the friction timescale is not a constant.The variable x ≡ |v − v K |/c s is the Mach number of the relative velocity between the sBH and the background gas, while the coefficient is a characteristic damping timescale determined by the mass of the sBH (m), the gas density ρ, and the sound speed c s .Because we consider sBHs on closely packed orbits, the dynamical evolution of the sBHs is confined within a small extent where the local disk profile is approximately constant.Therefore, we assume a constant τ 0 (but m-dependent) for each sBH and a constant disk aspect ratio h = c s /|v K | for each of our simulations (i.e., τ 0 is independent of the sBH position).
The function f (x) is defined piecewise in the subsonic and supersonic regimes as (Ostriker 1999) In practice, to achieve fast convergence of numerical integrations, we adopt two approximations.First, when x → 0, we note that both f (x) and x 3 approach 0 but x 3 /f (x) → 3. Thus, we set τ(0) = 3τ 0 and adopt a linear interpolation to calculate τ(x) between x = 0 and 0.01.Second, to avoid the discontinuity in f (x) at x = 1, we linearly interpolate τ(x) between x = 1 − δ and 1 + δ based on the exact values of τ at x = 1 − δ and 1 + δ.To ensure the smoothness of τ(x), we choose δ = 0.04.The smoothed function τ(x)/τ 0 decreases from 3 at x = 0 to ∼0.5 at x ∼ 1 and increases to ∼8 at x ∼ 3.

N-body Simulation Results: Binary Formation with Gas Friction Model 1
In this section, we run N-body simulations for the "SMBH +2sBHs" systems with the gas frictional force model 1.
Figure 1 shows the results of sBHs in five example runs.Each row represents one simulation: the leftmost panel shows the relative trajectory of the sBHs, and the second left, second right, and rightmost panels depict the specific angular momentum, semimajor axis, and eccentricity of the sBHs' relative orbit, respectively.
These examples show two apparently distinct outcomes: the sBHs either form (first and third through fifth rows) or do not form (second row) bound binaries.Among the formed cases, both prograde (first, third, and fifth rows) and retrograde (fourth row) binaries are observed.
Clearly, the outcomes of the close encounters may depend on the initial conditions (K, Δf 0 ), the masses of sBHs m 1,2 , and the gas damping timescale τ.Therefore, we conduct parameter studies to investigate the effect of varying these parameters.

Criteria for Binary Formation
Before showing the results of our parameter study, we first discuss our prescription to identify binary formation and determine the sense of rotation (prograde/retrograde) of the formed binaries.This is needed to efficiently categorize the simulation results in the parameter survey.
Given our choice of 150P K as the integration time, the two BHs would have sufficient interaction time after their first close encounter.We consider the binaries to have formed if the following two conditions are both met.
(i) The separation between the two BHs remains smaller than R H until the end of integration, after reaching R H at some time.
(ii) The relative energy of m 1 , m 2 becomes negative at some time and remains negative until the end of integration, i.e., where μ = m 1 m 2 /m 12 is the reduced mass and r 1,2 and v 1,2 are the position and velocity vectors of m 1,2 .
In practice, the orbital integration becomes more timeconsuming when the separation |r 1 − r 2 | gets too small.To avoid such a situation, we use a third condition to terminate the integration and directly conclude that the binary has formed once the condition is satisfied.
(iii) The relative semimajor axis of m 1 , m 2 is positive and less than 0.1R H , i.e., This condition implies that the relative apocenter separation of the BHs is less than 0.2R H .We have verified that the above conditions provide good criteria for discerning whether the two BHs form stable binaries in our numerical integrations.
We determine the orientation of the binaries by computing the specific relative angular momentum of the sBHs, With z being the unit vector along the sBHs' initial orbital angular momentum around the SMBH, a binary is considered prograde if (i.e., the binary's relative orbit is in the same direction as its center-of-mass orbit about the SMBH) and retrograde if In some simulations, the sign of ℓ rel might change a few times before stabilizing (see the second left panel of the bottom row of Figure 1.)We determine the orientation of a binary based on the direction of l rel at the end of the integration time.

Dependence on K and τ
We run a suite of simulations with the initial K in the range [1, 2.5], 1/τ in the range [0, 0.4Ω K ], and the initial Δf 0 = π/2.Figure 2 shows the result.Each dot represents one integration with a different combination of K and τ; a black dot means a binary is not formed, a blue dot means a prograde binary is formed, and a red dot means a retrograde binary is formed.
This result suggests that the formation outcome strongly depends on the initial K and the friction timescale τ.BH binaries are only formed when 1/τ exceeds a critical value, given by All binaries are retrograde when 1/τ is approximately between 0.03Ω K and 0.13Ω K .Prograde binaries appear when 1/τ  0.13Ω K .Most binary formations are formed in systems with the initial K ä [1.7, 2.4].For initial K  1.7, binaries are only formed in a few systems with some special combination of K and τ (see Section 3.3).

Dependence on Δf 0
To study the Δf 0 dependence of our results, we run a suite of simulations with initial Δf 0 ä [10R H /a 1 , 2π − 10R H /a 1 ], 1/τ between [0, 0.4Ω K ], and initial K = 2. Figure 3 shows the result, which suggests that the formation outcome depends weakly on Δf 0 .
This result matches the expectation.For sBHs with initial Δf 0 ä [10R H /a 1 , 2π-10R H /a 1 ], their motion remains approximately Keplerian until the mutual separation decreases to a few R H (i.e., when Δf(t)  10R H /a 1 ).Therefore, Δf 0 should only affect the timing of the sBH-sBH close encounters but not their outcomes.

Dependence on m 1,2 /M
We check the dependence of our results on the masses of the sBHs by repeating the simulations in Section 3.2.1 and Figure 2 with m 1 = 10 −3 M and m 2 = 0.5 × 10 −3 M. The result is shown in Figure 4.
Comparing Figures 2 and 4, we find that the two are similar: the distribution of formed binaries is only slightly shifted to smaller K when m 1 and m 2 increase by 3 orders of magnitude.In Section 5, we demonstrate that the simulation results can be reproduced by a set of equations of motion that are independent of sBH masses (see Equations ( 14) and (15)).Indeed, in the shearing box local approximation, which is accurate to the first order in R H /a 1 , the binary formation outcome is independent of the sBH-to-SMBH mass ratio.Hence, in our setup, the  4)).Each dot represents an integration.The black, blue, and red dots represent the distribution of simulations that do not form binaries, form prograde binaries, and form retrograde binaries, respectively.difference in the formation outcomes for different m 1,2 /M can be expected to be small.

Multiple Close Encounters
In Figure 2, most binary formation occurs in systems with initial K  1.7.When K  1.7, binary formation can only happen inside the narrow fingerlike regions around some particular values of K.
We find that binaries with initial K  1.7 are formed during the first BH-BH close encounter, while binaries in the fingerlike regions are formed via multiple encounters.Figure 5 shows four examples in the latter case.During their first encounter, the two sBHs avoid each other on a horseshoe orbit (see the left panels of Figure 5). 4The two BHs are then quickly circularized by the gas drag on two new orbits with a large orbital separation (see the right panels of Figure 5).Therefore, the sBHs will have a larger "initial" K during their second close encounter and hence can form a binary just like those cases with a true initial K  1.7.This also explains why the substructures of each fingerlike region are qualitatively similar to the substructures in the main binary formation region at large K.We expect similar fingerlike regions to also appear at K < 1 (outside of our parameter space), corresponding to binary formations via even more encounters.However, at smaller K values, those formation regions are expected to be narrower, and the formation timescales are expected to be longer due to the longer synodic period.Hence, binary formations via even more encounters may only have a negligible effect on the total binary formation rate.

N-body Simulation Results: Binary Formation with Gas
Friction Model 2 In this section, we conduct a series of parameter studies that are similar to those in Section 3.2 but now using the GDF (force model 2).Since the results in Section 3.2 suggest that the effects of varying m 1,2 and Δf 0 are weak, we focus on exploring how the formation outcome depends on K and τ(x).
To make fair comparisons with the results obtained with force model 1 (Section 3.2 and Figure 2), the simulations in this section also use m 1 = 10 −6 M, m 2 = 0.5 × 10 −6 M, initial Δf 0 = π/2, and initial K between [1, 2.5].The characteristic damping timescale of the GDF τ 0 (see Equation ( 7)) is chosen such that 1/τ eff ä [0, 0.4Ω K ], where When the two sBHs enter their mutual Hill sphere, their relative velocity is of order Ω K R H . Thus, we expect that τ eff is qualitatively equivalent to τ in model 1.
Figure 6 shows the results of the parameter studies using model 2 (GDF).Each panel of Figure 6 adopts a different h = c s /v K , the disk aspect ratio, but they all demonstrate similar distributions of the formation outcomes.These results are qualitatively similar to the model 1 results (Figure 2) in terms of the shapes of the formation regimes, the distribution of the prograde and retrograde cases, and the required values of τ and τ eff for binary formation.
In all cases, BH binaries can be formed when K this is similar to the result for model 1 (see Equation (12)).Higher formation rates require shorter damping timescales.Figure 7 shows the τ 0 profiles in four different AGN disk models.The two profiles in the upper panel, originally from Sirko & Goodman (2003) and Thompson et al. (2005), are for SMBHs with M = 10 8 M ☉ , while the lower panel represents the two models from Dittmann & Miller (2020), with M = 4 × 10 6 M ☉ .
In the Sirko & Goodman (2003) disk model, the damping timescale is around a few W - K 1 at radius R (or the semimajor axis of the sBHs) between [4 × 10 2 R g , 10 4 R g ], where R g = GM/c 2 is the gravitational radius of the SMBH.Binaries may form within this region.However, the formation rate may be low except near R = 10 3 R g , where τ 0 Ω K  1.This high-rate location corresponds to the radius where the disk has the smallest aspect ratio h and the largest surface density Σ = 2ρhR.
The Thompson et al. (2005) profile has Thompson et al. (2005) is less dense than in Sirko & Goodman (2003), it also has an aspect ratio h ∼ 10 −3 for R ä [3 × 10 2 R g , 10 4 R g ], which is about 10 times smaller than in Sirko & Goodman (2003).Therefore, the damping timescale τ 0 can be very short in this thin middle region, leading to higher binary formation efficiency.
Dittmann & Miller (2020) consider disks around low-mass AGNs at high redshifts.Based on their models, τ 0 Ω K are less than unity from R ∼ 10 3 R g to ∼10 6 R g .BH binaries can be formed in the middle region of the disk, as well as in the outer part of the disk that connects with the star-forming region.This may allow more progenitor BHs to participate in dynamical encounters and produce a larger population of binaries.

Equation of Motion in Shearing Box
In addition to using direct N-body simulations, an alternative method to study the sBH-sBH close encounters and binary formation is to use the shearing box model.This model applies to the situation when the separation between the BHs is much smaller than the distance to the SMBH (i.e., |r 1 − r 2 | = a 1 , a 2 ), and it solves their relative motion in the corotating frame anchored on the center of mass of m 1 and m 2 .
We assume that the m 1 -m 2 center of mass is circulating around the SMBH at a constant angular velocity 3 , which is approximately true when |r 2 − r 1 | is small and m 1,2 = M. Keeping terms up to the first order of |r 2 − r 1 |/a 1 , we obtain the shearing box equations of motion for a BH binary subjected to the gas force: where x and y are the radial and azimuthal components of (r 2 − r 1 )/R H , ˆ t t t = W W com K is the dimensionless dissipation timescale, and the dot represents the normalized timederivative operator d/d(Ω K t).In this section, we consider the gas force model 1, so that t is constant.Equations ( 14) and (15) are the same as the Hill equations but with the additional terms from the gas force.A similar set of equations has also been used by Schlichting & Sari (2008a) to study the formation of Kuiper Belt binaries.
The initial condition for the "SMBH+2sBHs" systems that we adopt in the previous sections (see Section 2.1) can be implemented in the shearing box model as Equations ( 14)-( 19) do not depend explicitly on m 1 , m 2 , or M, which is consistent with the N-body simulation results in Section 3.2.3.Nevertheless, we are aware that this shearing box model is only accurate when R H is at least 2 orders of magnitude smaller than a 1 , i.e., R H  10 −2 a 1 .This is because (i) the initial velocities of m 1 and m 2 are Keplerian only when the separation between them is at least a few times R H , i.e., |r|  10R H , and (ii) the shearing box approximation requires |r| = a 1 .To meet both requirements, the BHs need to have 2 .This is equivalent to requiring m 12 /M  10 −6 , which is accordant with the canonical BH masses in our N-body simulations (see Sections 3 and 4).
To apply the binary formation criteria (see Section 3.1) in the shearing box model, we compute the relative energy of m 1 and m 2 as where r and v are the relative position and velocity vectors of m 1 and m 2 in the corotating (shearing box) frame.(Note that v is different from v 2 − v 1 in Equation (9).) Equation ( 20) can be further simplified by realizing that  W W z com K is perpendicular to the plane in which r lies.In dimensionless form, the relative energy is ( )

. Parameter Study
We run a suite of shearing box calculations with the initial condition , 10, 0, 3 2 22 0 for K between [1, 2.5] and t-1 in the range of [0, 0.4].These ranges for K and t are to mimic the parameter sweepings in Section 3.2.1.The value y 0 = 10 corresponds to Δf 0 = 10R H /a 1 , which ensures 10R H  |r 2 − r 1 | = a 1 initially.We fix the value of y 0 in all of our shearing box calculations.This is because of the weak dependence on Δf 0 seen in our N-body results; we expect our shearing box model results to be approximately independent of y 0 .
Figure 8 shows the results.The blue, red, and black dots represent the formation of prograde binaries, the formation of retrograde binaries, and no binary formation, respectively.Their distribution qualitatively reproduces the results of Nbody simulations.The only difference is that the shearing box does not allow any binary formation at K  1.4 (i.e., no fingerlike regions as seen in Figure 2) because the shearing box only describes the first close encounters between the BHs.

General Description
The distribution of formed binaries in the K-versusdissipation-rate parameter space can be divided into four regions, corresponding to four types of qualitatively different formation trajectories.We label the four regions in Figure 8. Region I corresponds to the peninsula on the left at K ≈ 1.5.This region has a lower limit at around t » 1 0.25and is filled with prograde binaries.The rest of the formed binaries are in region II, which can be further divided into three subregions: region IIb refers to the wide band of retrograde binaries at around K ä [1.8, 2.2], and regions IIa and IIc are the left and right bands of mostly prograde binaries adjacent to region IIb.When t 1 is between 0.05 and 0.15, all binaries inside region II are retrograde.
We analyze two questions regarding these regions: (i) what determines the binary orientation (prograde versus retrograde) in each region and (ii) why there are only retrograde binaries for ˆ t 1 0.15.

Binary Orientation in Each Region
Figure 9 shows some typical trajectories of the BHs in our shearing box integrations.In the gas-free regime ( t = 1 0; see the left panel of Figure 9), the dynamical evolution of the BHs is described by the well-known "satellite-encounter" problem (e.g., Petit & Henon 1986).14) and ( 15)) with different initial orbital separations K and frictional damping timescales t.All calculations adopt y 0 = 10.The black, blue, and red dots represent the systems that do not form binaries, form prograde binaries, and form retrograde binaries, respectively.The distribution of the outcomes can be divided into regions I, IIa, IIb, and IIc.BHs evolving from regions I, IIa, and IIc eventually become prograde binaries, while those from region IIb form retrograde binaries.
Figure 9. Left: typical trajectories in our shearing box calculations in the gas-free regime ( t = 1 0).See the text for the descriptions of the six types of trajectories.We do not include boundary cases (e.g., some messy trajectories between the blue and green curves) for conciseness.Right: examples from the four binary-forming parameter-space regions depicted in Figure 8 (with t = 1 0.4).The incoming trajectories (i.e., the precapture portion of the curves) are similar to those in the gas-free cases (highlighted in the left panel: blue → region I, green → region IIa, red → region IIb, and black → region IIc).
1.For K  1.6 (gray) trajectories, m 2 moves in a horseshoe orbit without getting very close to m 1 .2. For K around 1.65 (blue), m 2 is attracted toward m 1 in the counterclockwise direction (prograde), then moves through the left portion of the Hill sphere.3.For K around 1.85 (green), m 2 is scattered by m 1 from the top right; the close encounter between the two BHs is almost head-on.4. For K around 2.0 (red), m 2 spirals toward m 1 from the right, then completes a full clockwise rotation (retrograde).5.For K around 2.27 (black), m 2 initially passes m 1 from the right but then makes a U-shaped turn toward m 1 due to the tidal force and the Coriolis force.Eventually, m 2 flies by m 1 at a short distance in the counterclockwise direction (prograde).6.For K  2.35 (gray), m 2 passes m 1 on an outer orbit without any close encounters.
In all cases above (with no friction), m 2 eventually leaves m 1 , and no binaries are formed.However, the directions of these close encounters can determine the orientation of the binaries (prograde versus retrograde) when the gas force is added.
The right panels of Figure 9 show the dynamical evolution of the BHs in the four parameter-space regions (with t = 1 0.4).The BHs approach each other on trajectories similar to those in the gas-free cases.As they become tightly bound due to the gas, BHs from regions I, IIb, and IIc inherit the directions of their encounters as the orientations of the new binaries (i.e., prograde, retrograde, and prograde, respectively).The BHs from region IIa experience a near head-on encounter, forming a binary with a very small initial angular momentum.Before the gas fully stabilizes the binary, m 2 first travels to the far-off apocenter on the right.It is then deflected by the Coriolis effect to the lower right quadrant, where the binary relative energy is dissipated by the tidal force.Hence, as m 2 returns, it orbits m 1 in the prograde direction as a bound binary.
In summary, the orientation of a binary is determined by the type of its encounter trajectory, which is controlled by the initial orbital spacing K.We also note that strong gas damping can increase the "effective K" between two BHs, so the four regions shift to smaller K as t 1 increases.

Dissipation Timescale Requirement for Prograde and Retrograde Binaries
From Figure 8, we see that essentially all5 binaries are retrograde when ˆ t 1 0.15.Similar results have been reported by Schlichting & Sari (2008b) for the formation of Kuiper Belt binaries due to (particle) dynamical friction. 6chlichting & Sari (2008b) suggest that this is because the retrograde binaries are more stable than the prograde ones.One possible explanation for the stability difference is the Coriolis force; during the binary formation process, the Coriolis force always points outward for prograde binaries and inward for retrograde binaries, and this tends to destabilize and stabilize their orbits, respectively (see, e.g., Hamilton & Burns 1991;Hamilton & Krivov 1997).This difference in the stability could imply that, for BH pairs that are trying to form binaries, those undergoing retrograde encounters are more likely to succeed.
To gain some insight on prograde and retrograde binary formation, we rerun our shearing box calculations for 2, 2.2 0.05, 0.15 , where all bound binaries are retrograde; however, at some point during the binary formation process, we reverse the relative velocity of the BHs from ( ) x y , and check whether they can form prograde binaries instead.Our sample contains 822 systems that would form stable retrograde binaries without the velocity change.When we reverse the relative velocity of the BHs during their first pericenter passage, we find that 407 systems remain bound as wider retrograde binaries, and the rest become unbound.If we wait for the BHs to pass the pericenter completely and then reverse their velocity when their separation + x y 2 2 becomes less than 0.25, 6 out of 822 systems can form stable prograde binaries, and the rest become unbound.If we reverse the velocity at the first apocenter, all systems fall apart.This experiment shows that the prograde binaries are indeed less stable than the retrograde ones.Hence, prograde binaries require stronger gas friction to form, while retrograde ones can still form when the friction is weaker.

GW-emission-assisted Formation of Binaries
In previous sections, we have considered binary formation when the only dissipative process is the gas drag.However, when the sBHs reach a sufficiently small separation, GW emission becomes important.Two encountering sBHs can form a tight binary via GW bremsstrahlung when their pericenter distance r p becomes less than a critical value r GW (Li et al. 2022a), where η is a constant of order unity.For our fiducial sBHs at a 1 ∼ 10 3 R g , at which the gas-assisted binary formation is efficient (see Figure 7), r GW ∼ 2 × 10 −5 R H .

Method
We rerun our fiducial shearing box calculations incorporating the GW-assisted binary formation; while we simulate encountering sBHs with Equations ( 14) and (15), we immediately consider them as a bound binary once their mutual separation becomes smaller than r GW = 10 −5 R g .For simulations with (r 2 − r 1 )/R H > r GW at all times, we use the same binary formation criteria as in the fiducial shearing box calculation and the N-body studies (see Section 3.1).

Results
Figure 10 shows the shearing box calculation results when we include binary formation by GW emission.In addition to the blue, red, and black dots defined in the same way as in Figure 8, the deep green and light green dots represent prograde and retrograde binaries formed due to GW emission, respectively.Most of these cases exhibit highly chaotic orbital trajectories, and the sBHs can sometimes undergo multiple pericenter passages before they enter the GW-capture regime.
In the absence of gas friction ( t = 1 0), binary formation by GW emission occurs only when K is around two narrow ranges centered at K ; 1.87 and 2.2, corresponding to the two cases where the two sBHs experience close encounters with nearly zero relative angular momentum (see the left panel of Figure 9).Boekholt et al. (2023) have shown that within these narrow ranges, the character of the trajectory depends sensitively on the precise value of K and exhibits a fractal-like structure.
In Figure 10, the cases where the sBHs enter the GWdominated regime (r p  10 −5 R H ) during their first pericenter passage are represented by two thin green lines passing through ( ˆ) ( ) t » K, 1 1.9, 0 and (2.2, 0).However, most GW-induced binary formation occurs during the second or later pericenter passages.These binaries are distributed at the boundaries between the gas-dominated prograde and retrograde regions.One can think of them as binaries that are initially "precaptured" by the gas drag (like their prograde and retrograde neighbors) and then enter the GW regime due to their low angular momentum.

Binary Formation Rate
In the previous sections, we have obtained the probability of binary formation as a function of K and the friction timescale τ.These results can be applied to a "realistic" situation where many sBHs populate the AGN disk and undergo two-body close encounters.Let n be the number density (per unit area) of sBHs in a certain local region of the AGN disk.The binary formation rate per unit area is  Similar equations can be used to estimate the rates of forming prograde and retrograde binaries.
The left panel of Figure 11 shows the formation rate estimated by Equation (26) using the F(K, τ) from Figure 2 (i.e., N-body model 1), while the right panel shows the formation rate estimated from Figure 8 (i.e., shearing box results).The two panels are similar, since the shearing box results are the same as the model 1 N-body results except for the binaries formed through multiple close encounters.Retrograde binaries begin to form at around 1/τ ; 0.03Ω K (see Equation ( 12)).Their formation rate reaches around W n R 0.5 K 2 H 2 at 1/τ = 0.1Ω K and saturates at around for larger 1/τ.Prograde binaries appear at 1/τ  0.13Ω K ; their formation rate grows as the gas force increases and surpasses the retrograde rate at around 1/τ  0.32Ω K .The total formation rate is  8 except that here we include binary formation by GW emission in addition to the gas drag.The blue and red dots represent binaries formed without triggering GW emission, i.e., the only dissipative process being gas drag.The deep green and light green dots represent prograde and retrograde binaries formed due to GW emission, respectively.The labels ending with "PP" on the green dots indicate during which episode of pericenter passage the system reaches |r 1 − r 2 | < r GW = 2 × 10 −5 R H and forms a bound binary due to GW emission.roughly a linear function of 1/τ.We also calculate these rates based on Figure 6 (gas force model 2) and show the results in Figure 12.All these results are qualitatively similar, except that the h = 0.01 case has a slightly lower overall formation rate than the others.
We noted that, as depicted in Figures 11 and 12, the binary formation rate becomes negligibly small when the gas friction rate 1/τ drops below the critical value ;0.03ΩK (which depends on the gas drag models).As discussed in Section 6, the rate is not zero even in the limit of 1/τ → 0 since binary formation can be facilitated by GW emission.We can estimate the "gas-free" binary formation rate as follows.When two sBHs have the initial orbital separation K centered around the special value K 0 ; 1.87 or 2.2, they enter the Hill sphere and approach each other with nearly zero relative angular momentum.Since the typical velocity inside the Hill sphere is of order H , we can estimate the effective width ΔK (around K 0 ) within which the two sBHs enter the GW-dominated regime as 1 2 .Thus, the dimensionless binary formation rate in the gas-free limit is This is consistent with the scaling result obtained by Li et al. (2022a).Figure 2).Right: same as the left panel, except using the shearing box calculation.Note that the jaggedness of the curves (especially in the left panel) results from the finite number of grid points in the (K-τ −1 ) parameter space of our calculations and the sensitive dependence of binary formation on the value of K, especially in the borderline region (see, e.g., Figure 2).

Summary
In this paper, we have studied the formation of sBH binaries in AGN disks due to the frictional forces from the gas.We have focused on the "SMBH+2sBH" systems, where the sBHs are initialized on circular and coplanar orbits separated by K times the Hill radius R H .When K is small (2.5), the sBHs quickly experience a close encounter, during which they may capture each other and form a tightly bound binary due to the energy dissipation caused by gas.The key results of this paper are summarized as follows.
(i) Using N-body simulations with a linear frictional force, we systematically investigate the conditions of binary formation and the orientation distributions (prograde versus retrograde) of the newly formed binaries (Section 3).We show that the formation probability depends strongly on the initial orbital separation K and the frictional damping timescale τ (Figure 2).In particular, retrograde binaries begin to form when τ −1  0.03Ω K , and prograde binaries begin to form when τ −1  0.1Ω K .The results depend weakly on the initial orbital phase difference Δf 0 (Figure 3) and are qualitatively the same for sBHs with different masses m 1,2 /M (Figure 4).These Nbody results can be reproduced by shearing box calculations (Section 5), except that N-body calculations allow for binary formation during two or more close encounters.Multiple features of the distribution of the binary orientations in the K − τ −1 parameter space are identified (regions I, IIa, IIb, and IIc in Figure 8); these are associated with four groups of qualitatively different trajectories in the gas-free regime (see Figure 9).
(ii) We rerun some simulations using the GDF forces (Section 4).The results are similar for the two gas force models when τ (of the linear friction; see Equation (4)) and τ eff (of GDF; see Equation (13)) are similar.Based on Equation (7) and the AGN disk models in the literature, we estimate the possible disk locations where BBHs may form (Figure 7).
(iii) We also rerun our shearing box calculations that include the possibility of binary formation via GW emission at very small binary separation (r  r GW = R H ; see Equation (23)).These GW-assisted binary formations occupy specific domains in the K − τ −1 parameter space (see Figure 10), depending on the number of pericenter passages the sBHs undergo before GW emission dominates, which leads to binary capture.The binary orientations can be either prograde or retrograde.
(iv) Based on our numerical results for the formation probability as a function of K and τ, we obtain the rate of binary formation in AGN disks as a function of the sBH number density n and the gas damping timescale τ (or τ 0 ; Section 7).The results show that the binary formation rate increases approximately linearly with 1/τ (Figures 11 and 12) when 1/τ is above a threshold (;0.03ΩK ); below this threshold, GW-assisted capture dominates.The ratio between the prograde and retrograde binaries depends on τ or τ 0 ; when the gas friction is weak (small 1/τ), only retrograde binaries may be assembled, and when the gas friction is strong, prograde binaries may become equally likely as retrograde binaries.

Discussion
We mention several caveats of our study.First, all our calculations assume that the two sBHs initially have circular orbits and zero inclination relative to the disk (so that their motion remains coplanar).In reality, the sBHs may have a small mutual inclination and orbital eccentricities, and these may affect some aspects of our results (see Rom et al. 2023).Our preliminary studies suggest that an inclination of ∼1°has a negligible effect on the gas-dominated formation of binaries.This is because, to form a binary, the gas damping timescale τ needs to be less than about 5P K , which is considerably shorter than the synodic period between the sBHs.This implies that the gas drag is able to damp most of the out-of-plane motion of the sBHs before the BH-BH encounter, so that the formation outcomes are the same as in the coplanar cases.However, binary formation through GW emission may be sensitive to the residual inclinations, because r GW /R H is very small.Future studies need to examine this effect in detail (see Li et al. 2022a, for the gas-free scenario).
Second, we have only studied two-body encounters in this paper.Multiple sBHs may reside in the same region of the AGN disk (e.g., due to migration trap).Interactions between three or more sBHs may lead to new features of binary BH formation.
Third, to link our results with the observed GW events, one must also further study the follow-up evolution of the BBHs after their formation.For example, the orientation of the BBHs may still change after their formation (e.g., through BH-gas interaction or binary-single scatterings).Recent hydrodynamical simulations by Dittmann et al. (2023) suggest that binaries that are not perfectly aligned with the AGN disk can be driven into alignment during their long-term evolution.
Finally, the most important caveat of our study is our treatment of the gas effect.Although the two gas friction models examined in this paper yield qualitatively similar results in terms of the binary formation outcomes in the K − τ −1 parameter space, both models are highly simplified compared to reality.Hydrodynamical simulations show that when two sBHs undergo close encounters, the dominant interaction involves the collision of the minidisks around each BH and the postcollision drag (Li et al. 2023;Whitehead et al. 2023).More studies using hydrodynamics simulations are needed to obtain a more definitive understanding of the gas effects on binary BH formation in AGN disks.

Figure 1 .
Figure 1.Selected results of the N-body simulations with the linear frictional force (model 1; Equation (4)).Each row represents one simulation: from left to right, the panels show the relative trajectory of the sBHs in the corotating frame of m 1 (with X and Y being the radial and azimuthal components of r 2 − r 1 ), the specific relative angular momentum of the sBHs ℓ rel = l rel • z (Equation (11)), and the osculating orbital elements a rel (Equation (10)) and ( ) = e ℓ G m a 1 rel rel 2 12 rel .The simulations adopt /t = W 1 0.3 K , m 1 = 10 −6 M, m 2 = 0.5 × 10 −6 M, and Δf 0 = π/2, and the values of the initial K are shown in the left panels.Prograde binaries (with ℓ rel > 0) are formed in the first, third, and fifth rows; retrograde binary (with ℓ rel < 0) are formed in the fourth row; and no binary formation is observed in the second row.

Figure 2 .
Figure 2. Outcomes of N-body simulations with different initial orbital separations K and frictional damping rate τ −1 (in units of Ω K , the initial orbital frequency of m 1 ).All simulations adopt m 1 = 10 −6 M, m 2 = 0.5 × 10 −6 M, Δf 0 = π/2, and the model 1 gas force (Equation (4)).Each dot represents an integration.The black, blue, and red dots represent the distribution of simulations that do not form binaries, form prograde binaries, and form retrograde binaries, respectively.

Figure 3 .
Figure 3. Same as Figure 2, except for simulations with different initial angular separations Δf 0 .We set K = 2 in all simulations.

Figure 5 .
Figure 5. Four examples of retrograde/prograde BH binaries that form through consecutive close encounters.The left panels show the relative trajectories in the corotating frame of m 1 (with X and Y being the radial and azimuthal component of r 2 − r 1 ) during the close encounters, while the right panels show the time evolution of the BH orbital radii r 1,2 (distance to the SMBH).All cases adopt 1/τ = 0.3Ω 1 , Δf 0 = π/2, m 1 = 10 −6 M, m 2 = 0.5 × 10 −6 M, and the values of the initial K as shown in the left panels.

Figure 6 .
Figure6.Same as Figure2, except for simulations with the canonical GDF (model 2; Equation (5)).Since τ 0 is mass-dependent, the two vertical axes show timescales of the gas drag acting on the inner sBH, m 1 .The second vertical axis on the right indicates the effective dissipation timescale τ eff (Equation (13)) to compare with the model 1 result shown in Figure2.

Figure 8 .
Figure 8. Outcomes of our shearing box model (Equations (14) and (15)) with different initial orbital separations K and frictional damping timescales t.All calculations adopt y 0 = 10.The black, blue, and red dots represent the systems that do not form binaries, form prograde binaries, and form retrograde binaries, respectively.The distribution of the outcomes can be divided into regions I, IIa, IIb, and IIc.BHs evolving from regions I, IIa, and IIc eventually become prograde binaries, while those from region IIb form retrograde binaries.

Figure 10 .
Figure10.Same as Figure8except that here we include binary formation by GW emission in addition to the gas drag.The blue and red dots represent binaries formed without triggering GW emission, i.e., the only dissipative process being gas drag.The deep green and light green dots represent prograde and retrograde binaries formed due to GW emission, respectively.The labels ending with "PP" on the green dots indicate during which episode of pericenter passage the system reaches |r 1 − r 2 | < r GW = 2 × 10 −5 R H and forms a bound binary due to GW emission.

Figure 11 .
Figure11.Left: prograde (blue), retrograde (red), and total (green) binary formations calculated based on Equation (26) and the model 1 results (linear gas friction; Figure2).Right: same as the left panel, except using the shearing box calculation.Note that the jaggedness of the curves (especially in the left panel) results from the finite number of grid points in the (K-τ −1 ) parameter space of our calculations and the sensitive dependence of binary formation on the value of K, especially in the borderline region (see, e.g., Figure2).