Formation and Evolution of Compact Object Binaries in AGN Disks

The astrophysical origin of gravitational wave (GW) events discovered by LIGO/VIRGO remains an outstanding puzzle. In active galactic nuclei (AGN), compact-object binaries form, evolve, and interact with a dense star cluster and a gas disk. While most binaries are soft and get disrupted by binary-single interactions, an important question is whether and how binaries can merge in these environments. To address this question, we have performed one-dimensional $N$-body simulations combined with a semi-analytical model which includes the formation, disruption, and evolution of binaries self-consistently. We point out that binaries can form in single-single interactions by the dissipation of kinetic energy in a gaseous medium. This"gas capture"binary formation channel contributes up to $97\,\%$ of gas-driven mergers and leads to a high merger rate in AGN disks even without pre-existing binaries. We find the merger rate to be in the range $\sim 0.02-60\,\mathrm{Gpc^{-3}yr^{-1}}$, whose high end corresponds to top heavy initial mass functions, highly anisotropic initial BH distribution, short AGN lifetime, and large AGN disk sizes. The results are insensitive to the assumptions on gaseous hardening processes: we find that once they are formed, binaries merge efficiently via binary-single interactions even if these gaseous processes are neglected. We find that the average number of mergers per BH is $0.4$, and the probability for repeated mergers in 30 Myr is $\sim 0.21-0.45$. High BH masses due to repeated mergers, high eccentricities, and a significant Doppler drift of GWs are promising signatures which distinguish this merger channel from others. Furthermore, we find that gas-capture binaries reproduce the distribution of LMXBs in the Galactic center, including an outer cutoff at $\sim1$ pc due to the competition between migration and hardening by gas torques.

Twelve low-mass X-ray binary (LMXB) candidates have also recently been discovered by Hailey et al. (2018) within a distance r 1 pc of the Galactic center, with a density distribution of ∝ r −1.5±0.3 , by observing Xray sources within ∼ 4 pc presented by Muno et al. (2009). Generozov et al. (2018) proposed that the density profile of these hard binaries can be explained by the tidal capture mechanism and stellar relaxation processes. Although this model predicts the radial distribution of LMXBs to be ∝ r −(0.9−1.4) , the outer cut-off at ∼ 1 pc remains unexplained.
Galactic nuclei are the densest environments of stars and compact objects in the Universe (Walcher et al. 2005;Merritt 2010;Norris et al. 2014;Gallego-Cano et al. 2018). In the active phase of a galactic nucleus, a highdensity gas disk forms within 0.1-10 pc (Burtscher et al. 2013) of a central supermassive BH (SMBH). In such environments, binaries form and evolve through interaction with densely populated stars and gas. Baruteau et al. (2011) showed that even when a binary is so mas- sive that it opens a gap within the accretion disk around a SMBH, it is efficiently hardened via gas dynamical friction. Bartos et al. (2017) have proposed a pathway for BH-BH mergers in active galactic nucleus (AGN) disks, in which binaries are captured in an accretion disk within ∼ 0.01 pc of the central SMBH due to linear momentum exchange, and after that, binaries are hardened by gas dynamical friction from the AGN disk and by type I/II torques from a circumbinary disk. Stone et al. (2017) proposed another pathway, in which in-situ formed binaries at ∼pc scale evolve via the effects of binary-single interactions with a stellar disk and via type I/II torques from a circumbinary disk. Leigh et al. (2018) showed that fewer than ten binary-single interactions are sufficient to drive hard binaries with a binary separation of s 10 AU to merger. McKernan et al. (2018McKernan et al. ( , 2019 estimated the mass and spin distributions of the merged BHs, and the merger rate in AGN disks. Bellovary et al. (2016) suggested that BHs accumulate and merge with each other in migration traps at 20 − 300 Schwarzschild radii, where the torque from the AGN disk changes sign. Secunda et al. (2019) and Yang et al. (2019a) modeled the formation of binaries within the migration traps. Just et al. (2012), Kennedy et al. (2016) and Panamarev et al. (2018) discussed the capture of stars in an AGN disk and their subsequent migration toward the central SMBH due to the ram pressure from an AGN disk.
Previous studies of compact object mergers in AGN disks have focused on the role of gas in driving binary mergers assuming pre-existing binaries in the nuclear star cluster (Bartos et al. 2017) or in the disk itself (Stone et al. 2017;McKernan et al. 2018), or assuming that binaries form at migration traps (Bellovary et al. 2016;Secunda et al. 2019;Yang et al. 2019a,b;Gayathri et al. 2019). In the present study, we examine the formation of binaries during close two-body encounters in a gaseous medium, where the gas absorbs some of the initial kinetic energy of the two objects. Goldreich et al. (2002) have proposed that planetesimal binaries can form in the Kuiper belt due to the dissipation of the relative velocity between the planetesimals by dynamical friction in the environment of a background population of smaller solid bodies. They showed that this leads to efficient binary formation in the Kuiper belt, but to our knowledge, the analogous mechanism of binary formation in an AGN disk has not been previously explored. Here we include this "gas capture" binary formation mechanism, and find that it supplies the majority of binaries in AGN disks. We also examine binary formation in dynamical three-body encounters (e.g. Aarseth & Heggie 1976;Binney & Tremaine 2008), which have also been neglected in previous studies of mergers in AGN disks, but find this mechanism to be less important.
More generally, in this paper we investigate whether and how binaries form and merge in AGN disks. We combine N -body simulations with a semi-analytical model, which incorporates the effects of gas dynamical friction, type I/II migration torques, GW radiation, and several different types of stellar interaction. We simulate the evolution of both single and binary objects, and follow their radial position from the central SMBH, as well as their velocities in time, together with the evolution of the binaries' separation. Our flexible model allows us to test previous assumptions on whether and how efficiently binaries may merge in AGN disks (e.g. McKernan et al. 2019, and references above), and to examine the dependence of the merger rate on the model parameters of the AGN disk and the stellar cluster.
The rest of this paper is organized as follows. In § 2, we give an executive summary of our method, and describe the numerical scheme and the setup of simulations in detail in § 3. We present our main results in § 4. We discuss the implications for the spin and eccentricity distribution, the merger rate, and a comparison with LMXBs observed in the Galactic nucleus in § 5. We summarize our conclusions in § 6. For clarity, the variables used in this paper are listed and defined in Table 3 in the Appendix.
2. OVERVIEW OF SIMULATIONS We consider a system describing a galactic nucleus, consisting of the following five components: (1) a central SMBH, (2) a gaseous accretion disk ("AGN disk"), (3) a spherical stellar cluster, (4) a flattened cluster of BHs, and (5) stars and BHs inside the AGN disk, referred to as the "disk stellar" and "disk BH" components ( §3.1.2). Figure 1 illustrates our setup. In our fiducial model, we adopt the SMBH mass and the distribution of stars from observations of the quiescent central region of the Milky Way at present ( § 3.1.1), which does not have an AGN. Also we generate the BH mass distribution using the results of population synthesis models and accounting for an initially mass segregated radial distribution ( § 3.1.1). On the other hand, our model with an AGN disk represents the conditions during the active phase which is believed to have existed at an earlier time in the Milky Way's history (e.g. Zubovas et al. 2011).
We employ the AGN accretion disk model proposed by Thompson et al. (2005) ( § 3.1.3), as adopted in earlier N-body BHs/binaries Radial distance, velocity, and binary separation evolve via semi-analytical prescription Gaseous interaction (type I/II migration torque, gas dynamical friction, accretion torque) -in the AGN disk -in a circumbinary disk Interaction with single stars/BHs and BH binaries due to (binary-single encounter, weak scattering) -a spherical stellar component -a disk BH component -a disk stellar component GW emission (GW torque, post-merger recoil) Binary formation and disruption (three-body encounter, two-body gascapture, binary-single interaction) Update the number density, velocity dispersion, average mass, and typical orbital height for each component in each radial cell Figure 2. Schematic diagram for following the evolution of the BH population in our model, including both single and binary BHs. The N-body simulation keeps track of these individual objects (starting with 2×10 4 BHs and 1.5×10 3 binaries in our fiducial model). Single BHs are characterized by their radial position (r i ) from the central SMBH and by their velocity (v i ) relative to the local Keplerian value. Binaries are similarly characterized by their center-of-mass position (r j ) and velocity (v j ), and additionally by their orbital separation (s j ). These variables are updated via semianalytic prescriptions in each "N-body" time-step, due to multiple processes as listed in the diagram.
work by Stone et al. (2017). This represents a Shakura-Sunyaev α-disk with a constant viscosity parameter α and accretion rate in the region where it is not selfgravitating. The model describes a radiatively efficient, geometrically thin, and optically thick disk and extends the disk to pc scales with a constant Toomre parameter in the self-gravitating regime (see Fig. 4 below), assuming that it is stabilized by radiation pressure and supernovae from in-situ star formation. We assume that stars and BHs form in the disk at the rate required to stabilize the AGN disk, and some fraction of BHs are initially formed in binaries ( § 3.1.3).
To follow the time-evolution of the BHs in this system, focusing on their capture by the disk, and the formation and disruption of BH binaries in the disk, we run N -body simulations combined with a semi-analytical method. Binaries form in the disk either due to gas dynamical fric- tion ( § 3.3.8) or due to three-body encounters ( § 3.3.7), and are disrupted by soft binary-single interactions. We assume that binaries are disrupted when the binary separation becomes larger than the Hill radius of a binary with respect to the SMBH.
We model the evolution of the orbital separation (s j ), the radial position (r j ), and the magnitude of a random velocity relative to the local Keplerian AGN disk motion (v j ) for all binaries labelled with the binary index j. We also track the radial position (r i ) and random velocity (v i ) of single compact objects, which represent stellarmass BHs. These quantities evolve according to analytical formulas as summarized in Figure 2 and illustrated visually in Figure 3.
Our model includes physical processes due both to the presence of gas and to multi-body dynamical interactions, as follows. For the interaction with gas, the radial positions of all BHs evolve in response to (i) type I/II migration torques by the AGN disk ( §3.3.2), and the velocities of all BHs relative to the local AGN disk decrease due to (ii) the accretion torque ( §3.3.4), and (iii) gas dynamical friction in the AGN disk ( §3.3.3). For binaries of stellar-mass BHs, the separation evolves due to gas dynamical friction from the AGN disk and to type I/II migration torque from a small circumbinary disk that forms within the Hill sphere of the binary. We also account for dynamical interactions with single stars and BHs and BH binaries: the binaries' separations and velocities evolve due to binary-single interactions ( §3.3.5), and the velocities of all BHs additionally evolve due to scattering ( §3.3.6). For dynamical interactions, we neglect the BHs in the spherical cluster, since they are greatly outnumbered by stars in the cluster, but we include both stars and BHs in the disk component. We also account for GW emission, which reduces the binary separation rapidly once the binary is sufficiently tight. For simplicity, eccentricity evolution is ignored and orbits around the SMBH and binary orbits are both assumed to be circular.
To model the orbits of merged BHs, a recoil velocity is added to the BH remnant due to anisotropic GW radia-tion ( §3.3.10). The small mass loss during mergers due to GW radiation is taken into account assuming zero BH spins.
In this study, we ignored several processes for simplicity. These include exchange of binary components at binary-single interaction, the formation and evolution of stellar binaries, radial migration of stars due to torque from the AGN disk, evolution of compact objects other than BHs, the Kozai-Lidov effect of the SMBH or a third stellar-mass object on binaries, dynamical relaxation processes, counter rotating BHs or stars in the AGN disk (Ivanov et al. 2015;Sánchez-Salcedo et al. 2018), stellar evolution, supernova feedback, binary mass transfer, and possible presence of massive perturbers like an SMBH companion and/or IMBHs. A few IMBHs, if present, may efficiently disrupt most BH binaries which may greatly reduce the merger rates (Deme et al. 2019).
The above model allows us to describe the timeevolution of the binary BH population in a self-consistent and flexible way. It extends the simplified prescriptions of previous studies of stellar-mass BH binary mergers in AGN disks, and creates a self-consistent semi-analytical N -body simulation that includes the time-dependent formation, disruption, and evolution of binaries in AGN. We use this method to estimate the contribution of binaries formed during AGN phase. We confirm previous suggestions (Secunda et al. 2019;Yang et al. 2019b) that repeated mergers are frequent in AGN disks, although in difference from these previous works, in our models repeated mergers occur due to efficient binary formation and evolution processes well outside the "migration traps".
3. METHOD Here we describe in detail the method and the initial conditions adopted in this study.

Stellar mass BHs, stars and AGN disk
In this section we describe the initial condition in the calculations.
3.1.1. Initial BH and stellar distribution We simulate the evolution of N -body particles representing stellar-mass BHs. We assume these are initially distributed according to where N BH,ini (r) labels the total initial number of BHs within distance r from the central SMBH, and γ ρ is a power-law index. Theoretically, γ ρ is expected to be between ∼ 0.5 and 1.25 for plausible mass functions for spherically symmetric systems (Hopman & Alexander 2006;Freitag et al. 2006;O'Leary et al. 2009;Keshet et al. 2009;). In our fiducial model, we adopt γ ρ = 1 between r in,BH ≤ r ≤ r out,BH where r in,BH = 10 −4 pc and r out,BH = 3 pc. We set the total stellar mass within 3 pc to be M star,3pc = 10 7 M (2)  as the fiducial value. The minimum and maximum mass for progenitor stars are assumed to be 0.1 and 140 M , respectively. The BH mass is determined through the relations between the progenitor mass (m star,i ) and the BH mass (m BH,i ) of which roughly matches population synthesis simulation results in Belczynski et al. (2010) for their model with solar metallicity and a weak wind. Since observational studies (Bartko et al. 2010;Lu et al. 2013) suggest a top-heavy initial mass function (IMF) for stars in the Galactic center region, we investigate IMFs with δ IMF in the range 1.7 − 2.35. We set the fiducial value to be δ IMF = 2.35, yielding an average stellar mass m star = 0.36 M and initial number of BHs N ini,BH = 2.0 × 10 4 . For 1.7 ≤ δ IMF ≤ 2.35,m star and N ini,BH vary between 0.36 − 1.78 M and 2.0 × 10 4 − 1.0 × 10 5 . The simulation tracks the velocity of particles relative to the local Keplerian AGN disk in the plane of the disk v xy,k and perpendicular to it v z,k at the point where the orbit crosses the equatorial plane, where k is the particle index. The direction of v xy,k is assumed to be axisymmetrically random in the xy plane. The x−, y−, and z− components of the velocity of each BH relative to the local disk are initially drawn randomly from a Gaussian distribution with dispersion of β v v Kep (r)/ √ 3 and zero mean. Here v Kep (r) = {G[M SMBH + M star (< r)]/r} 1/2 is the Keplerian orbital velocity at the distance r from the central SMBH, M star (< r) is the stellar mass within r, and β v is a parameter that determines the initial velocity dispersion of BHs. In our fiducial model we set β v = 0.2. Since this is with respect to the comoving Keplerian frame, it corresponds to a net rotation for the pre-existing BHs component, which is consistent with observational (e.g. Trippe et al. 2008;Yelda et al. 2014;Feldmeier et al. 2014Feldmeier et al. , 2015 and theoretical (Kocsis & Tremaine 2011;Szolgyen & Kocsis 2018) suggestions that massive stars in the central region of 1 pc have some degree of net rotation.
We assume that some fraction of BHs are initially in binaries as follows. Spectroscopic observations show that the binary fraction of O stars in the Galactic field is ∼ 0.7 (Sana et al. 2012), but the binary fraction of OB/WR stars in the Galactic center is estimated to be only ∼ 0.3 (Pfuhl et al. 2014). We define the corresponding initial "pre-existing" BH binary fraction f pre as the number of BH binaries over the total number of (single+binary) BHs. Belczynski et al. (2004) found that if the binary fraction of progenitor stars is 50% (2/3 of stars are in binaries), the binary fraction of BHs is ∼ 10 % as a result of stellar evolution due to supernova kicks and mergers during the common envelope and Roche-lobe overflow phases. In the Galactic center, binary disruption due to soft binary-single interactions may further decrease the BH binary fraction (see also Stephan et al. 2016). We adopt f pre = 0.15 in our fiducial model.
We draw the initial separation of pre-existing binaries randomly from a log-flat distribution following Abt (1983). The minimum separation R min has large uncertainties since R min is determined by common envelope evolution, which is not well understood (Ivanova et al. 2013). For the fiducial value, we set R min by the sum of the radii of the progenitor binary components. We compute the stellar radius as ) We set the maximum binary separation R max to 10 5 R following binary evolution models (e.g. Belczynski et al. 2008;Kinugawa et al. 2014). We further assume that binaries that are soft compared to the local spherical stellar component ( §3.3.5) are promptly disrupted prior to the AGN phase. The timescale of the disruption of binaries is ∼ 300 Myr at r j 0.01 pc and ∼ 20 Myr at r j = 3 pc (Eq. 7.173 of Binney & Tremaine 2008), in which we use the velocity dispersion, the density, and the average stellar mass for the spherical stellar component in the fiducial setting of our simulations (Eq. 6), the Coulomb logarithm is assumed to be 10, each BH binary component is assumed to have a mass of 5 M , and the binary separation is the maximum of the hard-soft boundary (Eq. 85 below) and the separation at which a binary merges within the Hubble time (Eq. 86 below). Due to the disruption of binaries prior to the AGN phase, the binary fraction at the beginning of the simulation is reduced to ∼ 7%. Note that this value varies according to the initial mass function, the total stellar mass, and the mass of the SMBH (see § 4.4 below).

Stellar and BH components
We categorize stars and BHs by whether they reside within or orbit outside of the AGN disk. These components are referred to as the disk stellar, the spherical stellar, the disk BH, and the anisotropic BH components ( Figure 1). We assume that stars are initially spherically distributed, and the velocity of stars follow a Maxwell-Boltzman distribution with no net rotation, while BHs are initially distributed with some degree of net rotation (see justification below). Due to the interaction with the AGN disk and star formation in the outer regions of the AGN disk, the number of BHs and stars in the AGN disk gradually increase. As a result, the BH and stellar density within the AGN disk rise. Thus the BH and stellar disk components form during the AGN phase ( Figure 1). Since the mass outside of the AGN disk is dominated by stars (i.e. the spherical stellar component), we ignored the interaction of BHs/binaries with BHs orbiting outside of the AGN disk (i.e the anisotropic BH component) for simplicity.
To compute the rates of various density and velocitydependent processes (see § 3.3.5-3.3.8 below), we calculate the number density (n c ), the velocity dispersion (σ c ), the average mass (m c ), and the typical orbital height (h c ) for each component in each radial cell in a grid. The spherical radial grid extends from r in = 10 −4 to r out =5 pc, and is divided into N cell = 120 cells uniformly on a log scale in the fiducial model. The dependence of results on N cell is discussed in §4.4. We neglect the pos-sible effects related to migration traps, which may exist at 9 × 10 −5 (M SMBH /4 × 10 6 M ) pc for the model by Thompson et al. (2005) (Bellovary et al. 2016), just outside of the simulated domain.
We set the time-independent physical quantities for the spherical stellar component in each cell to be wherem star is the average stellar mass in the spherical stellar component determined from the assumed IMF Eq. (4), v Kep,l is the Keplerian velocity at the radius r l of the geometric center of a cell l, and ρ Sstar (r) is the density profile of the spherical stellar component. We adopt chosen to match the observed stellar surface density distribution in the Galactic nucleus (Merritt 2010;Feldmeier et al. 2014), where M star,3pc is given by Eq. (2). The typical orbital height for the spherical stellar component (h c ) Sstar,l is set to r l / √ 2 considering a uniform distribution in each spherical shell.
On the other hand, the quantities for the disk BH component in each cell are initialized as which evolve with time, where DBH,l refers to the BH components which are within the AGN disk in the lth cell, N DBH,l is the number of single BHs and BH binaries within the AGN disk in cell l, V l = 4πr 2 l ∆r l is the spatial volume of the lth spherical shell so that h l V l /r l is the spatial volume of the AGN disk in the lth cell, and ∆r l is the radial width of the lth cell. Note that we include the number of BH binaries in N DBH,l . For the binaries, m k and v k refer to the total mass and center of mass velocity of the k th binary relative to the Keplerian velocity in Eq. (8). The average mass m c,DBH,l and the velocity dispersion σ c,DBH,l of the disk BH component in a cell l are given by the the average of m k and the root mean square of v k / √ 3 of BHs, respectively. We assume that when h k < h AGN,l , the k th object is embedded in the AGN disk, where h k = v z,k r k /v Kep is the typical height of orbital motion for the k th object and h AGN,l is the height of the AGN disk at r l (h AGN,l is derived in § 3.1.3). We assume that the disk BH component rotates in the same sense as the AGN disk, so σ c = 0 means that BHs in the disk co-rotate with the Keplerian gas.
The statistical quantities of the disk stellar component in each cell are calculated as (m c , n c , σ c , h c ) Dstar,l =(m star , n Dstar,l , σ c,DBH,l h c,DBH,l ) where Dstar, l refers to the stellar component embedded within the AGN disk in the lth cell. We ignore the accretion and migration of stars captured by the AGN disk for simplicity. This is a conservative assumption for the merger fraction since it neglects the possibility for migration to increase the stellar density in the inner regions which would facilitate BH mergers by frequent binary-single interactions. Since we do not calculate the evolution of stars, we assume that their velocity dispersion and scale height match that of the BH disk component, given by Eq. (8). The density of the disk stellar component is calculated from the number of the stars by simply assuming that the stars reside in the same volume as the disk BH component (V l h c,DBH,l /r l ). The number of stars in the disk stellar component is calculated considering three factors. First we assume that stars form with rateΣ * in the outer regions of the AGN disk ( Figure 4, § 3.1.3). We ignore the evolution of newly formed stars. Second we assume that spherically distributed stars are captured in the AGN disk at the rate estimated in §3 of Bartos et al. (2017). Bartos et al. (2017) estimated the timescale on which objects are captured in the AGN disk based on the torque due to Bondi-Hoyle-Lyttleton accretion during crossing the AGN disk. In our simulation, we calculate the critical inclination angle of stellar orbits with respect to the AGN disk at which the alignment timescale (Eq. 11 of Bartos et al. 2017) becomes the same as the elapsed time. We assume that the inclination of stars cos i is distributed uniformly between -1 and 1 as in Bartos et al. (2017), and derive the fraction of stars whose inclination is smaller than the critical inclination angle. In this way we calculate the number of stars captured in the AGN disk in each step. Third, we reduce the number of stars in the disk stellar component in the lth cell hosting the binary by one for each BH binary that experiences a hard binary-stellar single interaction with objects in the disk stellar component. This reflects the fact that the recoil kick following a hard BH binary-single interaction with the typical low-mass stars is so large that the interacting stars are usually kicked out from the AGN disk in the vertical direction.
3.1.3. AGN disk We employ the AGN accretion disk model proposed by Thompson et al. (2005). In the fiducial model, we set the SMBH mass to M SMBH = 4 × 10 6 M , and the accretion rate from outer boundary of r out = 5 pc to beṀ out = 0.1Ṁ Edd , whereṀ Edd = L Edd /(c 2 η c ) is the Eddington accretion rate, L Edd is the Eddington luminosity, c is the light speed, and η c is the radiative efficiency. We adopt η c = 0.1 assuming a standard thin disk model. We adopt the opacity model given by Bell & Lin (1994) which gives the opacity as a function of temperature and density. Following the fiducial values in Thompson et al. (2005), we assume the pressure ratio parameter ξ = 1. We improve the calculation of the   Figure 4. Physical quantities for the adopted disk model as a function of the distance r from the SMBH. The colored lines represent the disk height over the distance h disk /r (black), the temperature (red), the gas density ρ AGN (blue), the background stellar density (orange), the star formation surface densityΣ * (green), and the accretion rateṀ d (brown). Units are M , Myr, • K and pc.
conversion efficiency of star formation to radiation in Thompson et al. (2005) by taking into account the limitations due to the AGN and stellar lifetimes (see Appendix § B for details). This reduces by a factor ∼ 4 to = 1.5 × 10 −4 for the fiducial model, and varies between 1.5 × 10 −4 and 7.7 × 10 −4 according to the IMF exponent (−2.35 ≤ −δ IMF ≤ −1.7). We assume that the efficiency of angular momentum transport due to global torque in the outer region is m AM ∼ 0.1 − 0.2 as suggested by Thompson et al. (2005). In the inner region, we adopt the α-model for angular momentum transport, in which the alpha parameter is α SS = 0.1 (King et al. 2007;Bai & Stone 2013). In the transition between the inner and outer regions, the degree of angular momentum transport is adjusted to keep the Toomre parameter at Q = 1. We assume a locally isothermal equation of state to calculate the sound speed of the AGN disk as c s = (p/ρ gas ) 1/2 , where p is the total gas+radiation pressure, and ρ gas is the local density of a gas disk. The viscosity is given by ν = rv r /(dlnΩ/dlnr), where v r is the radial velocity, and Ω is the angular velocity of the gas disk, which is rotating around the enclosed mass of the SMBH and stars (Eq. 7). In the α-prescription, the viscosity ν is assumed to be proportional to the total pressure, ν = α SS c 2 s /Ω. The surface density of the gas disk is calculated using the radial velocity and the gas inflow rate (Ṁ ) as Σ disk =Ṁ /(2πrv r ).
The AGN disk properties are shown in Figure 4. The outer region ( 1 pc) is stabilized by radiation pressure and supernovae from in-situ formed stars. Therefore the star formation rate is determined by the disk model, which depends on the mass of the SMBH, the accretion rate, and efficiency parameters through Eqs. (C1)-(C12) of Thompson et al. (2005). The green line in Figure 4 represents the surface density of the star formation ratė Σ * . The star formation rate and BH formation rate surface densities are given by f star ×Σ * and f BH ×Σ * , respectively, where f star and f BH are the ratio of the mass of stars with the mass less than 20 M and all stellarmass BHs to the mass of all stars at formation, respectively. BHs continuously form at this rate in our simulations, and the number of BHs formed within 100 Myr is ∼ 7% of the initial number of BHs in the fiducial model. The accretion disk model with large size and star formation is motivated by several observations, which are early enhancement of metallicity in disks (Artymowicz et al. 1993;Xu et al. 2018;Novak et al. 2019), long-timescale transients in AGN (Graham et al. 2017), and supernovae found in a vicinity of an AGN (Pérez- Torres et al. 2010).
Based on Eq.
(3), f BH varies from 0.016 to 0.092 according to the initial mass function (1.7 ≤ δ IMF ≤ 2.35), and f BH = 0.016 in the fiducial model. We set the velocity of newly formed BHs relative to the local AGN motion v k to be the sonic velocity of the AGN disk at their location. This is motivated by numerical simulations suggesting that the scale height of newly born stars within an AGN disk is roughly similar to the thickness of the AGN disk (Nayakshin et al. 2007). In this study, BHs form immediately during star formation neglecting the lifetime of the progenitor stars for simplicity. We assume that the mass distribution and the binary fraction of in-situ formed BHs are the same at formation as those of the pre-existing BHs described above.
3.2. Formation, destruction, and orbital evolution of binaries In this study we do not follow the evolution of individual stars, but only track their average statistical properties in a grid as explained in § 3.1.2 for both the spherical and the disk stellar components. However we follow the orbital parameters of individual BHs and BH binaries as follows.
For clarity we use the indices i and j in this paper exclusively to label a single BH and a BH binary, respectively. We use the index k to denote either a single BH or a BH binary. For single BHs, we characterise their orbits with their orbital radius around the SMBH r i , their magnitude of z-direction velocity v z,i , and the xy-direction velocity v xy,i at z = 0. Here v i = (v 2 xy,i + v 2 z,i ) 1/2 is the velocity relative to the local Keplerian velocity of a circular orbit with z = 0. Given v z,i , the maximum height of the orbit is h i = r i v z,i /v Kep (r i ). For binaries, we follow their radial position r j , their center of mass velocity components v xy,j and v z,j , and their binary separation s j . The position, velocity, and separation of BHs evolve via interacting with gas and/or stellar and BH components ( Figure 3). In this paper, we do not consider the evolution of the binary eccentricity (e j ), and assume e j = 0 for all binaries.
We incorporate the effects of gas dynamical friction, GW radiation, weak gravitational scattering, binarysingle interaction, type I/II migration torque, binary formation via three-body encounter and via gas-capture mechanism, binary disruption, and star formation. The velocity v i = (v xy,i , v z,i ) of a single BH evolves via the equation of motion (Papaloizou & Larwood 2000) dv where a acc,i is the acceleration due to the accretion torque ( §3.3.4), gas dynamical friction a GDF,i ( §3.3.3), and weak gravitational scattering a WS,i ( §3.3.6). For BH binaries, the center of mass velocity v j changes due to all of these processes and additionally due to binary-single interactions as dv j dt = a acc,j + a GDF,j + a WS,j + a BS,j , where a BS,j is the acceleration due to binary-single interactions (a BS,j . Bartos et al. (2017) proposed that the binary separation decreases due to gas dynamical friction from an AGN disk, type I/II torque from a circumbinary disk, and GW radiation after the capture of close-in binaries in the smooth disk in ≤ 0.01 pc. Stone et al. (2017) considered mergers from binaries formed in-situ at the unstable part of the disk at larger radii of ∼ pc. Stone et al. (2017) find that binary-single interaction with the disk stellar component also harden the binary separation. In this study, we incorporate these effects in the evolution of the binary separation as where the four terms on the right hand side are the evolution rates for the binary separation due to GW radiation ( §3.3.1), gaseous torque ( §3.3.2 and §3.3.3), and binarysingle interaction, respectively. Single and binary BHs migrate radially toward the SMBH due to the gaseous torque of the AGN disk, given by type I or type II torque formulae ( §3.3.2).
Furthermore, binaries form and are disrupted according to dN bin dt = where P 3bbf,i is the binary formation rate by threebody encounter ( §3.3.7), P gas,i is the binary formation rate by gas-capture mechanism ( §3.3.8), and K dis = dN bin /dt| dis < 0 is the binary disruption rate. Binaries are disrupted in the simulation when the binary separation s j becomes larger than the Hill radius of a binary with respect to the SMBH The terms in Eqs. (10)-(14) are described in the next section, §3.3. These equations are calculated separately using the local statistical quantities describing the stellar environment in each shell l.

Individual Processes
Here we describe in detail the prescription adopted for each of the mechanisms included in our simulations.

Gravitational wave radiation
When the binary separation is small, GW radiation strongly decreases the binary separation. The hardening rate via GW radiation is given as assuming zero eccentricity (Peters 1964), where m j1 and m j2 are the masses of the primary and secondary BH in the j th binary, respectively.
3.3.2. Type I and II migration Objects in a gaseous disk interact gravitationally with nearby gas, resulting in radial migration. When the gravitational torque exerted by an object within the disk exceeds the viscous torque of gas, a gap opens in the disk around the object (e.g. Ward 1997;Crida et al. 2006). When a gap does not open around the object, the object migrates due to torques from the Lindblad and corotation resonances on the type I migration timescale (e.g. Ward 1997;Tanaka et al. 2002;Paardekooper et al. 2010;Baruteau et al. 2011) of where f mig is a dimensionless factor depending on the local temperature and density profiles (see Paardekooper et al. 2010; Baruteau et al. 2011), h disk is the half thickness of the disk, M cen and M sat are the central and satellite object mass, respectively, and Σ disk is the surface density of the disk. 1 We calculate the evolution in the range r ≥ r in ∼ 10 3 R g ) where f mig remains positive and its variation is not significant. We set f mig = 2 in the fiducial model, which is the typical value found numerically by Kanagawa et al. (2018). If a gap opens in a disk around an object, the object migrates due to the torque of the gas approaching the gap boundary on a timescale related to the viscosity. This process is the so-called type II migration (e.g. Lin & Papaloizou 1986;Ward 1997;Ida & Lin 2004;Edgar 2007;Haiman et al. 2009;Duffell et al. 2014). Recent hydrodynamic simulations Durmann & Kley 2015;Kanagawa et al. 2018) have shown that even when the torque from a BH exceeds the viscous torque of gas, the gas is able to pass through the gap. Duffell et al. (2014) and Kanagawa et al. (2018) show that the migration timescale even for massive migrator is given by the type I migration timescale with a reduced gas surface density in the gap (Σ disk,min ), Fung et al. (2014) and Kanagawa et al. (2015) show that where We calculate the migration rate of stellar-mass BHs within the AGN disk as dr where p disk,k is the fraction of time that the k th object spends in the disk along its orbit around the SMBH, which we calculate as (21) assuming that the time spent inside the AGN disk is approximated by the ratio of the scale height of the disk and the BH's orbit. To calculate the migration rate, we substitute h AGN,l into h disk in Eq. (17). To reduce the computational cost, we use quantities at the center of each cell r l in Eqs. (17)-(18). For comparison, we also investigate cases in which migration by the type I/II torque does not operate, since the migration of BHs is not well understood due to the complexity of the effects of Nbody migrators (Broz et al. 2018), feedback from BHs Due to the torques exerted by the BHs, the gas density is reduced near each BH according to Eq. (19). We take into account this reduction when we use gas dynamical friction ( §3.3.3), gas accretion ( §3.3.4), and the gas-capture mechanism ( §3.3.8). In these mechanisms, when a BH is in the AGN disk, we use the local gas density of the AGN disk around each BH as where ρ AGN,l is the unperturbed density of the AGN disk in the lth cell hosting a BH ( §3.1.3). When a BH orbits outside of the AGN disk (h z,k > h AGN,l ), we use ρ gas = ρ AGN,l = Σ disk,l /(2h disk,l ). After gas is captured within the Hill's sphere of a binary, a circumbinary disk forms. This disk exerts a torque on the binary and changes its separation similar to Type I/II migration, although a consensus is beginning to emerge in recent hydrodynamical simulations that an equal-mass binary is softened, rather than hardened by the presence of a circumbinary gas disk (Miranda et al. 2017;Tang et al. 2017;Muñoz et al. 2019;Moody et al. 2019;Duffell et al. 2019). We assume that the binary can be hardened by the torque of a circumbinary disk, and the hardening timescale is given as Eq. (20) by substituting s j into r k , Here, we substitute the BH binary component masses for M cen and M sat , and use the angular velocity of the binary in Ω in Eq. (20). If the Toomre parameter of the circumbinary disk satisfies Q > 1, it is stable against gravitational fragmentation. In this case, we calculate the surface density of a circumbinary disk using Eq. (14) in Goodman & Tan (2004), in which we assign the opacity consistently following Bell & Lin (1994), and the disk temperature is given by the maximum of Eq. (13) in Goodman & Tan (2004) and the temperature of the AGN disk at the position r k . We assume that the gas pressure dominates over the radiation pressure for the circumbinary disk, which is a valid approximation for the radial range of the disk ( ∼ 10 9 cm) we are interested in. When Q < 1, we reduce the surface density of gas disk around a binary to satisfy Q = 1. We assume the accretion rate onto a binary or a single BH is given by the minimum of the Bondi-Hoyle-Lyttleton rate (Eq. 29) times p disk,k (Eq. 21) and Eddington limited accretion, where Γ Edd,cir is the Eddington ratio, L Edd,k is the Eddington luminosity for a binary, and η c is the energy conversion efficiency. We set Γ Edd,cir = 1 and η c = 0.1. The Eddington limited accretion rate is motivated by regulation due to strong radiation pressure acting on dust grains for gas of around solar metallicity (Toyouchi et al. 2019), by a bipolar jet (Regan et al. 2019), and also by inefficient angular momentum transfer in circumbinary disks (Sugimura et al. 2018;Inayoshi et al. 2018). We assume that the type I/II migration torque generated by a circumbinary disk surrounding the binary operates to shrink the binary separation only when the Bondi-Hoyle-Lyttleton radius r BHL,k is larger than the binary separation and otherwise the circumbinary disk does not exert a torque on the binary. Furthermore, we assume that a circumbinary disk is always aligned with the binary. This is justified since the alignment timescale of the disk with the binary is roughly given by the viscous timescale at the binary separation (e.g. Ivanov et al. 1999;Moody et al. 2019). This timescale is r 2 /ν ∼ α −1 SS (h disk /r) −2 Ω −1 ∼ 10 5 yr (Moody et al. 2019) for s j ∼ AU, which is shorter than the evolution timescale of binaries by type I/II torque from circumbinary disks ( Myr).

Gas dynamical friction
When an object has a non-zero velocity relative to the ambient gas, gas dynamical friction reduces the relative velocity. In this manner, gas dynamical friction hardens binaries (Escala et al. 2004;Kim & Kim 2007;Baruteau et al. 2011) and damps the velocity dispersion of BHs and stars (Papaloizou & Larwood 2000;Tanaka et al. 2002) in the disk. For simplicity, we adopted the formulation for deceleration by gas dynamical friction derived by Ostriker (1999) as where ln Λ gas is the Coulomb logarithm for gas. We set ln Λ gas = 3.1 referring the results in Chapon et al. (2013). For both binary hardening and velocity damping, we use the sonic velocity c s from the disk model of Thompson et al. (2005) ( §3.1.3) at the geometric center of the cell l hosting BHs. Although these formulae apply for linear motion, they remain approximately correct for circular motion, despite the strong curvature and possible interaction of the pair of "wakes" in the binary case (Kim & Kim 2007). For binary hardening, we assume that gas dynamical friction operates while the binary is captured to the AGN disk (h k < h AGN ). The AGN disk capture timescale on which the initial supersonic velocity of an object v ini,k decays due to crossing the disk due to gas dynamical friction is where we introduced the abbreviated labels Figure 4), where we assumed f ∼ 1 in a GDF for simplicity (see Eq. 25), and ignore the contribution of the stellar mass as v Kep ∼ (GM SMBH /r k ) 1/2 , which underestimates t capAGN by ∼ 2 at r k = pc. Similarly we define the gas dynamical friction hardening timescale of a binary as where v HS = (m star /m k ) 1/2 v Kep (r) is the orbital velocity of a binary of mass m k around its center of mass at which the binary is at the hard-soft boundary compared to the spherical stellar component. During the hardening of binaries, it is not obvious whether type I/II torques by a circumbinary disk or gas dynamical friction by the AGN disk is a better description. Baruteau et al. (2011) showed that a binary in an AGN disk is hardened roughly on the timescale of gas dynamical friction. Derdzinski et al. (2019) find that the gas captured by a rapidly migrating BH within an AGN disk has almost no rotation with respect to the migrator, which suggests that gas torques may operate in the manner of the dynamical friction (although this was demonstrated only for a single specific binary + disk configuration). Since this issue has not been settled yet, we investigate several different prescriptions. In the fiducial model, we investigate the cases in which the binary hardening rate due to gas interaction is given by (i) the maximum value of the gas dynamical friction and the type I/II torque, (ii) only the gas dynamical friction, (iii) only the type I/II torque, or (iv) zero (no hardening by gas interaction). Park & Bogdanovic (2017) have shown that gas dynamical friction does not operate when the Bondi-Hoyle-Lyttleton radius (r BHL,k = Gm k /(c 2 s + v 2 k )) is smaller than the size of a HII sphere (r HII,k = (3Q ion,k /4πα rec,B n gas ) 1/3 ), where α rec,B is the case-B recombination coefficient for H (evaluated at T = 10 4 K), and Q ion,k is the ionizing photon number flux from the k th BH. This is caused by radiation feedback from a BH which diminishes the wake created by gas dynamical friction. The condition r BHL,k > r HII,k can be rewritten as (28) if we assume the Eddington accretion rate with the radiative efficiency η c = 0.1, c s v k , and Q ion ∼ L k /hν p with hν p ∼ 13.6 eV, where h is Planck's constant, ν p is the average photon's frequency, and L k is the luminosity from the k th BH. Also when the velocity of a BH exceeds twice the sound velocity in the HII region (v k 50 km/s), gas dynamical friction is recovered (Park & Ricotti 2013;Park & Bogdanovic 2017). We set these conditions, r BHL,k > r HII,k or v k > 50 km/s, as criteria for gas dynamical friction to operate.

Accretion torque
As a BH crosses the AGN disk, it captures gas from the disk. The velocity of the BH decreases to satisfy conservation of momentum. The capture rate of gas on the k th BH during passing the AGN disk is given by the Bondi-Hoyle-Lyttleton ratė where r w,k = min(r BHL,k , r Hill,k , r shear,k ) and are the width and height of the gas bound to the k th BH, and is the typical radius over which gas motion changes due to the Keplerian shear of the AGN disk. Since v k represents the relative velocity compared to the local rotating motion of the AGN disk, the acceleration by gas accretion torque is Although gas is considered to be captured at the rate of Eq. (29), the accretion rate onto the BH may be smaller than this value due to radiation feedback and the inefficiency of angular momentum transport. When we calculate the type I/II torque from a circumbinary disk, the accretion rate onto the binary is limited by the Eddington accretion rate (Eq. 24).
We adopt the gas accretion rate onto single BHs and binary BHs using Eq. (24). On the other hand, when BHs are in binaries, we apportion the gas accretion rate between the binary components aṡ where λ is the ratio of the accretion rate onto the j th 2 BH over that onto j th 1 BH. We adopt λ given by the fitting formula in Eq. (1) of Kelley et al. (2019), based on the results of earlier hydrodynamical simulations ).

Binary-single interaction
After a binary-single interaction, the binary separation changes depending on the hardness of the binary (e.g. Heggie 1975;Hills 1975;Binney & Tremaine 2008). We make two types of prescriptions for binary-single interaction according to whether a binary is hard, which satisfies is the binding energy of the j th binary.
For soft binary-single interactions, we employ the prescription for the average softening rate derived in Gould (1991), where Λ = v 2 bin,j /σ 2 c is the Coulomb factor, and v bin,j is the relative velocity of components in the j th binary. We set v bin,j = (Gm j /s j ) 1/2 assuming e j = 0. This equation assumes an isotropic Maxwell-Boltzmann distribution for the velocity of the background objects, which is approximately justified for the interaction with the spherical stellar component. We do not account for soft binary-single interactions with the disk BHs and stellar components since most interactions with these components are not soft, due to the low relative velocity (σ c < v Kep h AGN /r Gm j /r Hill,j ).
On the other hand, for hard binary-single interactions, the hardening rate and the kick velocity after a binarysingle encounter are given following Leigh et al. (2018). The binary hardening rate is given by where E be,j and E af,j are the binding energies of the j th binary, i.e. Gm j1 m j2 /(2s j ), before and after the binarysingle interaction, and ∆K j and ∆K c are the changes in the kinetic energy of the j th binary and the escaping third body, respectively. Since ∆v BS,j and ∆v BS,c are on order of v bin,j and v j is typically smaller than v bin,j , we approximate ∆K j = 1 2 m j ∆v 2 BS,j and ∆K c = 1 2 m c ∆v 2 BS,c , where ∆v BS,j and ∆v BS,c are the kick velocities onto the j th binary and the third body, respectively. Due to the conservation of linear momentum, ∆v BS,j = −(m c /m j )∆v BS,c . Here ∆v BS,c is set to be the mode of the probability distribution for the kick velocity, which is determined by the energy and the total angular momentum L for the threebody system. Following Leigh et al. (2018), we sample L 2 uniformly from 0 to (11.5/18) L 2 max and set it to zero for the spherical and the disk components, respectively, where L max is the maximum angular momentum of the three-body system (Eq. 7.27 in Valtonen & Karttunen 2006). The assumption of the low angular momentum for the disk components is due to the low v rel,j , and this assumption increases ∆v BS,c for the disk components by a factor ∼ 1.6 compared to ∆v BS,c for the spherical cases with the typical angular momentum value of (2/9) L 2 max (Leigh et al. 2018). To determine the energy of the threebody system, we set r Hill,j and v rel,j to the initial third body position and velocity relative to the binary center of mass, respectively. Hence, when a hard binary-single interaction occurs, the binary is hardened as and a binary center of mass receives the kinetic energy ∆K j due to a recoil kick. For the interaction with the disk BH component, we randomly choose the third body k which is captured by the AGN disk (h k < h AGN ) and resides in the same cell with the binary j, and set m c = m k . We assign the third body a recoil kick given by ∆K c . Even when the k th object is binary, we treat as interaction with a single object with the mass m k . We assume that the direction of this kick velocity is random and isotropically distributed. Geller et al. (2019) verified that the binary evolution due to encounters calculated by this semi-analytical approach matches the results from direct N -body simulations.
The timescale for the occurrence of a binary-single encounter is given by (e.g. Binney & Tremaine 2008) where represents the typical velocity of the third body relative to the center of mass for the binary j, σ coll is the cross section, v rel,mig,j is the migration velocity of the binary relative to the third body, v shear,j is the shear velocity between the center of mass for the binary j and the third body, p c,j is the fraction of time that the j th object spends within the scale height of each component along its orbit around the SMBH. We set v shear,j = p uni v Kep r Hill,j /r j assuming that the difference between the SMBH and the binary j distance and the third object is p uni r Hill,j , where p uni is a random number uniformly distributed between 0 and 1. In Eq. (39), we define p c,j using For the interaction with the spherical stellar component, we always set p c,j = 1. We set the migration velocity of a third body in the disk BH and stellar components by the average of the migration velocity for BHs in a cell hosting BHs j, while a third body in the spherical stellar component has no migrating motion. Here, n c and σ c are the number density and velocity dispersion of objects in the cell of the binary, including the spherical stellar component, the disk BH component, and the disk stellar component. 2 The cross section is approximated by where b xy,j and b z,j are the effective maximum impact parameters of objects approaching from different directions such that they approach the binary center of mass within a binary separation s j at closest approach. This is given approximately by where the square root term accounts for gravitational focusing by the binary's center of mass calculated in the limit that it dominates over the gravity of the SMBH inside the binary's Hill sphere. We conservatively neglect binary-single interactions with objects that have a larger impact parameter than r Hill,j . 3 In Eq. (45), we also limit the maximum impact parameter due to the lack of objects moving at the elevation above h eff . When the binary is embedded in the AGN disk, it interacts with both the spherical stellar and the disk stellar and BH components, and otherwise interacts only with the spherical stellar component. We reduce the number of stars in the disk stellar component in the l th cell hosting the binary by one for each BH binary that experiences a hard binary-stellar single interaction with objects in the disk stellar component.

Weak gravitational scattering
Weak gravitational scattering is the velocity exchange due to encounters between single objects or between the scattering of the center of mass of a binary with a single object. Based on the Fokker-Planck approximation for an infinite homogeneous medium (Eq. 7.92 in Binney & Tremaine 2008), we assume that the mean acceleration due to weak gravitational scattering is where p c,k sets the fraction of time that the object spends within a scale height for each component, respectively (Eq. 41), the first term corresponds to dynamical friction, and the second term represents Brownian motion , v c is the mean velocity of the medium in the comoving Keplerian frame 6 , erf(X) is the error function, and ln Λ is the Coulomb logarithm. We set Λ = h c /b 90,k following Papaloizou & Larwood (2000), where b 90,k given by Eq. (44) is the impact parameter at which the direction of particles change by 90 degree during the encounter. Eq. (47) is valid when Λ > 1. On the other hand, when h c < b 90,k , the scattering is approximately confined within a two-dimensional plane. In this case, we assume that the direction of the acceleration is along the xy-plane, and set the acceleration according to the Fokker-Planck results for objects confined to an infinite homogeneous 2D medium (Kocsis, Tagawa, Haiman, in prep.) where n xy is a unit vector in a random direction in the xy (AGN disk) plane and the dynamical friction and Brownian motion terms are where x = v 2 c /(4σ 2 c ) and I 0 (·) and I 1 (·) are modified Bessel functions. In practice we resort to the approximations in Eqs. (C1)-(C2) of Appendix § C.
To simply incorporate the effect of the recoil on a third body due to dynamical friction by the disk BH component, we give the same kinetic energy ( k∈l K WS−DF,k /N DBH,l ) to the BHs which are in the AGN disk in the same cell as where K WS−DF,k is the kinetic energy added onto k th objects due to the dynamical friction by scattering with the disk BH component.

Binary formation via three-body interaction
If three bodies make a close encounter, a binary can form dynamically (e.g. Aarseth & Heggie 1976;Binney & Tremaine 2008). This process can be efficient in very dense or low velocity dispersion stellar environments.
For the formation rate of binaries we extend the rough derivation of Eq. (7.111) in Binney & Tremaine (2008) to account for the limitations of a disk geometry and the destructive effect of the SMBH. We assume that binaries can form only when the three bodies undergo a strong encounter. We require to have an encounter between an object i and a second object with impact parameter less than in a region where there is a third object. Here b 90,i is given by (Eq. 44). 7 For a given BH i, the mean rate of strong single-single scattering encounters is given by h eff,i is given by Eq. (46). In each encounter, the probability that a single third BH lies within a distance b i is of order Thus the mean rate for two single bodies to come within a distance b i to the i th BH is (Binney & Tremaine 2008). The factor of 1/3 in Eq. (54) compensates for triple-counting each three-body pair.
To get a rough understanding of the rate of binary formation per BH via three-body encounters, let us evaluate the order of magnitude of Eq. (54), assuming that the mass of objects in the nuclear star cluster is M NSC within an effective radius r NSC , with the density profile in Eq. (1). We approximate the BH density as where γ ρ = 0, f DBH is the fraction of BHs in the AGN disk 8 , v rel ∼ √ 3v Kep (r) h/r and for simplicity we approximated all height-related quantities with the disk scaleheight h assuming that h ≡ h c ∼ h z,i ∼ h eff,i in the radial cell of the object i. Furthermore, we use the empirical M SMBH − σ, M NSC − σ − r NSC relations given below in Eqs. (78), (77), (79) and the number of BHs per unit stellar mass is assumed to be η n,BH = 0.002 M −1 in Eq. (76). This leads to where the four cases are and typically h/r < h AGN /r ∼ 0.01-0.001, Eq. (59) represents the most typical case of the binary formation rate via three-body encounters after the BHs are captured by the AGN disk. Eqs. (57)-(60) show that the binary formation rate is high compared to the AGN lifetime; BHs are expected to form binaries outside of 0.01 pc.
We form a binary due to this dynamical three-body binary formation mechanism in the simulation with probability P 3bbf,i = Γ 3bbf,i ∆t using Eq. (54). We assume that the initial separation of the newly formed binary is s i ≡ b i (Eq. 53). This approximation is justified by the fact that in reality the distribution scales with s 9/2 for s b 90 and it is exponentially suppressed above b 90 (see Eq. 7.175 in Binney & Tremaine 2008).
When the binary formation criteria is satisfied for some BH i, we search for a binary counterpart i in the same cell l. The binary mass is given by m i + m i , and the binary has a velocity of v is the center of mass velocity of the three body system, and v kick = [m c /(m i +m i +m c )](Gm i m i /d)v kick is the kick velocity due to the binary formation, wheren andv kick are unit vectors drawn randomly from the isotropic distribution.
We also consider cases in which the third body is a member of the disk stellar component. In this case, we substitute (n c ) DBH,l (n c ) Dstar,l into n 2 c in Eq. (54), which further increases the binary formation rate compared to Eqs. (57)-(60). We assume that all newly formed binaries are BH-BH binaries, and due to the steep scaling of the formation rate with mass we do not consider the formation of BH-stellar binaries or stellar-stellar binaries in this study. Formation of BH-stellar or stellar-stellar binaries may further increase the rate of BH-BH mergers since such binaries are plausibly exchanged to BH-BH binaries after several binary-single interactions with BHs.
3.3.8. Gas-capture binary formation Generally, if any form of dissipation removes a sufficient amount of energy during the passage of two bodies inside their mutual Hill sphere, the two objects may not travel back to infinity following the encounter, and a binary forms. In the AGN disk, gas dynamical friction serves as a dissipation mechanism. Goldreich et al. (2002) confirmed that the fraction of two bodies passing within their Hill radius which becomes bound (P cap,i , C in Goldreich et al. 2002) is roughly coincident with the fractional decrease in the binding energy during a passage through the mutual Hill sphere (D in Goldreich et al. 2002). Assuming that the latter is approximated as t pass,i /t GDF,i , where is the crossing timescale across the Hill radius, and is the timescale of damping of the relative velocity between two bodies by gas dynamical friction (Eq. 25). Hence the probability for binary formation by the gascapture mechanism during a single BH encounter, is The binary formation timescale due to the gas-capture mechanism is where t −1 enc,i = p c,i n DBH r Hill,i z Hill,i v rel,i defines the rate of objects within the disk black hole component to enter the Hill sphere of BH i, where is the maximum height within the Hill's sphere, h eff,i is given by Eq. (46) Thus the binary formation rate for the i th BH is Γ gas,i = 1 2 t −1 cap,i = 1 2 t −1 enc,i P cap,i , where the 1/2 factor compensates for double-counting each two-body pair. The probability of binary formation for the i th BH within the simulation timestep is P gas,i = ∆t/(2t cap,i ). Since the gas-capture mechanism is caused by gas dynamical friction, we assume that it operates only within the AGN disk ( §3.3.3) and when radiation feedback is inefficient ( §3.3.3). We set the initial separation to be s j = r Hill,i . As in the case of the three-body binary formation, after the binary formation criteria are satisfied for some BH i, we search for a binary counterpart i in the same cell l in the AGN disk. The binary mass is given by m i + m i , and its velocity by v j = v i since the velocity of a counterpart i is dissipated during the passing time. Table 1 Fiducial values for the model parameters

Parameter
Fiducial value Mass of the central SMBH M SMBH = 4 × 10 6 M Gas accretion rate from the outer radiusṀout = 0.1Ṁ Edd Fraction of pre-existing binaries fpre = 0.15 Power-law exponent for the initial density profile for BHs γρ = 1 Parameter setting the initial velocity anisotropy for BHs βv = 0.2 Efficiency of angular momentum transport in the α-disk α SS = 0.1 Stellar mass within 3 pc M star,3pc = 10 7 M Initial mass function slope δ IMF = 2.35 Angular momentum transfer parameter in the outer disk m AM = 0.15 Accretion rate in Eddington units onto stellar-mass BHs Γ Edd,cir = 1 Numerical time-step parameter ηt = 0.1 Number of radial cells storing physical quantities N cell = 120 Maximum and minimum r for the initial BH distribution r out,BH = 3 pc, r in,BH = 10 −4 pc Let us estimate the rough timescale for gas-capture binary formation as in § 3.3.7. The capture fraction P cap,i = min(1, t pass,i /t GDF,i ) is ∼ 1 in the ranges we calculate (10 −4 pc < r < 5 pc) as where ρ gas,6 ≡ ρ gas /(10 6 M pc −3 ) ∝ r −3 1pc in the fiducial model according to Figure 4, and we assume f ∼ 1 in a GDF . The rate of gas-capture binary formation per BH is roughly estimated as where we use h < r Hill . By comparing Eq. (59) with Eq. (69), Γ 3bbf,i ∝ r 1/2 1pc , Γ cap,i ∝ r −1/2 1pc , and Γ cap,i ∼ 30 Γ 3bbf,i at r = 1 pc, we conclude that the rate of gas capture binary formation dominates the rate of threebody binary formation.
If the center of mass of a binary makes a close encounter with a third object, it may undergo a gas-capture interaction to form a hierarchical triple. In such a case, the Kozai-Lidov effect may facilitate the merger of the inner binary (e.g. Silsbee & Tremaine 2017;Liu & Lai 2018). Additionally, perturbations caused by density inhomogeneities may lead to a close binary-single encounter, which may contribute to the hardening of the binary. We conservatively neglect these merger pathways in the simulation. This subject merits future investigation.

Time-step
We use a shared time-step ∆t = η t min k s k ds/dt| GW,k , s k ds/dt| typeII,k , s k ds/dt| GDF,k , t BS,k , t GBS,k , where η t = 0.1 is a constant in the fiducial model, but is varied below.

Merger prescription
Since N -body particles represent stellar-mass BHs, a merger occurs when the separation of a binary becomes smaller than the sum of the innermost stable circular orbits (6Gm j /c 2 ) of the binary components. During a BH merger, the merged remnant recoils due to anisotropic GW emission. The GW recoil kick velocity depends on the mass ratio, spin magnitude and spin direction of the merged BHs (e.g. Baker et al. 2007;Campanelli et al. 2007;Koppitz et al. 2007;Herrmann et al. 2007). In the fiducial model, we assume that the BH spins are zero, and set the recoil velocity to v GW = 8, 830 km/s q 2 (1 − q) (1 + q) 5 (71) whose functional form and normalization are derived by Post-Newtonian predictions and fits to numerical relativity simulations, respectively, where q is the mass ratio of merged BHs (Baker et al. 2007). This formula gives a maximum kick velocity of v GW 160 km/s at q ∼ 0.4. However, since the kick velocity is very sensitive to the BH spins, we also investigate cases with constant kick velocity of v GW = 400, 600, and 1000 km/s in Models 22, 23, and 24, respectively. We add the recoil velocity to the velocity of the merged remnant v i , and the direction is assumed to be random and isotropically distributed. We set the mass loss due to GW radiation at mergers assuming zero spins for BHs as (Tichy & Marronetti 2008), where η = q/(q + 1) 2 is the symmetric mass ratio.
4. RESULTS We start this section by describing an illustrative example of a binary that formed and merged in our fiducial model ( § 4.1). We then present the demography and various physical properties of the entire binary population in the fiducial model ( § 4.2 and 4.3), and discuss how the most important quantities change when model parameters are varied (4.4). The cumulative change in the binary center of mass's kinetic energy due to dynamical friction by gas (dashed blue), gas accretion torque (dashed cyan), binary-single interactions with the disk BH (brown) and disk stellar (orange) components, and weak gravitational scattering with disk BHs (purple), disk stars (pink), and the spherical stellar (gray) components.

Formation and evolution of binaries
We first present the evolution of a binary in the fiducial model (labeled as Model 1 in Table 2, whose parameter values are listed in Table 1) in Figure 5. The binary in this figure forms at 0.74 Myr with an initial separation of 8.9 × 10 −3 pc at r j = 1.8 pc by the gascapture mechanism. The masses of binary components are 9.9 M and 9.1 M . In the early phase, the binary separation decreases significantly due to gas dynamical friction of the AGN disk (blue line in panel (b)). The binary also migrates toward the SMBH due to type I/II torque of the AGN disk (orange line in panel (a)). Affected by gas dynamical friction, accretion torque, and weak gravitational scattering, v j settles to a stochastic equilibrium between heating and damping (i.e. dashed blue, dashed cyan, purple, pink, and gray lines in panel (d)) and v j fluctuates stochastically around the equilibrium. The binary experiences several hard binary-single interactions with a disk BH and star (brown and orange lines in panels (b) and (d)). During the interaction, the binary receives recoil kicks due to the binary-single interaction in a random direction (black and blue lines in panel (c) and brown and orange lines in panel (d)). In some of interactions, this binary is chosen as a third object for the binary-single interaction with other binaries. For simplicity, we assume that this third object receives a recoil kick, even though it is a binary ( §3.3.5). Following a recoil kick, the binary's radial migration is delayed (orange line in panel (c)) since the binary moves out of the AGN disk. Then v j is damped by gas dynamical friction and the accretion torque of the AGN disk (dashed blue and cyan lines infrom the panel (d)). After the binary migrates to r j 6 × 10 −2 pc, and the separation reaches s j ∼ 10 −6 pc ∼ 0.2 AU binary-single interactions become frequent with disk BHs (brown line in panel (b)). After the binary is hardened to s j ∼ 1.8 × 10 −8 pc, GW radiation drives it to merge (purple line in panel (b) of Figure 5) at 7.9 Myr. This binary merges out-  (green), and in-situ formed BHs (gray), and BHs that migrated within the inner radius r in = 10 −4 pc (brown). The number of BHs is significantly reduced by frequent repeated mergers. The reduction of binaries due to disruption (orange dashed curve labeled "ionization") is relatively less important.
side of the AGN disk since the merger takes place soon after a binary-single interaction. During the evolution, 0.8 M and 0.9 M gas mass is accreted onto the primary and secondary BHs, respectively. In summary, this binary, formed via the gas-capture mechanism, it is initially mostly hardened by gas dynamical friction, and later by a series of binary-single interactions with other BHs in the AGN disk, and finally by GW radiation. Figure 6 shows the evolution of the number of BHs (dotted black), mergers (dotted red), all binaries (dotted purple), BHs (cyan) and binaries (green) in the AGN disk, BHs that migrated within the inner radius r in = 10 −4 pc (brown), and in-situ formed BHs (gray) in the fiducial model. Initially the number of binaries is 1.5 × 10 3 , the number of BHs and binaries in the AGN disk are 3.3 × 10 2 and 35, respectively. The number of BHs and binaries in the AGN disk increase for ∼ 3 Myr (cyan and green lines), and then these quantities decrease due to the reduction of the number of BHs due to mergers. The reduction due to binary disruption (orange dashed curve labeled "ionization") is relatively less important. The number of mergers and in-situ formed BHs keep increasing for 100 Myr (dashed red and gray lines). Up to 100 Myr, 1.2 × 10 4 mergers occur and 1.0 × 10 3 BHs are formed in-situ. This implies that some fraction of BHs merged several times within 100 Myr. Figure 7 shows the evolution of mergers among binaries formed by the different mechanisms in the simulation. Mergers among gas-capture binaries, pre-existing binaries, and dynamically formed binaries start to take place from 0.05, 0.2, and 0.3 Myr (solid cyan, purple, and orange lines), respectively. In the fiducial model, mergers among the binaries formed by the gascapture mechanism dominate over the other formation channels (cyan line). This highlights the importance of gas-capture binary formation when discussing compact object mergers in AGN disks. Mergers among in-situ formed BHs contribute only 4.1% of the total number of mergers at 100 Myr (dashed black line). A fraction of 0.85, 0.11, and 0.04 of mergers at 100 Myr are among gas-capture binaries, pre-existing binaries, and dynamically formed binaries (solid cyan, purple, and orange lines), and 0.64, 0.27, and 0.09 of mergers among in-situ formed BHs at 100 Myr are among gas-capture binaries, pre-existing binaries, and dynamically formed binaries (dashed cyan, purple, and orange lines), respectively. Figure 8 shows the distribution of binaries in binary separation (s j ) and distance from the SMBH (r j ). At t = 0 yr (panel (a)), 1462 pre-existing binaries are distributed according to the initial condition. The upper bound on s arise from the assumption that soft binaries are disrupted prior to the AGN phase, and the lower bound on s come from the assumption that the initial binary separation must be larger than the sum of the radii of the progenitor stars to survive without merging during the main sequence stellar phase.

Demography of the BH binary population
Since gas-capture binaries form within the AGN disk, they spend a large fraction of their time within the AGN disk. Both the radial position (r j ) and the binary separation (s j ) evolve simultaneously for such binaries (Figure 5). Figure 9 shows the dominant binary hardening mechanisms at two different times. Gas-capture binaries (cyan circles in Figure 8) are hardened mostly by gas dynamical friction (blue circles in Figure 9) before binaries migrate to r ∼ 10 −2 pc and their separation reduces to r [pc] Figure 8. The distribution of binaries in binary separation and distance from the SMBH at four different times. Cyan circles, purple squares, and orange triangles represent gas-capture binaries, pre-existing binaries, and dynamically formed binaries, respectively. Blue crosses, red pluses, and brown squares represent gas-capture binaries, pre-existing binaries, and dynamically formed binaries for in-situ formed BHs, respectively. Panels  are hardened mostly by type I/II torque from a circumbinary disk (cyan squares in Figure 9). Pre-existing binaries also merge after migrating to ∼ 10 −2 pc (purple line in panel (b) of Figure 10). Panels (a) and (b) of Figure 10 show the number of BHs, stars, and mergers as a function of the distance from the SMBH, respectively. Black, brown, and blue lines in panel (a) are all BHs, and stars and BHs in the AGN disk at 30 Myr, respectively. At r ∼ 0.01 − 0.1 pc, most of BHs are captured in the AGN disk due to the strong gaseous drag from the high-density AGN disk.
Most mergers occur in r 0.01 pc (panel (b) of Figure 10). In this region, the surface density of the AGN disk Σ disk,min is reduced by torques from stellar-mass BHs (Eq. 19) due to the low scale height for the AGN disk (blue line in Figure 4), and the migration timescale increases by a factor of a few (blue line in panel (c) of Figure 10).
Panel (c) of Figure 10 shows the relevant timescales affecting the evolution of a BH binary with m k = 10 M in the simulation. Black and gray lines show the timescale for the initial supersonic velocity of BHs relative to the The total number of mergers for all binaries (black), mergers of gas-capture formed binaries (cyan), pre-existing binaries (purple), and dynamically formed binaries (orange) at 30 Myr. (c): The relevant timescales for an m k = 10 M binary as a function of radial distance from the SMBH. Lines represent the decay timescale of the BH velocity while crossing the disk due to gas dynamical friction (t capAGN , Eq. 26) for v ini,k = 0.2 v Kep (black) and v ini,k = 0.3 v Kep (gray), the timescale for binary hardening by gas dynamical friction to the hard-soft boundary (t GDF,HS , Eq. 27, orange), the migration timescale (Eq. 18, blue), and the typical maximum lifetime of AGN disks (100 Myr, cyan). Binaries migrate to pc before they become hard binaries. . The total binding energy lost to different mechanisms as a function of the distance from the central SMBH at merger for all successfully merged binaries over 30 Myr in the fiducial model (Model 1). Points represent hardening by type I/II torque from a circumbinary disk (cyan), gas dynamical friction (blue), binarysingle interactions with the disk BH (brown) and stellar (orange) and spherical stellar (green) components. Binary-single interactions dominate the total energy lost during the binary hardening.
local AGN motion v ini,k to decay to due to gas dynamical friction and for the BH binary to be captured in the AGN disk, t capAGN , for v ini,k = 0.2 v Kep and v ini,k = 0.3 v Kep (Eq. (26)). The blue line in panel (c) of Figure 10 shows the migration timescale given by Eq. (18). The orange line in panel (c) of Figure 10 is the timescale on which the binary is hardened by gas dynamical friction, t GDF,HS , given by Eq. (27). The cyan line marks 100 Myr, which roughly corresponds to the upper limit for the typical lifetime of AGN disks (e.g. Martini 2004). We do not show the hardening time scale of binaries due to type I/II torque from a circumbinary disk surrounding the binary since it depends on the binary separation (Bartos et al. 2017). For example, the hardening timescale by type I/II torque from a circumbinary disk around a binary with m j1 = m j2 = 5 M and the Eddington accretion rate is ∼ 100 Myr at s j =1 AU and ∼ 200 Myr at s j =0.1 AU. In Figure 9 we have shown that binary hardening by type I/II migration dominates over gas dynamical friction for small binary separations at a fixed location in the AGN disk beyond 0.05 pc, but it is subdominant at nearly all separations for r 0.05 pc. Note that recently revised type I/II migration timescales (e.g. Duffell et al. 2014;Kanagawa et al. 2018) are longer than previously used type II migration timescale in which a clear gap is assumed. At around 1 pc, t mig and t GDF,HS intersect (orange and blue line in panel (c)). Most outer binaries in the AGN disk migrate within this radius before they are hardened by gas dynamical friction. Also outside of ∼1 pc, BHs are not captured by the AGN disk if the initial relative velocity of BHs v ini,k are higher than 0.2 − 0.3 v Kep .

Binary hardening mechanisms
To clarify the relative importance of different binary hardening mechanisms for the merging binary population in the fiducial model, Figure 11 shows the binding energy lost to several mechanisms as a function of the distance from the central SMBH (r j ) at merger for all successfully merging binaries (cf. Figure 9 showing the instantaneous dominant mechanism for all binaries). Most merged binaries are hardened mostly by binary-single interactions with the disk BH component (brown dots in Figure 11). Some fraction of mergers are hardened due to binarysingle interaction by 10 51 erg. The masses of these merged binaries are very high ∼ 10 2 − 10 3 M due to repeated mergers ( Figure 12 below and panel (b) of Figure 10). For comparison note that the separation as a function of binding energy E of a binary with total mass 10m j,10 M and symmetric mass ratio η j = q j (1+q j ) −2 is s j = 10 −6 pc m 2 j,10 (4η j )(E/10 48 erg) −1 , where 10 −6 pc = 0.2 AU.
Gas dynamical friction also contributes ∼ 10 48 − 10 49 erg to binary hardening. As seen in the Figure 9, gas dynamical friction hardens binaries during the early phases of their evolution.
For mergers at r j ∼ 0.01 pc from the SMBH, binarysingle interactions with the disk stellar component contribute 10 46 − 10 48 erg to binary hardening (orange dots in Figure 11). Since stars in the AGN disk are distributed at r 0.01 pc (panel (a) of Figure 10), interactions with the disk stellar component occur before binaries migrate to 0.01 pc. Table 2 shows the results for several different variations with respect to the fiducial model. The input represents the settings for parameters and mechanisms in each model. In the output columns, we show the properties of the system at 30 Myr, namely, the number of merged binaries (N mer ), surviving binaries (N bin ), the fraction of mergers among pre-existing binaries (f mer,pre ), gas-capture formed binaries (f mer,gas ), and dynamically formed binaries (f mer,dyn ) compared to the total number of mergers, the fraction of repeated mergers over total mergers (f mer,rep ), the number of BHs which migrate within the inner boundary r in (N acc ), and the number of mergers among BHs formed in-situ due to the fragmentation of the AGN disk (N mer,SF ).

Number and fraction of mergers
The fraction of the merger number (N mer ) over the initial BH number (N totBH,ini ) (f BH,mer ≡ N mer /N totBH,ini ) at t = 30 Myr is ∼ 0.02 − 0.8 for Models 1-34, where N totBH,ini = 2 × 10 4 in Models 1-14, 20-34 and N totBH,ini = 1.0 × 10 5 , 6.2 × 10 4 , 2.0 × 10 3 , 6.6 × 10 3 , and 6.0 × 10 3 for Models 15, 16, 17, 18, and 19, respectively. We find that the merger fraction is almost unaffected by hardening due to gaseous processes (Models 1, 3-5), the recoil kick velocity at mergers (Models 24-26), or the pre-existing binary fraction (Models 27-29). On the other hand, the merger fraction is lowest when migration does not operate (Model 2), and highest when we only consider the evolution for BHs in the inner regions (r out,BH = 0.3 pc, Model 17). This is because mergers require BHs to migrate to the inner regions of 0.01 pc where the hardening by binary-single interaction is efficient. In the cases of lower β v , which determines the dispersion of the initial BH velocity distribution, the number of mergers is higher (Models 1, 11-14) since a larger fraction of BHs is captured in the AGN disk where binary formation and hardening are efficient. Also the merger fraction is low in the high SMBH mass case (Model 10) because the high Keplerian velocity enhances the decay timescale of the BH velocity v k (see Eq. 26). Note that high-mass SMBHs tends to have larger AGN disks and nuclear star clusters, which also needs to be considered for the estimate of the merger rate (see § 5.5 below).
The number of mergers from in-situ formed BHs (N mer,SF ) depends strongly on the accretion rate of the AGN diskṀ out (Models 1, 7, 8) and the angular momentum transfer parameter m AM (Models 21-23). This is because the star formation rate depends onṀ out and m AM . The rate of mergers among in-situ formed BHs is estimated in § 5.6.

Convergence test
We checked whether the results change due to the time step parameter (Models 1, 33, 34) or the size of cells (Models 1, 31, 32) in which physical quantities for components (Eqs. 6,8,9), and the AGN disk ( Figure 4) are stored. The number of mergers (N mer ) is not significantly affected by these parameters at around the fiducial values (see N mer for Models 1 and 32 or Models 1 and 34). Hence the merger fraction, which we are interested in, is not influenced by the numerical resolution.

Repeated mergers
In the fiducial model, the fraction of repeated mergers is as high as ∼ 0.26 at 30 Myr, allowing massive BHs to form. Figure 12 shows the distribution of merged BHs in the binary mass vs binary mass ratio plane, and Figure 13 shows the cumulative distribution of the merged mass of binaries. Figure 12 and 13 show that BHs can grow to ∼ 10 2 − 10 4 M by repeated mergers in most models. In the fiducial model at 30 (10) Myr, the maximum BH mass is 1.4 × 10 4 M (8.6 × 10 3 M ), and 14 (9) and 1 (2) BHs are more massive than 10 2 M and 10 3 M , respectively. In Model 15 (δ IMF = −1.7) at 30 (10) Myr (panel (k) of Figure 12) in which the initial number of BHs is higher (N totBH,ini = 10 5 ), the maximum BH mass is 2.2 × 10 4 M (1.2 × 10 4 M ), and 128 (81) and 6 (4) BHs are more massive than 10 2 and 10 3 M , respectively. Thus we find that a lot of intermediate mass BHs of 10 2 M are reproduced by repeated mergers. Figure 14 shows the normalized detection rate in the binary mass vs mass ratio plane for several models in Table 2. The normalized detection rate is calculated as the product of the merger rate and the detection volume (e.g. Eq. 6 of The LIGO Scientific Collaboration & The Virgo Collaboration (2012)). We use the noise spectral density from the calibrated sensitivity spectra of aLIGO-Handford on 2018 November 8 (Kissel & Betzwieser 2015). Note that our simulations can only roughly estimate the total binary mass and mass ratio distribution of mergers since we do not take into account the exchange interactions during binary-single and binary-binary interactions. Table 2 Summary of results in different models. The first column indicates the model number. "DF", "type I/II", "max", and "No gas hard." in the "Gas" column represent models in which binaries are hardened by gas dynamical friction, type I/II torque, maximum of gas dynamical friction and type I/II torque, and not hardened by gas interaction, respectively. "No gas mig." and "No mass incr." in the "Gas" column mean models without the type I/II torque from the AGN disk and without the mass increase by gas accretion, respectively. In the "parameter" column, we indicate parameters that deviate from the fiducial model. In the output columns, we summarize the main results in each model: Nmer is the merged number, N bin is the binary number, fmer,pre, fmer,gas and f mer,dyn are, respectively, the fraction of mergers from pre-existing binaries, gas-capture binaries, and dynamically formed binaries compared to total mergers, fmer,rep is the fraction of repeated mergers over total mergers, Nacc is the number of BHs which migrate within r in , and N mer,SF is the merged number from in-situ formed BHs at t = 30 Myr. In Figure 14 panel (a), we have overlaid the mass distribution of the observed LIGO/VIRGO sources during the O1 and O2 observing runs (gray points). Interestingly, despite the fact that the expected maximum BH mass at birth is limited to 15 M due to the solar metallicity environment in galactic nuclei, our results in Figure 14 suggest that the total mass of the detectable merging binaries in AGN extends to masses of 250 M . This is beyond the mass of the detected sources announced to date.

BH binary parameter distributions
The hardening of the mass function of mergers in AGN migration traps, if they exist, was previously pointed out by Yang et al. (2019a). Our results confirm the assumption on the existence of a region similar to a migration at ∼0.01 pc. The hardening of the mass function seen in Figure 12 is more prominent than previously thought, due to the relatively high likelihood of multiple generations of mergers.

Impact of the GW recoil velocity
In the case that the recoil velocity at mergers due to anisotropic GW radiation is 600 km/s (Model 25), 19 (19) BHs grow to 10 2 M by repeated mergers at 30 (10) Myr in a single AGN, which is similar to our fiducial model assuming much lower GW kicks. To realize the high kick velocity of v GW 600 km/s, the BHs need to be rapidly spinning and the directions of BH spins must be misaligned with the angular momentum direction of a binary. Hence even if the recoil kick velocity is very high, repeated mergers build up binaries with total masses of 100 M , which contribute moderately to the detection rate.
Our provide a compelling case that mergers are facilitated by AGN disks, although masses in some cases can be significantly overestimated due to statistical noise fluctuations (Fishbach et al. 2019).
4.4.6. BH growth by gas accretion Figure 15 shows the ratio of the accreted mass over the BH mass for the fiducial model. Gas accretion contributes to the BH masses by less than several tens of percent. Thus gas accretion is not a dominant mechanism for the BH growth. The contribution of gas ac-cretion decreases as the BH mass increases. This suggests that BHs violently grow through repeated mergers. Figure 16 shows the merger number as a function of the generation of primary and secondary BHs for the fiducial model. The generation of a BH is defined as the number of mergers that the BH has experienced in the past plus one. Some BHs have experienced hundreds mergers in their past, during which the masses of BHs are increased due to repeated mergers also by a factor of about hundreds. Hence BHs grow mainly due to repeated mergers. Figure 17 shows the merger number as a function of the v GW =1000 km/s Figure 13. The distribution of binary masses at mergers. The layout of models in the different panels is the same as in Figure 12.
generation of primary BHs for the several generation of secondary BHs. Most mergers are between 1st generation BHs (orange histogram), while the growth of massive BHs (large N g,p ) is dominated by mergers with 2nd generation BHs (cyan histogram).
5. DISCUSSION 5.1. Spin of merging binaries The binary in Figure 5 merges outside of the AGN disk (h z,k > h AGN,l ). In the fiducial model, 81%, 80%, 97% and 19% of all mergers, mergers among gas-capture binaries, pre-existing binaries, and dynamically formed binaries, respectively, occur outside of the AGN disk (evaluated at 30 Myr). The direction of the internal orbital angular momentum of binaries orbiting outside of the AGN disk may be misaligned in a random direction following hard binary-single interactions for kicked binaries.
Similarly vector resonant relaxation and flyby encounters may reorient the orbital plane pre-existing binaries. In these cases, the correlation between the direction of BH spins and the angular momentum of binaries vanishes, which can produce mergers with low values of the effective spin parameter. GW observations may statistically distinguish a population of higher generation mergers based on the measurement of the effective spin (Gerosa & Berti 2017;Fishbach et al. 2017). Note that 91% of mergers among dynamically formed binaries merged in the AGN disk are mergers of 10 3 M BHs, which are not observed by LIGO.
Furthermore, 19% of mergers occur within the AGN disk in the fiducial model. In these mergers, the binary's orbital angular momentum is expected to be aligned or antialigned with the angular momentum direction of the  Figure 12, but showing the normalized detection rate of mergers in the binary mass vs. binary mass ratio plane in several Models. We use the noise spectral density of the ER13 (prior to O3) run of LIGO Hanford (Kissel & Betzwieser 2015). The merger rate is normalized by the maximum value in the plane. In panel (a), the GW events detected to date are overplotted (Abbott et al. 2019a). The mass and mass ratio distribution can well match the observed distribution in some parameter regions. Since our simulations do not take into account the exchange of binary components at interactions, this figure presents rough estimates of the mass ratio. Yang et al. 2019b;McKernan et al. 2019). Also previous studies suggest that a few to ten percent mass increase by gas accretion might be sufficient to align the spin directions of BHs in a binary with the angular momentum direction of a circumbinary disk (e.g. Scheuer & Feiler 1996;Natarajan & Pringle 1998;Ogilvie 1999;Hughes & Blandford 2003;King et al. 2005;Volonteri et al. 2007;Bogdanovic et al. 2007;Lodato & Gerosa 2013), and the angular momentum direction of a circumbinary disk is suggested to be the same as for the AGN disk (e.g. Lubow et al. 1999). GW sources with a moderate to high effective spin such as GW151216, GW170403 and GW170817A (Zackay et al. 2019b;Venumadhav et al. 2019;Zackay et al. 2019a) may represent mergers with multiple generation mergers in AGN disks (see also McKernan et al. 2018;Gayathri et al. 2019). Figure 11 and panel (b) of Figure 10 shows that most mergers occur between 0.001-0.02 pc from the SMBH. Due to the high acceleration in the gravitational field of the SMBH in these regions, the effect of the Doppler acceleration may be detectable in a multiyear detection campaign. The resulting frequency drift (or, change in apparent redshift) allows the binary's distance from the SMBH to be estimated from the LISA GW waveforms (Meiron et al. 2017;Inayoshi et al. 2017b;Wong et al. 2018). Inayoshi et al. (2017b) found that when the projected separation r sin I of a 10-10 M BH binary from a 4 × 10 6 SMBH is smaller than ∼ 0.2 pc (where I is the inclination of the orbital plane of the center of mass for a binary), the strain perturbation by Doppler acceleration is detectable. Wong et al. (2018) estimated that when r sin I is smaller than ∼ 0.01 − 0.03 pc for such binaries, the orbital period and velocity around the SMBH are determined with 10% uncertainties. In the fiducial model of a single AGN until 30 Myr, the expected number of mergers within 0.001, 0.003, 0.01, and 0.03 pc from the SMBH is 0.25%, 1.4%, 45%, and 100%, respectively. Even in the case that migration does not operate (Model 2), the fraction of all mergers within 0.003, 0.01, 0.03, 0.1, and 0.3 pc from the SMBH at 30 Myr is 0.50%, 5.0%, 25%, 78%, and 99%, respectively.

Doppler acceleration and gravitational lensing
Furthermore, for configurations in which the binary's  pected to be observable by future GW instruments, provided that the rates in this channel are sufficiently high for detections.
The binary eccentricity can be enhanced to typically around 0.7 during binary-single interactions if the eccentricity distribution approaches a thermal shape (e.g. Geller et al. 2019). 9 Once the binaries are driven by GWs, the eccentricity decreases (Peters 1964). To highlight an example, consider a 5 + 5 M BH binary, whose initial binary separation is 3 × 10 10 cm, and the initial binary eccentricity is 0.7. This setting corresponds to the typical binding energy (∼ 2 × 10 50 erg) down to which first-generation binaries in an AGN disk can be hardened by binary-single interactions. The orbital frequency of such a binary is 1.1 × 10 −3 Hz and the peak GW frequency is f GW = 1.2×10 −2 Hz (Wen 2003). Since LISA will be able to detect the eccentricity of 10 −3 at f GW ∼ 0.01 Hz (Nishizawa et al. 2016), a nonzero eccentricity is expected to be measurable if BH mergers originate from AGN disks. Following Peters (1964), the binary eccentricity subsequently evolves by GW radia- tion to ∼ 3 × 10 −4 at f GW = 10 Hz. The eccentricity is slightly higher than typical values for in-cluster mergers in globular clusters and somewhat smaller than for GW capture binaries in galactic nuclei (e.g. Rodriguez et al. 2018;Gondán et al. 2018;Zevin et al. 2019). Thus, the value of the eccentricity along with masses and spins may provide distinctive footprints to identify BH binaries in AGN. Figure 18 shows the dominant hardening mechanism as a function of the distance from the SMBH and the peak GW frequency f GW (Wen 2003) at two snapshots, 1 Myr (top panel) and 10 Myr (bottom panel) for the fiducial model assuming an eccentricity of 0.7. Most binaries are hardened predominantly by binary-single interactions in the LISA band at ∼ 10 −2 -10 −4 Hz peak frequencies before GWs drive the binaries to merger. 10 10 Note however, that we use the e = 0 approximation to simulate the evolution of binaries. Binaries may be hardened more efficiently by GWs at at a somewhat lower frequency for higher eccentricity than shown in Figure 18. Thus, LISA may detect the GWs of stellar mass BH and IMBH binaries embedded in AGN where astrophysical environmental effects are still significant (see also Barausse et al. 2015). For these binaries, the mean inspiral rate (df GW /dt, the "chirp") is higher than that 5.5. GW event rates 5.5.1. Stellar-mass BH-BH mergers in AGN Here we estimate the merger rate among stellar-mass BHs in AGN disks. We calculate the merger rate density R sBH for stellar mass BH mergers as

GW phase shift and population distribution
where N BH,cross is the number of BHs crossing AGN disks, t AGN is the average life time of AGN disks, n AGN is the average number density of AGNs in the Universe, f BH,mer is the merger fraction per BH in AGN given by Table 2, and M SMBH,min and M SMBH,max are the minimum and maximum SMBH masses we consider here.
To derive N BH,cross , we assume that the power law exponent of the probability distribution function of the BH orbital radii around the SMBH (Eqs. 1 and 55) is γ ρ = 0 on large scales (100 pc). Within pc scales, this value is motivated by theoretical studies of a relaxed masssegregated cluster (e.g. Hopman & Alexander 2006;Antonini 2014). Outside of the radius of a nuclear cluster, the slope of the one dimensional stellar density is observed to be γ ρ ∼ 0 up to ∼ 100 pc Schödel et al. 2014), although some fluctuation exists at each r. In the case of γ ρ = 0, the number of BHs which exist within r is given by N BH,NSC × r/r eff,NSC , where N BH,NSC is the number of BHs in a nuclear star cluster, and r eff,NSC is the effective (half-mass) radius of a nuclear star cluster. Also in this distribution (γ ρ = 0), the number of BHs crossing the AGN disk for a thermal eccentricity distribution is enhanced by a factor of ∼ten compared to the number of BHs on strictly circular orbits (see Eq. 7 in Bartos et al. 2017). Conservatively neglecting this enhancement due to eccentricity, we assume that the number of BHs crossing the AGN disk is where r AGN is the typical size of the AGN disk. The sizes of the AGN disks are highly uncertain. Radiation hydrodynamical simulations suggest that thin dense gas disks extend to pc scales, which are beyond the dust sublimation radius (e.g. Namekata & Umemura 2016;Wada et al. 2016;Williamson et al. 2019 (see Figure 36 in Burtscher et al. 2013), where L bol is the bolometric luminosity of an AGN disk. We set r AGN = of an isolated binary in vacuum, which leads to an astrophysical GW phase shift for non-stationary sources (see also Yunes et al. 2011;Derdzinski et al. 2019). For nearly stationary GW sources (with respect to the observation time), the f GW -distribution of binaries may be used to measure the residence time that the binaries spend at a particular frequency to infer the underlying astrophysical mechanism driving the evolution of the binary separation (Kocsis & Sesana 2011). r AGN,MW (M SMBH /4 × 10 6 M ) 1/2 assuming a fixed Eddington accretion rate as in Eq. (75), where r AGN,MW is the size of the AGN disk for M SMBH = 4 × 10 6 M .
We substitute the M SMBH dependence of r eff,NSC and N BH,NSC in Eq. (74) using empirical correlations as follows. The number of BHs may be expressed with the stellar mass of the nuclear star cluster as where the parameter η n,BH ∼ 0.002 M −1 represents the number of BHs per unit stellar mass for a Salpeter IMF (but see discussion below for BHs in NSCs). The mass of the nuclear star cluster follows (77) (Scott & Graham 2013), where the velocitiy dispersion of the bulge is given by (Kormendy & Ho 2013) . The radius of the nuclear star clusters in late-type galaxies is expressed as (Georgiev et al. 2016) Thus, r AGN /r eff,NSC in Eq. (74) increases from 0.01 to 0.1 as M SMBH increases from 10 5 to 10 9 M . Following Bartos et al. (2017), we use the log-normal fit to the observed AGN mass function in the local universe (Greene & Ho 2007 This mass function includes low-luminosity AGN with Eddington ratios of f Edd ≡ L/L Edd = 0.01. We include mergers in such low luminosity AGN since the rate of mergers from pre-existing BHs is not very sensitive to the accretion rate onto the SMBH (Models 1, 7, 8). Below f Edd ≈ 0.01 geometrically thin AGN disks transition to geometrically thick advection-dominated accretion flows (Narayan & McClintock 2008). In such low-density flows, it is not obvious if mergers proceed as we modeled in this paper. Here, to be conservative, we consider the rate of mergers only in AGN disks with f Edd 0.01. By integrating Eq. (73) between M SMBH,min = 10 5 M and M SMBH,max = 10 9 M using the assumptions above, we find The relative contribution of the mass ranges M SMBH = 10 5−6 , 10 6−7 , 10 7−8 , 10 8−9 M are 0.96 %, 34 %, 59 % and 6.1 %, respectively. The rate is dominated by M SMBH ≈ 10 7−8 M ; this is because the peak of the AGN mass function is at 10 6.7 M (Greene & Ho 2007) and the number of BHs crossing the AGN disk increases with M 0.83 SMBH . 5.5.2. Uncertainties in the merger rate estimate The merger rate estimate in AGN is parameterized in Eq. (81) with the merger fraction per AGN lifetime f BH,mer /t AGN , the radius of the AGN for a MW-sized galaxy, r AGN,MW , and the number of BHs per unit mass in the nuclear star cluster. The uncertainties in these parameters may be estimated as follows.
The radius of the AGN disk based on mid-infrared observations is given by Eq. (75). We assume that the allowed range of r AGN,MW = 0.06 − 0.2 pc for M SMBH = 4 × 10 6 M , considering that the merger rate is dominated by low-luminosity AGNs with f Edd = 0.01 − 0.1 (Kelly & Shen 2013). Table 2 shows the value of f BH,mer for the different models. For a stationary accretion disk with a fixed accretion rate and Eddington rate, the merger fraction per unit time (f BH,mer /t) is 1.7×10 −8 , 2.0×10 −8 , 1.9×10 −8 , 1.3 × 10 −8 , and 6.2 × 10 −9 yr −1 at 1, 3, 10, 30, and 100 Myr, respectively in the fiducial model. Compared to its value at 30 Myr, f BH,mer /t is higher by a factor of ∼ 1.6 at 3 Myr and lower by a factor of ∼ 2.1 at 100 Myr. Thus the merger fraction is correlated with the lifetime of the AGN disk (Figure 7), which reduces the uncertainties of the merger rate caused by the uncertain AGN disk lifetime. From Table 2, the maximum and minimum f BH,mer /t at 30 Myr are ∼ 2.5 × 10 −8 (Model 17) and ∼ 6.7 × 10 −10 yr −1 (Model 2), respectively. Based on these results, we consider a reasonable range of uncertainty to be f BH,mer /t AGN ∼ 3×10 −10 yr −1 -4×10 −8 yr −1 , so that 0.018 (f BH,mer /0.5)(t AGN /30Myr) −1 2.4 in Eq. (81).
In comparison, f BH,mer /t AGN has been poorly constrained in previous studies, with the allowed range as wide as ∼ 5 × 10 −13 − 3 × 10 −7 .
The uncertainty in f BH,mer comes from the uncertainties of the binary fraction, the capture fraction of BHs by AGN disks, and the merger fraction of binaries. The typical lifetime of AGN disks (t AGN ) also has large uncertainties of 1−100 Myr (e.g. Martini & Weinberg 2001;Haiman & Hui 2001;Marconi et al. 2004;Martini 2004). In the present study, we find that binaries are efficiently formed by gas-capture and dynamical mechanisms, and as a result, the merger fraction is relatively insensitive to the highly uncertain pre-existing binary fraction (Models 1, 27-29). The number of mergers is a factor 5 smaller than in the fiducial model if the initial BH population is isotropic (Model 14 β v = 1, cf. Model 1 with β v = 0.2.).
For a Salpeter initial mass function with δ IMF = −2.35, the number of BHs per unit mass η n,BH is ∼ 0.002 M −1 if we assume that the mass range of stars is 0.1 − 140 M and 20−140 M stars form BHs. On the other hand, the initial mass function in galactic centers is suggested to be top-heavy referring to observational (δ IMF = −1.7 ± 0.2 by Lu et al. 2013) and theoretical studies (Nayakshin et al. 2007). For example, η n,BH is ∼ 0.01 M −1 for δ IMF = −1.7. Furthermore, within parsec regions, numerical studies (Miralda-Escude & Gould 2000;Freitag et al. 2006;Hopman & Alexander 2006;Antonini 2014) show that the number of BHs per unit mass η n,BH is enhanced by a factor of ∼ 2 by mass segregation. We assume that the range of η n,BH is ∼ 0.002 − 0.02 M −1 al-lowing for the possibility of both a top-heavy initial mass function and mass segregation.

Extreme-mass-ratio inspirals
We next predict the rate of extreme-mass-ratio inspirals (EMRIs). EMRIs into SMBHs of ∼ 10 4 −10 7 M are promising targets for future low-frequency GW observations by LISA (Amaro-Seoane et al. 2007). In our simulation, the inner boundary is at 10 −4 pc, which corresponds to ∼ 530 R g for a 4×10 6 M SMBH. The distance from which a 5 M BHs on a circular orbit can migrate to the SMBHs within a Hubble time via GW radiation is ∼ 1000 R g . However, if BHs accumulate at migration traps ( 490 R g ) and repeatedly merge with one another (Bellovary et al. 2016;Secunda et al. 2019;Yang et al. 2019a,b;Gayathri et al. 2019), the number of BHs which can spiral into SMBHs may decrease significantly. Thus not all BHs which migrate within our inner boundary may spiral into the SMBH, and so N acc is an upper limit for the number of EMRIs. According to the results of our simulations (Table 2), N acc ∼ 1 − 5 at 30 Myr. If we assume that N acc is independent ofṀ SMBH and most AGNs are low-luminosity AGN with ∼ 0.01 f Edd (Kelly & Shen 2013), and N acc is roughly proportional to time (brown line in Figure 6), the number of BHs accreted onto the SMBH per accreted gas mass can be approximated as N acc /(Ṁ out t AGN ) ∼ (4−20)×10 −5 M −1 . Then we can roughly estimate the EMRI event rate density as whereρ SMBH is the total mass accretion rate onto all SMBHs in the local Universe, and we adoptρ SMBH = 3×10 3 M Gpc −3 yr −1 (Marconi et al. 2004 Derdzinski et al. 2019). Note that the detectable rate by LISA is enhanced by repeated mergers because of the increase of the detection volume due to the increased mass of migrators. In our simulation, all BHs migrate within the inner boundary R min are single (not binary). This is because binaries either merge or are disrupted by binary-single interactions before they migrate to R min .

Comparison with Previous Works
In this section, we compare our models and results with those in previous works on BH mergers in AGN disks by Bartos et al. (2017) and Stone et al. (2017). A major difference between the present study and esrlier works is that we here model the time-evolving system explicitly, enabling us to follow the formation and destruction of binaries, together with the evolution of their separation, their center-of-mass velocity, and their radial distance from the SMBH consistently. Bartos et al. (2017) considered the capture of binaries due to linear momentum exchange with an AGN disk within 0.01 pc from the SMBH, and the hardening of binaries by gas dynamical friction from the AGN disk and type I/II torques from a circumbinary disk. Bartos et al. (2017) found the merger rate to be ∼ 1.2 Gpc −3 yr −1 . Our study is a more detailed and extended version of Bartos et al. (2017). One difference between Bartos et al. (2017) and this work is the binary fraction. Bartos et al. (2017) considered the evolution of pre-existing binaries, whose fraction is assumed to be 0.3, while our study finds that a large fraction of single BHs captured within AGN disks later form binaries by gas-capture and dynamical mechanisms. A second difference is the assumption of the BH distribution. Bartos et al. (2017) assumed an isotropic velocity distribution, in which the merger rate is lower by a factor of ∼ 6 than with the anisotropic velocity distributions adopted here (see Models 1 and 11-14 in Table 2). Such anisotropic velocity distribution is predicted by theoretical (Kocsis & Tremaine 2011;Szolgyen & Kocsis 2018) and observational studies (Trippe et al. 2008;Yelda et al. 2014;Feldmeier et al. 2014Feldmeier et al. , 2015. Also Bartos et al. (2017) assumed the strongly masssegregated number density of γ ρ = 0.5, while our study assumes γ ρ = −0.5 − 0 (Models 1, 20). A third difference is the assumed size of the AGN disks. Bartos et al. (2017) consider a size of 0.01 pc, compared to the ∼ 0.03 − 0.1 pc (M SMBH /4 × 10 6 M ) 1/2 adopted here. The larger disk sizes are motivated by observations of AGN disks (Burtscher et al. 2013). Stone et al. (2017) considered the evolution of binaries formed in-situ in AGN disks at ∼pc from the SMBH. Both Stone et al. (2017) and our study use the disk model proposed by Thompson et al. (2005), and consider the hardening by type I/II torque and binary-single interactions. The differences are that we treat binarysingle interactions considering the evolution of the distribution of BHs, and we include a new mechanism of gas-capture to form new binaries. Stone et al. (2017) estimated the merger rate from in-situ formed binaries to be ∼ 3 Gpc −3 yr −1 , in which the binary fraction (the binary number over the BH number) and the merger fraction are assumed to be 0.56 and 1, respectively. In our simulation, the fraction of the number of mergers from in-situ formed BHs (N mer,SF ) over the number of in-situ formed BHs (N SF ) is ∼ 0.14 − 0.67 in Models 1, 7, 8, and 21-23. Note that these values are upper limits, since we assume that BHs form immediately at star formation. Following the estimate in Stone et al. (2017), we find the merger rate density for in-situ formed BHs to be wherem BH is the average BH mass, f BH is the mass fraction of BHs over stars,ρ SMBH is the total mass accretion rate onto all SMBHs in the local Universe, and f SF/AGN is the star formation rate within the AGN disk over the accretion rate onto the SMBH. In the estimate above, we use f SF/AGN = 1 andρ = 3 × 10 3 M Gpc −3 (Marconi et al. 2004) as adopted in Stone et al. (2017). Reflecting a top-heavy initial mass function with δ IMF = 1.7 − 2.35 (Lu et al. 2013), we find f BH = 0.016 − 0.092 and m BH = 8.4−9.2 M . Since these differences are relatively small, our estimated rate is consistent with that by Stone et al. (2017). We find that mergers between pre-existing BHs captured by AGN disks (∼ 0.08 − 300 Gpc −3 yr−1) likely dominate the mergers from BHs formed in-situ in AGN disks (∼ 0.7 − 22 Gpc −3 yr−1).

Spatial distribution of surviving binaries
The low-mass X-ray binary candidates found by Hailey et al. (2018) is a critical information to consider for binary formation and evolution in the Galactic center. Interestingly, Hailey et al. (2018) have shown that LMXB candidates are found only within ∼ 1 pc despite the fact that other X-ray sources have been observed out to ∼ 4 pc (Muno et al. 2009). A possible reason for the cutoff at ∼ 1 pc in the LMXB distribution has not been previously explained. Figure 19 shows the spatial distribution of hard binaries at 1, 3, 10, and 20 Myr in the fiducial model. Different colored lines show the distribution of binaries formed by different mechanisms. To select the binaries that would survive until today, we impose the following conditions: ((i)) The timescale of merger by GW radiation is longer than 10 Gyr.
((iii)) The timescale of binary-single interactions with the spherical stellar component at the separation s GW is longer than 10 Gyr, where s GW is the separation from which the binary merges in a Hubble time for the mass of the binary and assuming zero eccentricity (Peters 1964).
Condition ((i)) removes binaries whose separation is smaller than  Figure 19. The final binary distribution as a function of the distance from the SMBH r for the fiducial model at t =1, 3, 10, and 20 Myr for the fiducial model. Cyan, purple, and orange lines are the distribution for gas-capture formed binaries, pre-existing binaries, and dynamically formed binaries, respectively. Thick black line shows the observed distribution of X-ray binary candidates (Hailey et al. 2018) as a function of the projected distance from the Galactic center. These binaries satisfy the following conditions: (i) The timescale of merger by GW radiation is longer than 10 Gyr.
(ii) The binary is hard compared to the spherical stellar component (Gm j 1 m j 2 /(2s j ) > (1/2)mstarv 2 Kep ). (iii) The timescale of binary-single interaction with the spherical stellar component at the separation s GW is longer than 10 Gyr. Gas capture binaries can reproduce the distribution of the observed X-ray binaries. which corresponds to lower regions (s 10 −7 pc) in Figure 8. Condition ((ii)) accounts for the lack of the purple points in the upper left region in panel (a)). Condition ((iii)) produces the inner cut-off at ∼ 0.1 pc. Conditions ((i)), ((i))+((ii)), and ((i))+((ii))+((iii)) respectively exclude 1622, 1709, and 1880 binaries from the set of 2750 binaries at 10 Myr.
The outer cutoff seen in gas-capture binaries (cyan lines in Figure 19) can be explained by the competition between binary hardening due to gas dynamical friction and type I/II migration in the AGN disk, whose timescales are shown by the blue and orange lines in panel (c) of Figure 10, respectively. In Model 1, these timescales become equal at around ∼ pc. At pc, migration is faster than hardening. Thus, before gas capture binaries become hard, they migrate in to pc. The exact position of the outer cutoff for gas-capture binaries is influenced by the scale height of the AGN disk (black line in Figure 4) and the binary mass, which affect the ratio of the migration and hardening rates (Eqs. 17 and 27). As the scale height decreases or the binary mass decreases, the cutoff moves inner regions.
For pre-existing binaries, the outer cut-off is at ∼ 3 pc (purple lines in Figure 19), reflecting our assumed initial BH distribution. Since the BH distribution extends to larger radii (e.g. Freitag et al. 2006;Antonini 2014;Generozov et al. 2018;Schödel et al. 2014;Feldmeier et al. 2014), the cut-off at ∼ 1 pc is difficult to reproduce for pre-existing binaries. Although these binaries may hide the cutoff at ∼ 1 pc, the fraction of the pre-existing binaries in galactic nuclei is poorly constrained. Due to rapid hardening expected during the common envelope phase, the formation of LMXBs is puzzling (Podsiadlowski et al. 2003;Wiktorowicz et al. 2014;Li 2015;Wang et al. 2016). This is especially the case in galactic nuclei where only hard binaries are long-lived, so that surviving the common envelope phase without merging is even more difficult (Stephan et al. 2019).
These considerations suggest gas-capture binaries as a possible origin of LMXBs. Therefore, we compare the number of LMXBs between observations and our models. Including the undetected population of faint sources, the total number of LMXBs in the Galactic Center is estimated to be hundreds to a thousand (Hailey et al. 2018). In our simulation, several hundreds to thousands of binaries survive until the dissipation of the AGN disk (Table 2). Also the number of gas-capture binaries, which would survive until today, at 1, 3, 10, and 20 Myr in Figure 19 are 134, 441, 205, and 96 respectively. Note that we only consider the evolution of BH-BH binaries, leaving the corresponding number of LMXBs uncertain. Since stars outnumber BHs in the disk at 0.1 pc (Figure 10), stars are expected to reside frequently in binaries. However, to estimate the number of LMXBs, we need to consider the fact that massive objects more commonly reside in binaries after binary-single interactions (e.g. Heggie et al. 1996). Also after the dissipation of the AGN disk, the distribution of BHs may evolve due to two-body relaxation whose timescale is uncertain and comparable to the Hubble time (e.g. Merritt 2010;Emami & Loeb 2019). To understand the relation between the binary evolution in the AGN disk and observed LMXBs in detail, estimating the effect of the stellar relaxation processes would be necessary. Also, vector reso-nant relaxation redistributes the z-axis angular momenta of LMXBs (Kocsis & Tremaine 2011;Szolgyen & Kocsis 2018), which is roughly consistent with the 2 dimensional distribution of LMXBs (Hailey et al. 2018).
Finally, the apparent paucity of high-mass X-ray binaries (HMXBs) in the Galactic Center is an interesting additional piece of information for binary formation and evolution. Since the AGN phase in the Milky Way ended long time ago, the lack of HMXBs in our scenario can be explained by the short lifetime of massive stars if these HMXBs formed during such long-past AGN phases.
6. CONCLUSIONS In this paper, we investigated mergers of BHs in accretion disks during the active phase of galactic nuclei. We performed one-dimensional N -body simulations combined with a semi-analytical model. While simplified, this model allows us to incorporate the formation and disruption of binaries, binary-single interactions, weak gravitational scattering with ambient stars, gas dynamical friction in the AGN disk, binary hardening due to type I/II torques from circumbinary disks and gas dynamical friction, migration in the AGN disk, star formation, and gravitational wave radiation. Our main results can be summarized as follows: 1. Gas-capture binaries, which have been neglected in previous studies, contribute ∼ 58 − 97 % of all BH-BH mergers in the AGN disk, with dynamically formed binaries contributing up to ∼ 6 % ( Table 2).
2. After their formation, binaries in the AGN disk are hardened by gas dynamical friction, binary-single interactions with disk stellar and BH components, and finally by gravitational radiation ( Figure 5). Since binaries efficiently migrate to the inner regions of the disk ( 0.01 pc) before they merge, Doppler acceleration due to their center-of-mass motion around the SMBH is expected to be observable in many cases by the future gravitational wave detector LISA ( §5.2).
3. Due to the recoil kick at frequent binary-single interactions, ∼ 80% of mergers occur outside the AGN disk. In this case, the binaries' orbital planes at merger are randomly oriented.
4. On the other hand, ∼ 20% of mergers occur inside the AGN disk. In this case, due to accretion from a circumbinary disk, some degree of (anti-) alignment of the BH spins with the binary's internal orbital angular momentum is expected, as suggested in the high-effective spin merger candidates, GW151216, GW170403 and GW170817A ( §5.1).
5. Due to frequent binary-single interactions during the binary inspiral, the binary eccentricity is thermalized in the LISA band. LISA will detect such highly eccentric binaries ( §5.3).
6. The binary separation is driven predominantly by binary-single scattering interactions in the evolutionary phase preceding the GW-driven merger. The GW frequency is in the LISA band already in the binary-single interaction driven phase which may be identified by a significant GW phase shift for individual sources and/or the frequency distribution of a population of sources ( §5.4).
7. We explicitly compute the binary fraction, the capture fraction to AGN disks, the merger fraction of binaries, and the dependence of the merger fraction on the lifetime of the AGN disk. Accounting for uncertainties in these quantities, we find a volumetric BH merger rate of ∼ 0.02 − 60 Gpc −3 yr −1 , whose uncertainties are reduced by several orders of magnitude compared to prior works ).
8. Due to repeated mergers, this pathway naturally explains "heavy" BHs detected in existing GW observations, even if BHs are born with much smaller ( 15M ) masses. Our model also predicts that mergers of yet more massive BHs ( 10 2 M ) will be detected by LIGO in the future (Figure 14). IMBHs formed during repeated mergers in AGN in most of our models (one exception is Model 2 where gas-driven migration was turned off).
10. The distribution of surviving binaries formed by the gas-capture mechanism can reproduce the spatial distribution of LMXBs observed in the Galactic center, including their outer cutoff at ∼ 1 pc. This cut-off arises from the competition between binary hardening by gas dynamical friction and Type I/II migration in the disk ( §5.7). Binaries migrate inside 1 pc before they are hardened.
In this paper, we employed simplified prescriptions, and ignored the exchange of binary components during binary-single interactions and the evolution in the binary eccentricities, as well as other possibly important processes (see § 2). These issues will be further investigated in the future. Table 3 lists the definition of variables which appear in this paper.

CONVERSION EFFICIENCY
In the disk model of Thompson et al. (2005), the outer regions of the AGN disk are stabilized by radiation pressure and supernovae from in-situ formed stars. The model introduces a parameter defined as the conversion efficiency with which star formation converts mass into radiation, and uses this parameter to calculate the stellar contribution of radiation pressure and the radiative flux. We modify the calcalation of by accounting for the stellar and the AGN lifetime. The radiation flux of main sequence stars formed in the disk averaged over the two disk faces is where A = πr 2 is the disk's area of radius r, t star = 10 Gyr(m star /M )(L star /L ) −1 is the lifetime of a star (Hansen & Kawaler 1994),  The direction of the plane and the angular momentum of the AGN disk ρgas, ngas The mass and number density of gas t, ∆t The elapsed time and the time step in simulations Σ disk The surface density of a gas disk r, r k , r l The distance, the distance of the k th object, and the distance of the geometric center of a cell l from the SMBH h disk The scale height of a gas disk s j , v j , v bin,j The binary separation, the velocity of the center of mass of binary relative to the local AGN motion, and the relative rotation velocity of binary components in the j th BH binary cs The sound velocity of gas v k , v xy,k , v z,k The magnitude of the velocity and the xy and z-direction velocity of the center of mass for the k th object relative to the local AGN motion Ω The angular velocity of a gas disk m k , m j 1 , m j 2 The mass of the k th object, and the primary and secondary masses in the j th binary, respectively t AGN The typical lifetime of AGN disks h k The typical height of orbital motion for the k th object, h k = v z,k r k /v KepṀ out The gas accretion rate from the outer radius rout r Hill,k The Hill radius of the k th object with respect to the SMBH, r Hill = r k (m k /3M SMBH ) 1/3Ṁ

Edd
The Eddington accretion rate v rel,k The typical relative velocity of a third body relative to the center of mass for the k th object βv The parameter setting the initial velocity dispersion for BHs b 90,k The impact parameter of the k th object at which the direction of particles changes by 90 degree after an encounter, b 90,k = G(m k + mc)/v 2 rel,k δ IMF A pawer law exponent for the initial mass function of stars r BHL,k The Bondi-Hoyle-Lyttleton radius for the k th object, r BHL,k = m k G/(v 2 k + c 2 s ) 3/2 M star,3pc The stellar mass within 3 pc from the SMBḢ m BHL,k The Bondi-Hoyle-Lyttleton accretion rate for the k th object (Eq. 29) γρ A power law exponent of the radial profile of the initial BH distribution (Eq. 1) p disk,k , p c,k The time fraction that the k th object spends within the AGN disk and a background component, respectively, p disk,k = 2 π asin(h AGN,l /h k ), p c,k = 2 π asin(hc/h k ) Γ Edd,cir The Eddington accretion rate onto a BH binary from a circumbinary disk ρc, nc The mass and number density of a background component m AM The efficiency of angular momentum transport in outer regions of the AGN disk σc The one dimensional velocity dispersion of a background component v GW The recoil velocity on a merger remnant due to anisotropic GW radiation mc The average mass of a background component fpre The fraction of the number of preexisting binaries over the initial BHs number hc The average orbital height of a background component α SS The α parameter which gives the efficiency of angular momentum transport in standard thin disks ρ AGN,l The mass density of the AGN disk at r l ηt The time step parameter h AGN,l The height of the AGN disk at r l N cell The number of the cells storing physical quantities N totBH,ini , Nmer, N bin , Nacc, N mer,SF The number of initial BHs, mergers, binaries, migrator within r in , mergers among in-situ formed BHs r in,BH , r out,BH The inner and outer boundaries for r within which BHs are initially distributed fmer,pre, fmer,gas, f mer,dyn , fmer,rep The fraction of mergers among preexisting binaries, gas-capture binaries, dynamically formed binaries, and repeated mergers over total mergers r in , rout The inner and outer boundaries for r within which we calculate f BH,mer The number of merger over the initial number of BHs, f BH,mer = Nmer/N ini,BH f BH , f BH,n The fraction of the mass and the number of all stellar-mass BHs over the mass of all stars