Stellar/BH Population in AGN Disks: Direct Binary Formation from Capture Objects in Nuclei Clusters

The Active Galatic Nuclei(AGN) disk has been proposed as a potential channel for the merger of binary black holes. The population of massive stars and black holes in AGN disks captured from the nuclei cluster plays a crucial role in determining the efficiency of binary formation and final merger rate within the AGN disks. In this paper, we investigate the capture process using analytical and numerical approaches. We discover a new constant integral of motion for one object's capture process. Applying this result to the whole population of the nuclei cluster captured by the AGN disk, we find that the population of captured objects depends on the angular density and eccentricity distribution of the nuclei clusters and is effectively independent of the radial density profile of the nuclei cluster and disk models. An isotropic nuclei cluster with thermal eccentricity distribution predicts a captured profile $d N/d r \propto r^{-1/4}$. The captured objects are found to be dynamically crowded within the disk. Direct binary formation right after the capture would be promising, especially for stars. The conventional migration traps that help pile up single objects in AGN disks for black hole mergers might not be required.


INTRODUCTION
Black hole (BH) mergers within the active galactic nucleus (AGN) disks have garnered considerable attention as a mechanism potentially explaining the existence of massive merging BHs (Vokrouhlicky & Karas 1998;Cuadra et al. 2009;McKernan et al. 2012McKernan et al. , 2014;;Bartos et al. 2017;McKernan et al. 2018;Hoang et al. 2018;Secunda et al. 2019;Yang et al. 2019b,a;McKernan et al. 2020;Tagawa et al. 2020;Li et al. 2022a;Bhaskar et al. 2022).Nuclear stellar clusters (NSCs), as the densest environments of stars and BHs, coexist with most of the supermassive black holes (SMBHs) (Paumard et al. 2006;Merritt 2010;Genzel et al. 2010;Kormendy & Ho 2013).For binaries around the SMBH, the typical binary semi-major axis is much larger than the critical size that gravitational wave radiation can drive the binary to merge within the Hubble time.However, the secular perturbations from the SMBH can excite the eccentricity of binaries to near unity, resulting in high gravitational wave radiation efficiency (Wen 2003;Antonini et al. 2015;Antonini & Rasio 2016;Stephan et al. 2016;VanLandingham et al. 2016;Petrovich & Antonini 2017;Liu & Lai 2018;Hoang et al. 2018;Liu et al. 2019a,b;Bhaskar et al. 2022).These eccentric binaries can shrink their semi-major axis by gravitational wave radiation much faster than circular binaries and could eventually merge within the Hubble time.However, due to the relatively large velocity dispersion and the deep gravitational potential of the SMBH, it is not easy to form binaries within the NSC.Unlike the galactic field with a binary fraction of approximately 50%, the binary fraction in star clusters is believed to be less than ★ yihan.wang@unlv.edu10% (Ivanova et al. 2005;Hurley et al. 2007;Sollima et al. 2007).Moreover, most of the binaries in NSCs are located in the outskirts of the NSC, far away from the SMBH, where the secular perturbations from the SMBH are weak.The timescale for the SMBH to excite the eccentricity of binaries located at larger distance to the SMBH is very long.Therefore, a direct binary BH merger from the secular effects of the SMBH may not be efficient.
In the active phase of a galactic nucleus, gas is funneled into the nuclear region, leading to the assembly of accretion disks.The accretion process associated with high luminosity makes the SMBHs visible as active galactic nuclei.Once an AGN disk is formed, some of the stars and BHs in NSC may eventually be trapped within the disk through star/BH-disk interactions (Artymowicz et al. 1993;Rauch 1995;Kennedy et al. 2016;Panamarev et al. 2018;Macleod & Lin 2020;Fabj et al. 2020;Davies & Lin 2020;Generozov & Perets 2023;Nasim et al. 2023).The moving trapped stars/BHs excite density waves in the in AGN disk.The Lindblad resonance exerts torques on those embedded objects, pushing them to migrate within the disk (Tanaka et al. 2002;Baruteau & Lin 2010;Baruteau et al. 2011).Previous studies (Bellovary et al. 2016;Peng & Chen 2021) have suggested that migration traps may be required to accumulate single BHs such that there's a dense region for binary to be formed through three-body/multi-body scatterings (Hills 1975;Aarseth & Heggie 1976;Heggie et al. 1996;Leigh et al. 2018;Zevin et al. 2019) or gravitational wave captures (O'Leary et al. 2009;Kocsis & Levin 2012;Gondán et al. 2018;Samsing et al. 2020;Li et al. 2022b).However, the hypothetical existence of the migration traps ignores multi-body interactions during the migration process, stochastic dif-fusion due to chaotic torque in turbulent disks (Wu et al. 2024), and relies on special idealized disk models.
In addition to migration traps, Li et al. (2023); Rowan et al. (2022); Rozner et al. (2023); DeLaurentiis et al. (2023) also investigated direct binary formation in AGN using gas dissipation with arbitrary initial BH populations.However, accurately estimating the rate of BH binary formation in AGN disks has proven challenging due to uncertainties in the density profile, initial mass function, binary fraction in NSC, and poorly constrained AGN disk models.
It is particularly challenging to accurately predict the BH binary merger rate within the AGN disk due to the unclear path from binary formation to the mergers.In the case of disk-embedded binaries, there are two primary mechanisms for reducing the semi-major axis to the regime of gravitational wave radiation: dynamical encounters (Tagawa et al. 2021a,b;Wang et al. 2021b;Samsing et al. 2022) and gas dissipation (Baruteau et al. 2011;Li et al. 2022a;Li & Lai 2023a,b;Kaaz et al. 2023).
Regarding dynamical encounters, if these embedded binaries exist in a denser dynamical environment, close interactions between these binaries and other single or binary systems can occur frequently.When the binaries are hard, characterized by significant environmental velocity dispersion, the scattering process statistically hardens the binary system.Notably, in this impulse regime, the scatterings tend to increase the eccentricity of the binaries, thereby accelerating the hardening process through enhanced gravitational wave radiation.However, this process is sensitive to the local populations of BHs and stars.The presence of hard binaries and an adequate number density are required for the hardening process to take effect.
In the case of gas-assisted mergers, the hardening process is more dependent on the binary orientation within the AGN disks.Retrograde binaries, with a greater relative velocity compared to the local disk velocity, efficiently dissipate orbital energy (Li & Lai 2023b).However, for prograde binaries, there is no definitive conclusion on whether the gas can sufficiently reduce the binary's semi-major axis to reach the regime of gravitational wave radiation.Furthermore, the population of prograde/retrograde binaries, which relies on the initial conditions of binary formation, remains poorly constrained in current research.
All these subsequent steps (migration, binary formation, hardening) leading to BH merger depend on the capture process of stars and BHs by the AGN disks.When the relative velocity greatly exceeds the surface escape velocity, e.g. for stars on most orbits intersecting an AGN disk, the dominant force is aerodynamic drag.Conversely, when the relative velocity is considerably lower than the surface escape velocity of the object, e.g. for black holes and neutron stars on practically all orbits intersecting an AGN disk, the dominant force is gas dynamical friction.In this study, we investigate the capture of stars and BHs from the nuclei cluster by AGN disks, focusing on estimating the timescales associated with semi-major axis damping, inclination decrease, and eccentricity excitation/damping.By analyzing the dependence of these timescales on the initial orbital properties and their relative configuration with respect to the AGN disk, we aim to accurately predict the captured star/BH population in AGN disks.
Our paper is organized as follows.In Section 2, we investigate the disk-star/BH interactions, derive the equation of motion for this capture process, and obtain the timescales of semi-major axis, eccentricity, and inclination evolutions.In Section 3, we perform N-body simulations to test the physics picture we obtained from Section 2. In Section 4, based on the timescales and derived integral of motion, we calculate the population of the captured stars/BHs via Monte Carlo simulations and provide analytical expressions for stars/BHs population in AGN disks.

DISK-STAR/BH INTERACTION
Figure 1 shows the schematics of the disk-star/BH interactions.The celestial body in orbit around the supermassive BH traverses the disk's mid-plane at two specific points known as the ascending node and descending node.Given the axis symmetry of the system, it is possible to establish Cartesian coordinates in such a way that the ascending node coincides with the positive x-axis.Consequently, the positional vector of the ascending and descending interacting nodes along the x-axis can be mathematically expressed as follows: where  is the argument of periapsis,  = (1 −  2 ) is the semi-latus rectum with semi-major axis  and eccentricity .The corresponding velocities of the object (v * )and the Keplerian rotating disk (v  ) at the crossing point are where  tot =  SMBH +  is the total mass of the supermassive black hole and crossing object, and  is the argument of the periapsis.The relative velocity between the object and the disk can be expressed as And the magnitude of the relative velocity is where we use Upon the celestial object crosses the disk, it becomes subject to a drag force.The magnitude of this force depends on the ratio between the relative velocity of the object with respect to the disk and the surface escape velocity of the object.Two distinct mechanisms can dominate the drag force: aerodynamic drag force and gas dynamical friction.
When the relative velocity greatly exceeds the surface escape velocity of the object, the dominant force is the aerodynamic drag force, which is directly proportional to the square of the relative velocity ( 2 rel ).Conversely, when the relative velocity is considerably lower than the surface escape velocity of the object, the dominant force is gas dynamical friction, which is inversely proportional to the square of the relative velocity ( −2 rel ) (Ostriker 1999).If the aerodynamic drag force dominates the overall drag force, it can be described by the following equation: where   is the gas density of the disk, v rel is the relative velocity between the crossing object and the disk, and  is the effective radius of the crossing object.The effective radius is defined as the maximum value between the physical radius and gravitational radius / 2 rel .If the dynamical friction dominates the overall drag force, the drag force can be described as, where I is a function of Mach number M.
with Coulomb logarithm Λ = log( max / eff ), where  max is the typical size of the medium.In the subsonic regime where M < 1, F dyn is asymptotic to M/3 and in the supersonic regime where M > 1, F dyn is asymptotic to M −21 .The criteria remain essentially the same for cases where the gravitational radius dominates over the physical radius or when gas dynamical friction dominates over aerodynamic drag (considering the 1 Note that in the subsonic regime, the gas dynamical friction is typically much smaller than the dynamical friction in collisionless medium, e.g.background stars, where 2 with velocity dispersion of the collisionless medium .physical radius), where For main sequence stars, BHs, and neutron stars in circular orbits, the criteria writes as, For compact objects, such as black holes or neutron stars, the gravitational radius generally dominates over the physical radius.In addition, gas dynamical friction tends to dominate over aerodynamic drag when considering the physical radius (for BHs, the event horizon).On the other hand, for main sequence stars, the physical radius typically outweighs the gravitational radius in the majority of parameter space, except in cases where the inclination is sufficiently low.Consequently, there exists a critical inclination for main sequence stars, below which the gravitational radius becomes significant compared to the physical radius, and gas dynamical friction becomes dominant over aerodynamic drag, Also see (Grishin & Perets 2015) for a comparison between gas drag and gas dynamical friction for planetesimals moving in protoplanet disk.

Aerodynamic drag
If the disk is assumed to be thin (h≪1), and the drag force is dominant by the aerodynamic drag with geometry cross-section, the specific momentum and specific angular momentum change per disk cross are, where Δ ∼ / rel,⊥ and  =   2 * Σ   rel  rel,z .The specific energy, specific angular momentum, and the relative velocity are, where  = √︁   tot  and  is the orbit inclination.Based on Equation 20, the specific angular momentum and energy change per disk cross are, ΔL =  (0, sin , 1 Assuming that the relative angular momentum change per disk crossing, denoted as , is small ( ≪ 1), we can use Equation 23and the relationship cos  =   / to derive the linear equation of motion for the angular momentum and inclination, where We can also get d Σ at the two crossing points  ± is assumed to be nearly constant.This assumption holds true when  cos  is not close to unity, indicating that the orbit is not extremely eccentric with cos  ≈ ±1.If  cos  is close to unity, it implies that  − is orders of magnitude smaller than  + , causing the surface density Σ( − ) to be much larger than Σ( + ).Consequently, the dissipation is primarily influenced by the crossing point  − , and terms involving , , and  with  + and  + can be disregarded.

Inclination damping and eccentricity evolution
Equation 27 indicates that d cos /d is always greater than zero, implying that the inclination evolves towards cos  → 1, or in other words, the inclination tends to align with the disk plane, which is consistent with direct 3-D hydrodynamical simulations (Rein 2012;Arzamasskiy et al. 2018;Zhu 2019).For objects in nuclei clusters with inclined orbits, the aerodynamic drag will gradually align the orbits with the disk.If the time is long enough, all objects with different orbital configurations will be captured by the disk in prograde orbits.However, in reality, the lifetime of the AGN disk might be short compared to the capture timescale.Therefore, only a fraction of the objects in the nuclei cluster can be captured.
In Equation 26, we have d 2 /d ∝ cos  − .Therefore, if cos  > , d 2 /d will be positive, indicating the angular momentum will increase, otherwise, it will decrease.From Figure 2, we can observe that  is always smaller than one.The criterion for angular momentum damping is cos  < .If the inclination is low enough and the periapsis/apoapsis are closer to the disk mid-plane ( cos  → ±1), the angular momentum can increase.In this configuration, at the apoapsis, the local disk Keplerian velocity is much higher than the velocity of the orbiter.Thus, the disk can accelerate the orbiter and transfer angular momentum to the orbiter.From Equation 39, due to ᾱ > 0 and γ > 1, the semi-major axis always decays.Therefore, the angular momentum increase indicates fast eccentricity damping.
Note that this angular momentum increase criterion also requires the aerodynamic drag to dominate over the gas dynamical friction.Therefore, the inclination must be larger than the critical inclination indicated by Equation 18.Therefore, the parameter space for angular momentum increase also requires cos  d > .Otherwise, the required inclination for angular momentum increase is below the critical inclination where aerodynamic drag dominates the total drag force.
Regarding the eccentricity evolution, one can prove that Equation 40 is always negative if  ≠ 0. So the eccentricity will always decrease.Figure 3 shows the value of  as a function of cos  for different  and .As shown in the upper panel and bottom panel, for cos  = 0, where the two crossing points have the same distance to the SMBH, or for cos  = ±1, where the periapsis and apoapsis lie within the disk,  is always negative for various initial eccentricities and inclinations.Therefore, the eccentricity will be damped during the capture process.

Timescales
From Equation 26to 40, we can obtain the timescales for semi-lotus rectum, semi-major axis, inclination and eccentricity evolution, For the aerodynamic drag-dominated regime, the timescales for ,  and  evolution are at the order of magnitude of Σ * Σ .For main- where the periapsis and apoapsis lie within the disk.In general, the semi-major damping timescale is much shorter than the inclination damping timescale and the inclination damping timescale is shorter than the eccentricity excitation/damping timescale.
sequence stars, if we use the mass-radius relation we can get So, for massive stars, the timescales for semi-major axis, eccentricity, and inclination evolution are shorter.
Figure 4 shows the timescales of semi-major axis, eccentricity, and inclination evolution in a unit of  0 =  Σ * Σ for different  and .The upper panel shows the case cos  = 0, where the two crossing points on disk mid-plane have the same distance to the central SMBH.
For low eccentricity, we can see the timescale for semi-major axis damping is much shorter than the timescale of inclination damping with high inclinations and is much longer than the timescale of inclination damping with low inclinations.For high eccentricities, the semi-major axis damping timescale is always shorter than the other two timescales.The timescale for eccentricity is always shorter than the inclination damping timescale for both low and high eccentricities.Therefore, for high-inclination orbits, the sem-major axis and eccentricity decrease fast, with nearly unchanged inclination.As the inclination decreases and the orbit becomes circularized, the semimajor axis damping timescale becomes longer and longer, thus the orbit will eventually stall at a certain semi-major axis.
The bottom panel shows another case, where the periapsis and apoapsis lie on the disk mid-plane, the general picture is very similar to the case with cos  = 0, except for the inclination damping timescale is more sensitive to the eccentricity.

Integral of motion
As indicated by Figure 2, the  is nearly constant within one order of magnitude and is close to unity for almost all eccentricities and .Therefore, we can integrate Equation 28 to be with approximation  ∼ 1, this expression can be further simplified to, From the previous discussion we learn that if the capture timescale  I,aero is shorter than the lifetime of the disk, the objects in the nuclei cluster will be eventually captured by the disk.As a result, the final inclination ( f ) will be effectively zero.The final captured angular momentum can be calculated by From the previous discussion, we know that as the inclination is below the critical inclination  d , the dynamical friction takes over the aerodynamical drag.Later, we will show that for gas dynamic friction, this integral of motion holds almost true as well, and the eccentricity will always be circularized once the inclination is low enough.Therefore, the captured objects will effectively have zero eccentricity.In this way, the final captured semi-major axis can be well described as This integral of motion is effectively independent of , and thus is independent of the surface density of the disk.Since  tot =  SMBH +  ∼  SMBH , this constant is also essentially independent of the mass of the crossing objects.Therefore, even if the crossing objects are accreting during the capture process, the integral of motion still holds true.The captured objects will shrink the semi-major axis by a factor of . Therefore, the initial retrograde obiters with high eccentricity will be captured in the innermost region of the disk.For orbiters with an initial inclination larger than the critical inclination where  TDE is the tidal disruption radius of the SMBH on orbits, they will be tidally disrupted by the SMBH.There's also another critical angle   , if the initial inclination is greater than this angle, the stars will be tidally disrupted in a retrograde orbit out of the disk.Following the integral of motion, we can obtain the critical angle The general picture for the stellar capture is shown in Figure 5, where for a given semi-major axis, two critical angles   and  TDE divide the parameter space into two different regimes.For very high inclination orbits, the semi-major axis shrinks very fast by a factor of . This contributes to the TDE by disk captures.
For retrograde orbits not associated with extremely high inclination, their semi-major axis damping would be slower and associated with eccentricity decrease.

Gas dynamical friction
If the dynamical radius  2  2 / 2 rel is larger than the geometry radius  * or the gas dynamical friction dominant over the aerodynamic drag, the angular momentum equation should be where Similar to Section 2.1, we can write the equation of motion of the angular momentum, inclination, semi-major axis and eccentricity, Figure 5. Carton shows the capture of stars with different initial inclinations at a given fixed semi-major axis.High-inclination stars will end up with smaller semi-major axis while low-inclination stars tend to have relatively larger final semi-major axis.Orbits with inclination  smaller than  TDE will eventually be captured by the AGN disk while orbits with inclination  greater than  TDE will be tidally disrupted by the SMBH.The trajectories and critical angles are not to scale, especially the semi-major axis that shrinks fast than the inclination damping and eccentricity evolution.

Capture, eccentricity damping, and excitation
Equation 58 indicates that d cos /d is always positive.Therefore, similar to aerodynamic drag, the gas dynamical friction decreases the inclination of the orbiters.For Equation 57, if cos  −  is negative, the angular momentum will decrease while if cos  −  is positive, the angular momentum will increase.The reason for the angular momentum increase/decrease is the same as the aerodynamic drag that has been discussed.Figure 6 shows  as functions of  and .Similar to Figure 2,  is close to unity but exhibits a more pronounced curvature as  increases.
For gas dynamical friction, the criterion for eccentricity damping is It's eas to calculate the critical inclination  e for BH eccentricity damping, cos Figure 7 illustrates the  as a function of inclination for different  and .Similar to Figure 3, the upper panel/bottom panel shows the case for cos  = 0 and ±1, respectively.For cos  = 0, similar to the aerodynamical drag, the eccentricity always decays, even if the orbit is purely retrograde.As cos  e becomes larger than -1, retrograde orbits with inclination cos  < cos  e will undergo eccentricity exci- tation.For eccentricity damping to appear, we need cos  e > −1, such that any inclination  <  e will undergo eccentricity damping.The criterion for cos  e > −1 is | cos | ∼> 0.1 for various eccentricities.Therefore, there's only a small region in the parameter space where eccentricity damping never appears.

Timescales
Similar to Section 2.1, we can get the semi-latus rectum, inclination, semi-major axis, and eccentricity evolution timescales, Figure 8 shows the timescale of semi-major axis, eccentricity, and Σ   2 .The upper panel shows the case for cos  = 0. Similar to aerodynamical drag, the timescale for the semi-major axis is much shorter than the timescales of inclination damping at high inclination and becomes longer as the inclination decreases.The timescales for  and  are insensitive to eccentricity.The bottom panel shows the case for cos  = ±1.For low eccentricity orbits that start from high inclination as indicated by blue lines, the eccentricity excitation timescale is much shorter than the semi-major axis/inclination damping.As the eccentricity becomes higher, the timescale for inclination damping becomes shorter, then the orbit inclination quickly damps, until the orbit becomes prograde, once the inclination is below the critical inclination given by Equation 69, the eccentricity starts to decay.As eccentricity decreases, the timescale for inclination damping becomes longer.The semi-major axis damping takes over the inclination damping.As the inclination decreases further, the semi-major axis decay stops.Then the inclination and eccentricity timescales become shorter again at a very low inclination.A fully circularized orbit will be obtained at the end of the capture process.

Integral of motion
For Equation 61, similar to Section 2.1, we can obtain the integral of motion, assumed that  is nearly constant within one order of magnitude as shown in Figure 6.The difference is that for the same cos  and , the value of  is slightly smaller than .
The general picture for BHs is very similar to the stars as shown in Figure 5.However, there's no parameter space for TDE and the eccentricity excitation timescale is much shorter than the stars.Therefore, the retrograde BHs will undergo significant eccentricity excitation.Meanwhile, the inclination-damping timescale for the BH is much more sensitive to the eccentricity, unlike the star the inclinationdamping timescale is nearly constant.High eccentricity orbits decrease inclination much faster than circular orbits.Therefore, all the BHs will be eventually captured by the disk in prograde orbits.Since all the orbits will enter the fast eccentricity damping regime, they will eventually be fully circularized.
Figure 9 illustrates the BH capture process.For initially retrograde BHs with an inclination higher than   , their eccentricity becomes excited as the inclination dampens to   .Once their orbital inclination is below   , their eccentricity starts to decrease.

N-BODY EXAMPLES
To validate the calculations presented in Section 2, we conducted N-body simulations using the software tool SpaceHub(Wang et al. 2021a).Two sets of simulations were performed to examine the capture process involving a black hole (BH) and a main sequence star, where gas dynamical friction and aerodynamic drag dominate, respectively.In both cases, the mass of the supermassive black hole (SMBH) was set to 10 8  ⊙ , and the mass of the star or BH was set to 30 ⊙ .The initial semi-major axis was fixed at 0.1 pc for stars and 0.01 pc for BHs, and initial eccentricity was set to thermal average value ∼ 0.672 .The initial inclinations were evenly distributed between 5 and 175 degrees, with an interval of 17 degrees.Each case was simulated with two different values of cos : 0 and ±1.
For the star, the force was implemented according to Equation 12, while for the BH, the drag force was implemented as described in Equation 13.The Coulomb logarithm Λ was kept constant at 3 for gas dynamical friction.

Disk models
We adopted two different disk models for the capture process, a thin -disk model that can be well described by simple equations and a thick Sirko-Goodman (SG) disk model that has larger surface density around the active capture radius.

𝛼-disk
The  disk model can be characterized by three parameters: the accretion rate efficiency ( d ), the Toomre parameter ( d ), and the viscosity parameter ( d ).The accretion rate can be approximated as follows: where,  represents the accretion rate,  Edd denotes the Eddington accretion rate,  SMHB represents the mass of the SMBH,  Edd is the Eddington luminosity, and  d is a constant with a value of 0.1.The disk surface density Σ can be described by, where,  represents the distance from the SMBH,  is the scale height, and Ω d denotes the orbital frequency given by √︁   SMBH / 3 .The parameters  d ,  d , and  d are assumed to be constant values set to 1.

Sirko-Goodman (SG) disk
The solution of the SG model (Sirko & Goodman 2003) can only be solved numerically.There are no simple equations that can be used directly to describe the SG disk profile.
Figure 10 shows the surface density and specific scale height of the -disk and SG disk around a 10 8  ⊙ SMBH. the surface density and specific scale height of the SG model with  > 1.2 × 10 3   is roughly 3 times of the surface density of an -disk model with  = 1.In the inner region, the surface density of the -disk is much higher than the SG disk since the inner disk in the SG model is gravitationally stable.The surface density and specific scale height of the SG disk in the outer disk regime are similar to an -disk with  = 0.01.

Stellar trajectories
Figure 11 depicts the trajectories of a 30 solar mass star orbiting an SMBH with a fixed initial thermal average eccentricity of 2/3 ∼ 0.67 and various initial inclinations.All orbits are assumed to have cos  = ±1 initially.
The upper middle panel of the figure demonstrates that for orbits with high inclinations, the semi-major axis experiences faster decay, as described by Equation 42and Figure 4.The inclination damping timescale is much longer than the semi-major axis damping timescale, as observed in the bottom middle panel.Consequently, stellar orbits shrink rapidly while inclination damping occurs at a slower rate.
As the inclination decreases, the semi-major axis damping timescale increases significantly until it surpasses the inclination damping timescale.This results in the stalling of the semi-major axis decay and a pronounced decrease in inclination, as shown in the bottom panel of Figure 4.
The upper right panel displays the eccentricity evolution.In agreement with Equation 40 and the bottom panel of Figure 4, both prograde and retrograde orbits experience eccentricity damping.As illustrated in the bottom right panel, orbits with very high inclinations exhibit a significantly shorter timescale for eccentricity excitation compared to the inclination-damping timescale.However, as the Figure 9. Carton shows the capture of BHs with different initial inclinations at a given fixed semi-major axis.For orbiters with  >  e , the orbital eccentricity will be excited.However, as the inclination decreases, the orbiter will enter a regime where eccentricity damping is fast.The trajectories and critical angles are not to scale, especially the semi-major axis that shrinks fast than the inclination damping and eccentricity evolution.
semi-major axis rapidly decreases, the aerodynamic drag becomes very efficient due to the high disk surface density.The quasi-static assumption that the orbital parameter changes per period are small, which we adopted in the derivation, is no longer true.Thus, for retrograde orbits with high inclinations, we observe some unpredicted eccentricity excitations.
The upper left panel presents the inclination evolution in the presence of semi-major axis damping.As the semi-major axis decreases, the disk's surface density increases substantially, leading to a much shorter inclination-damping timescale.Since high-inclination orbits have shorter semi-major axis damping timescales, orbits with high inclinations are captured more rapidly by the disk.
The trajectories with other initial eccentricities and cos  are very similar to Figure 11 and can be well described by timescales in Section 2.1, Figure 4 and Figure 5.

BH trajectories
Figure 12 shows the trajectories of 30  ⊙ BHs starting with different inclinations.For large initial inclinations, because the eccentricity timescale is much shorter than the semi-major axis and inclination damping timescale as indicated by Figure 8, and the  is positive, therefore, the eccentricity of highly inclined orbits grows fast.As the eccentricity becomes high, the inclination damping timescale of those eccentric orbits becomes very short, leading to fast inclination decreases.Once the orbital inclination is below the critical angle  e , the eccentricity starts to decrease.For retrograde BHs, they all will undergo thus an eccentricity excitation and damping process during the interaction with the disk.Unlike the stars, whose inclination damping timescale is insensitive to the eccentricity, the inclination damping timescale for BH is very sensitive to the eccentricity.Therefore for retrograde BHs, the capture process can be divided into three stages: fast eccentricity growth, quick inclination damping and slow eccentricity damping.The semi-major axis keeps shrinking along the three stages.Since the eccentricity damping may not be as efficient as inclination damping, there might be some eccentricity residual on captured BHs.
For initially prograde BHs, the picture is very similar to the stars, except the inclination and eccentricity evolution timescale is relatively shorter than the semi-major axis damping timescale.Due to this reason and the residual eccentricity in the integral of motion (1 −  2 ) cos 4 (/2), the semi-major axis of captured BHs may not shrink as much as stars.

Test the integral of motion
To verify the integral of motion we obtained in Section 2.1.3,we set up Monte Carlo (MC) simulations with different initial orbital conditions and tracked the orbital parameters of those orbits during the disk capture process.The mass of the black hole (BH) and star is kept constant at 50 ⊙ .The semi-major axes are drawn from a distribution such that the density profile of the nuclear star cluster (NSC) follows ∝  −1.5 from 1  to 0.1 pc.The eccentricities are drawn from a thermal distribution with a probability density function () = 2, and the inclinations are drawn from a distribution where cos  is uniformly distributed between -1 and 1.We use the -disk model for the MC simulations.The MC simulations are terminated once the inclination of the orbiter is below the specific scale height of the disk.Figure 13 shows the integral of motion at the beginning and end of the MC simulations.The upper panel shows the results for stars, while the bottom panel displays the results for black holes (BHs).As depicted in Figures 2 and 6,  has a much larger range than .Therefore, according to Equation 47 and 74, it is expected that  cos 2 (/2) will deviate more for the BH.One should note that this integral of motion only holds true for the capture process.Once the objects become fully embedded within the disk, this integral of motion is no longer a constant value.

POPULATIONS OF CAPTURED OBJECTS
Due to the significant observational challenges, the profile, compact object fraction, and binary fraction in the nuclear star cluster (NSC) are poorly constrained.Additionally, the capture process in AGN disks can be intrinsically complex with numerous uncertainties due to the variety of AGN disk models.Consequently, establishing the population of captured objects in AGN disks can be a challenging task that highly depends on the adopted NSC and disk models.
However, in this section, we will demonstrate that certain properties of the captured objects remain insensitive to the NSC and disk models.This suggests that it is possible to establish a robust population of captured objects in AGN disks.
As discussed earlier, due to the short timescale of eccentricity damping, both stars and black holes (BHs) will be captured in nearly circular orbits.Using the integral of motion, the final captured semimajor axis can be well described as  f =  0 (1 −  2 0 ) cos 4 ( 0 /2).Furthermore, since the inclination damping timescale,   , is insensitive to the inclination, eccentricity, and argument of periapsis of the orbit, we can obtain the distribution of captured objects in AGN disks at any given time for a given disk model and NSC profile.

Density power-law index for captured objects
We employ the same disk models as in Section 3.1, namely an -disk with reasonable disk parameters  d = 1,  d = 1, and  d = 0.1, as well as a SG disk model.For the NSC model, we consider two different masses: 1  ⊙ and 50  ⊙ for stars, and 5  ⊙ and 50  ⊙ for black holes (BHs).The number density/density profile of the NSC follows a power-law distribution, where the number density is given by (Merritt 2013): Here,  represents the mass of the star/BH in the NSC, and  m corresponds to the gravitational influence radius of the supermassive black hole (SMBH).We adopt three power-law index  NSC , 1, 1.5 and 2. Within this radius, the total enclosed masss is 2 SMBH and the dynamics of orbits is predominantly influenced by the SMBH's gravity, allowing us to neglect the effects of the dark matter halo and other stars/BHs.The value of  m can be determined using the following equation: where  NSC represents the velocity dispersion of the NSC, given by (Kormendy & Ho 2013): To obtain the population of captured objects in AGN disks at any given time , we assume that objects in the NSC will be captured by the disk after, the inclination damping timescale ( =  I ).Once objects are captured by the disk, we populate the captured objects at the disk mid-plane.Based on the integral of motion during the capture process, we can calculate the final semi-major axis from  f =  0 (1 −  2 0 ) cos 4 ( 0 /2). Figure 14 depicts the 1-D number density profiles of captured 50  ⊙ stars at different times for various initial NSC profiles and disk models.The upper panels present the number density profiles of captured disk stars for the -disk, while the lower panels display the profiles for the SG disk.Each column corresponds to different NSC profiles.Surprisingly, the number density profiles of the captured disk stars exhibit remarkable consistency across different disk models and NSC profiles.
The captured disk stars follow a density profile proportional to  −1/4 from 10 5 years, which represents the typical lifetime of shortlived AGN disks, to 10 8 years for long-lived AGN disks.The capture process causes an accumulation of stars in the inner region of the disk.The profile of the disk stars is not sensitive to the NSC radial profiles because the final captured semi-major axis is more sensitive to the inclination and eccentricity distribution of the orbits within the NSC, owing to the angular dependence in the integral of motion  cos 2 (/2).As we assume the NSC to be nearly isotropic and the eccentricity is isothermal, the captured disk star profiles exhibit a high level of consistency.Indeed, in the asymptotic regime, if we assume all the orbits can be captured by the disk, the final semimajor axis of the capture objects is Then the probability function for  f can be obtained from where ( 0 ) is the power-law probability function of the initial semimajor axis distribution in NSC and  max =  m is the outer boundary of the NSC.The dashed lines in Figure 14 and 15 indicate that this probability function that follows nearly ∝  −1/4 (asymptotic to ln 2 () 10 6 10 5 10 4 103 10 2 10 1 10 0 10 1 10 2 L 0 cos 2 (I 0 /2) 10 6 in the inner region) is a very good approximation for the final semimajor axis distribution.The captured star/BH number density profile thus depends on the angular distribution of the star/BH in NSC rather than the radial distribution.An isotropic star/BH distribution in NSC gives d/d ∝∼  −1/4 .Figure 15 displays the disk star profile for the 50  ⊙ black holes (BHs).Since the integral of motion for BHs deviates slightly from  cos 2 (/2), the captured BH profile follows a constant distribution.In contrast to star capture, the density profiles for BHs extend to further out of the AGN disks.This indicates that BHs can be captured in the outer region of the disk up to 10 6   .
Figure 16 illustrates the power-law index of the captured disk stars/BHs as a function of time for different disk models and NSC density profiles.In the upper panel, the power-law index of the 3D number density profile for captured disk stars exhibits remarkable consistency and converges to approximately -2.25, indicating / ∝  −1/4 .

Mass filling function of captured objects
Besides the initially embedded objects, the total captured mass can For captured BHs, the total captured mass can also be obtained as Figure 17 shows the captured total mass as a function of time for different disk models and NSC density profiles from our MC simulations.The upper panel shows the cases for stars with  * = 1 and the bottom panel shows the cases for BHs with  BH = 1.As indicated, the mass filling function for stars M * ∝  1− NSC and the mass filling function for BHs M BH ∝ .

Direct binary formation from capture
The formation of binary black holes (BHs) in AGN disks relies on migration traps, where single BHs can be captured and accumulated.In the migration trap, the net torque resulting from Lindblad resonance is zero, allowing BHs to be trapped and accumulate in that region.However, migration traps strongly depend on the specific disk models.Furthermore, the concept of migration traps is based on the single migrator model, where the gravitational effects of other migrators can be ignored.In reality, in an N-body migration system, resonances from other migrators come into play and can significantly influence the migration of individual BHs.Considering that the population of captured stars and BHs can be well described by a power-law distribution regardless of the NSC density profiles and disk models, we can examine how dynamically crowded the captured objects become once they are captured.We employ a quantity called  Hill , which represents the number of objects within the volume of a ring with a radius  and a cross-section radius  Hill =  (   SMBH ) 1/3 .This number indicates how many objects can enter the Hill's radius of a single orbiter over time.If  Hill is larger than two, it suggests that a binary could potentially form through encounters aided by gas dissipation in AGN disks.Conversely, if  Hill is too small, it indicates that the embedded objects are dynamically distant from each other. Hill for an NSC consisting of 1  ⊙ stars (with higher number density), while dashed lines represent 50  ⊙ stars (with lower number density).
Initially,  Hill in the NSC is generally smaller than two, except for the outer layer at a distance of approximately 10 4   from the SMBH.However, the capture process reduces the inclination of orbits at a given  and leads to the accumulation of captured stars in the radial direction, resulting in a dynamically crowded population.In the case of an -like disk with a relatively high surface density in the inner disk region and a relatively low surface density in the outer disk region, along with a not too steep NSC profile ( NSC = 1.5),only NSCs with very high number density consisting of 1  ⊙ stars can be accumulated by AGN disks.If the number density of the NSC is low, only a long-lived AGN disk with a lifetime of approximately 10 8 years can capture stars with  Hill > 2. Notably,  Hill increases significantly as the NSC density profile becomes steeper and insensitive to the disk models.
Figure 19 presents the results for BHs, similar to Figure 18.Since the capture timescale (Equation 71) is much longer for BHs compared to stars (Equation 43), on average, more time is required for the disk to accumulate BHs.For NSC with mild power-law index  NSC = 1.5, only a long-lived AGN disk makes it possible to directly form binary For more concentrated NSCs, the direct BH binary formation is more promising.For a normal lifetime AGN disk, direct BH binary formation occurs around hundreds of   .
Although, in general, captured BHs are not dynamically crowded with  Hill > 2. However, the captured stars are over-dense upon being captured by the AGN disk.There will be frequent binary star formations in the AGN disk around hundreds of   .Those binaries could eventually evolve into binary BHs through stellar evolution.The stellar evolution in the AGN disk is far more complicated that involves non-isotropic accretion, and internal stellar structure, which is beyond the scope of this paper (Cantiello et al. 2021;Jermyn et al. 2022;Ali-Dib & Lin 2023;Huang et al. 2023).However, we believe, in this over-dynamical-dense region, there will be binary BH remaining throughout the stellar evolution.Those binary BH could be the source of merging BHs in AGN disks and bypass the traping phase with lots of uncertainties from N-body migration.

Figure 1 .
Figure1.Schematics of the disk capture.The coordinates system is built in a way such that the ascending node of the star/CO orbit is on the x-axis.

Figure 2 .
Figure 2.  (, ) as a function of cos  for different eccentricities.cos  <  indicates angular momentum damping.

Figure 3 .
Figure 3.  as a function of inclination for different .The upper panel shows the case with cos  = 0 where the two crossing points have the same distance to the SMBH.The bottom panel shows the case with cos  = ±1 where the periapsis and apoapsis lie within the disk.Eccentricity will always decrease for all inclinations.

Figure 4 .
Figure 4. Timescales of semi-major axis, eccentricity and inclination evolution in a unit of  0 = Σ * /Σ for stars.The upper panel shows the case for cos  = 0 where the two crossing points have the same distance to the SMBH.The bottom panel shows the case with cos  = ±1 where the periapsis and apoapsis lie within the disk.In general, the semi-major damping timescale is much shorter than the inclination damping timescale and the inclination damping timescale is shorter than the eccentricity excitation/damping timescale.

Figure 6 .
Figure 6. (, ) as a function of cos  for different eccentricities.cos  <  indicates angular momentum damping.

Figure 7 .
Figure 7. Similar to Figure 3, but for gas dynamical friction.Unlike the aerodynamical drag, there is parameter space for eccentricity excitation.The vertical dashed lines show the critical inclination for eccentricity excitation/damping indicated by Equation 69.

Figure 8 .
Figure 8. Timescales of semi-major axis, eccentricity and inclination evolution in a unit of  0 =   2 tot Σ   2 for BHs.

Figure 10 .
Figure 10.Surface density and specific scale height of -disk and SG disk.

Figure 11 .
Figure 11.Trajectories of stars for different initial inclinations.Upper panels: inclination, semi-major axis and eccentricity as a function of time.Bottom panels:  as a function of ,  as a function of  and  as a function of .Different colors indicate different initial inclinations that can be seen in the upper left panel.

Figure 13 .
Figure 13.The initial value verses the final value of the integral of motion for the numerical simulations.The upper panel shows the cases for star and the bottom panel shows the cases for BH.
2 (cos ) () () 2 ddd cos (85)    where  * is the mass fraction of the star in NSC, () = 2 is the probability function for eccentricity distribution, (cos ) = 1/2 is the angular distribution function and () is the probability function for density distribution in NSC that can be obtained by multiplying  to Equation80.For captured stars, the total captured mass by the AGN disk can be approximated as M * () ∼ 2 SMBH  * and Σ  are the orbital period and surface density at the gravitational influence radius   , respectively 3 .The captured stars distribute between  TDE and  max, * with a power-law distribution d/d ∝  −1/4 , where  max, * ∼   2 SMBH  * indicates all stars within   (with the total mass of 2 SMBH ) are captured by the disk.The corresponding time  =   Σ * /Σ  is exactly the capture timescale for a circular orbit at   , the final captured orbit of the NSC.

Figure 14 .Figure 15 .
Figure 14.d /d as a function of  for the initial NSC consists of 50  ⊙ stars and captured stars in the disk at different times.The upper panels show the results for -disk (with  = 0.1) capture and the bottom panels show the results for SG disk capture.Different columns indicate different initial NSC density profiles.The NSCs are set to be isotropic with number density profile dn/dr ∝  − NSC .Regardless of the NSC density profile and disk model, the density profile of the captured star follows d /d ∝  −1/4 .

Figure 18 Figure 16 .
Figure18depicts  Hill as a function of  for different disk models and NSC density profiles around a 10 8  ⊙ SMBH.The upper panels show the results for an -disk with  = 0.1, while the lower panels display the results for an SG disk with  = 0.01.Solid lines represent

Figure 17 .
Figure 17.Total mass of captured objects as a function of time for different disk models and NSC density profile.The upper panel shows the results for the captured star and the bottom panel shows the results for captured BHs