Orbital Migration of Interacting Stellar Mass Black Holes in Disks around Supermassive Black Holes

The merger rate of stellar-mass black hole binaries (sBHBs) inferred by the Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) suggests the need for an efficient source of sBHB formation. Active galactic nucleus (AGN) disks are a promising location for the formation of these sBHBs, as well as binaries of other compact objects, because of powerful torques exerted by the gas disk. These gas torques cause orbiting compact objects to migrate toward regions in the disk where inward and outward torques cancel, known as migration traps. We simulate the migration of stellar mass black holes in an example of a model AGN disk, using an augmented N-body code that includes analytic approximations to migration torques, stochastic gravitational forces exerted by turbulent density fluctuations in the disk, and inclination and eccentricity dampening produced by passages through the gas disk, in addition to the standard gravitational forces between objects. We find that sBHBs form rapidly in our model disk as stellar-mass black holes migrate toward the migration trap. These sBHBs are likely to subsequently merge on short timescales. The process continues, leading to the build-up of a population of over-massive stellar-mass black holes. The formation of sBHBs in AGN disks could contribute significantly to the sBHB merger rate inferred by LIGO.


Introduction
The Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) has detected the merger of stellar-mass black holes (sBHs) more massive than those previously inferred from electromagnetic observations in our own Galaxy. Additionally, while isolated binary evolution could potentially account for the high sBH merger rate inferred from LIGO detections, -+ 52.9 27.0 55.6 Gpc −3 yr −1 (Belczynski et al. 2016; The LIGO Scientific Collaboration & The Virgo Collaboration 2018), an additional mechanism of sBH mergers in the local universe would ease several of the assumptions necessary in these models.
It has been suggested that over-massive sBHs are most likely to form in galactic nuclear star clusters (Hopman & Alexander 2006;O'Leary et al. 2009;Antonini & Rasio 2016;Rodriguez et al. 2016). The gas disks in active galactic nuclei (AGNs) are particularly promising locations for the formation and merger of over-massive sBHs. As McKernan et al. (2014McKernan et al. ( , 2018 point out, these gas disks will act to decrease the inclination of intersecting orbiters and harden existing binaries, already making them interesting possible locations for LIGO detections of merging sBHs. The recent discovery of a possible black hole (BH) cusp in the core of our own Galaxy (Bahcall & Wolf 1976;Hailey et al. 2018) lends further weight to this possibility.
Orbiters in a gas disk exchange angular momentum with the surrounding gas, leading to a change in semimajor axis known as migration. Migration of objects embedded within the disk provides opportunities for sBHs to form binaries if they encounter each other at small relative velocities, in particular at far smaller relative velocities than in gas-free star clusters (McKernan et al. 2012Leigh et al. 2018). If a gas disk is locally isothermal, the gas torques cause all isolated orbiters to migrate inward (Goldreich & Tremaine 1979;Ward 1997;Tanaka et al. 2002). However in the more realistic case of a disk with an adiabatic midplane, for some values of the radial density and temperature gradients the torque from the disk can also lead to outward migration (Paardekooper & Mellema 2006). Paardekooper et al. (2010) used analytic arguments and numerical simulations to model the sign and strength of migration, and found that there are regions of gas disks where outward and inward torques cancel out, leading to a region of zero net torque where migration halts. Lyra et al. (2010) showed that such regions of zero net torque, or migration traps, are predicted by standard models of protoplanetary disks, and Horn et al. (2012) showed that the migration of protoplanets toward these migration traps can lead to the rapid collisional build-up of giant planet cores.
While Paardekooper et al. (2010) considered only fully unsaturated torques (where the angular momentum of the corotational region is continuously replenished by viscous mixing, thus continuously driving migration), updated work showed migration rates including saturation (Paardekooper et al. 2011). The basic change due to saturation is twofold. First, only larger orbiters with mass ratio q10 −5 will experience sustained outward migration. For lower masses, the width of the horseshoe region is small enough that diffusion saturates the torques; rapid inward migration occurs for planets outside of a narrow range in mass. Second, inclusion of saturation introduces a mass dependency to both the location and existence of convergence zones (Hellary & Nelson 2012;Coleman & Nelson 2014;Dittkrist et al. 2014). The theory of planet migration continues to be refined, with a dynamic corotation torque found (Paardekooper 2014;Pierens 2015) dependent on the migration rate and viscosity, stemming from an asymmetry in the co-orbital region as the planet moves. This torque can stall inward and boost outward migration, taking planets away from the convergence zone and essentially enlarging the region of outward migration. Its action requires a Shakura-Sunyaev (Shakura & Sunyaev 1973) viscosity parameter α10 −2 for orbiters of mass ratio ≈10 −5 (see the Appendix). These torques were included in N-body calculations by Sasaki & Ebisuzaki (2016), who found that they helped form cores of giant planets. Finally, a heating torque was found by Benítez-Llambay et al. (2015), resulting from the protoplanetʼs accretional luminosity, and found to counteract inward migration. The theory of this heating torque has been further developed by Masset (2017) and Eklund & Masset (2017), showing that it can lead to significant eccentricity and inclination pumping. A torque formula for inclusion in evolutionary simulations has been extracted by Jiménez & Masset (2017). The state of the art in the application of these models for planet population synthesis calculations is discussed in Mordasini et al. (2017).
It is entirely plausible that migration models will undergo significant modifications in the future, driven by advances stemming from the unabated rate of exoplanet discoveries. Yet, some of the differences between AGNs and protoplanetary disks give rise to a pause, first and foremost the fact that the latter are relatively cold and thus poorly ionized, with large swaths not unstable to the magnetorotational instability (MRI) (Blaes & Balbus 1994;Gammie 1996;Wardle 1999;Bai & Stone 2010;Lesur et al. 2014;Lyra & Umurhan 2018). Application of planet migration theory to AGN disks should thus focus on results for high-viscosity and turbulent gas. In this respect, dynamical torques, requiring α10 −2 , should probably not be too relevant (see the Appendix). The heating torque, on the other hand, should also exist for black hole orbiters in AGN disks: even though they do not have a surface to heat via accretional shocks, the accretion disks they develop are hot and luminous and should heat up the surrounding AGN gas. We defer exploring this sBH feedback effect to a future publication.
In this work, as in Horn et al. (2012), we prefer to work with the unsaturated torque because Nelson & Papaloizou (2004), Baruteau & Lin (2010), Uribe et al. (2011), and Baruteau et al. (2011) find that the co-rotational torques in turbulent disks are subject to stochastic turbulent fluctuations that keep the corotational torque unsaturated even in locally isothermal simulations. The result has been corroborated by more recent simulations (Guilet et al. 2013;Uribe et al. 2015;Comins et al. 2016); yet, because they could not resolve the width of the corotational region for smaller objects, saturation remains a possibility if the turbulent fluctuations are strong enough to wipe out their horseshoe turns. McKernan et al. (2012) drew on the work of Lyra et al. (2010) and Horn et al. (2012) to develop a model describing a BH merger hierarchy in the AGN disk. McKernan et al. (2014) explored the consequences of this model and predicted that LIGO should detect gravitational waves from a previously unconsidered population of merging overweight sBHs in AGN disks. Bellovary et al. (2016) explored this analogy applying the Paardekooper et al. (2010) migration torque model to two steady-state analytic supermassive black hole (SMBH) accretion disk models derived by Sirko & Goodman (2003) and Thompson et al. (2005). Bellovary et al. (2016) showed that migration traps do exist in both AGN disk models.
Here we build on Bellovary et al. (2016), by using a modified version of the N-body code described by Sándor et al. (2011) and Horn et al. (2012) which implements several manifestations of the gravity of the gas disk around the SMBH in addition to the standard gravitational forces between particles. The additional effects include migration torques, a stochastic gravitational force exerted by turbulent density fluctuations in the disk, and inclination and eccentricity dampening produced by passages through the gas disk on inclined orbits. In order to explore the dynamical behavior of multiple interacting sBHs approaching a migration trap, we take as an example the migration rates and other disk parameters derived from the analytic AGN disk model of Sirko & Goodman (2003).
Embedded sBHs will migrate toward the migration traps modeled in Bellovary et al. (2016) and, due to this migration, sBHs on prograde orbits encounter each other at low relative velocities. These encounters provide favorable conditions for fast sBHB formation and evolution, resulting in frequent mergers detectable by LIGO. Future constraints from LIGO on this merger channel (e.g., from spins or rates) will allow us to constrain AGN disk physics better than present spectroscopic modeling efforts (see McKernan et al. 2018 for a discussion of which parameters can be best constrained by LIGO).

Methods
In this section we describe in detail our modified N-body simulations. Our simulations neglect forces exerted by sBHs on the gas disk aside from those implicitly modeled by the migration torques, the effects of accretion onto either the central SMBH or orbiting sBHs, and general relativistic effects. We also only consider sBHs on prograde orbits and ignore sBHs on retrograde orbits around the central object. We defer detailed modeling of retrograde objects until the torques on them have been derived from work in progress.

Disk Models
The Sirko & Goodman (2003) model is a modification of the classic Keplerian viscous disk model (Shakura & Sunyaev 1973), with a constant high accretion rate fixed at Eddington ratio 0.5. The disk is assumed to be marginally stable to gravitational fragmentation; however the model does not directly take into account magnetic fields or general relativistic effects. The Sirko & Goodman (2003) model assumes some additional unspecified heating mechanism in the outer disk in order to maintain the stability of the disk and prevent fragmentation. Sirko & Goodman (2003) use the opacity models from Iglesias & Rogers (1996) and Alexander & Ferguson (1994) for the high-and low-temperature regimes, respectively. The inner disk is optically thick due to a high rate of Thompson scattering from electrons produced by the ionization of hydrogen. The intermediate region of the disk has a lower electron density, and is therefore less optically thick and cooler.
We use an SMBH mass of M å =10 8 M e . The total mass of the disk integrated out to 2×10 5 au is 3.7×10 7 M e . The midplane temperature, surface density, scale height, optical depth, and Toomre Q as a function of radius in this model are plotted in Figure 1.

Torque Model
We model the disk torque on the sBHs using the analytical prescription of Paardekooper et al. (2010) which incorporates the effects of non-isothermal co-rotation torques. For the azimuthally isothermal case the normalized torque is a b G G = ---( ) 0.85 0.9 , 1 iso 0 while for the purely adiabatic case the normalized torque is g a b xg G G = ---+ ( ) 0.85 1.7 7.9 . 2 ad 0 The adiabatic index γ=5/3, and the variables α, β, and ξ represent the negative local gradients of density, temperature, and entropy, respectively, and are defined as The torques are normalized by where q is the mass ratio of the migrator to the SMBH, h is the aspect ratio of the disk, and Ω is the rotational velocity. The effective torque is interpolated between the isothermal and adiabatic torque models using where Θ is the ratio of the radiative timescale to the dynamical timescale. Lyra et al. (2010) show that Θ depends on the local disk properties as t ps where c v is the thermodynamic constant at constant volume, σ is the Stefan-Boltzmann constant, and the effective optical depth taken at the midplane is (Hubeny 1990;Kley & Crida 2008) The true optical depth τ is given by where κ is the opacity used in the Sirko & Goodman (2003) models (see Section 2.1). Each component of the torque depends on the local disk gradients of density, temperature, and entropy. These torques are implemented into our N-body code as forces on the particles with vector dependence AGN disks are sufficiently ionized (certainly in the inner regions) that the MRI will drive turbulence. We use a model for turbulence developed by Laughlin & Bodenheimer (1994) and further modified by Ogihara et al. (2007) that gives the gravitational forces exerted by turbulent density fluctuations as = -F ( ) F C , 1 0 turb where C is a scaling factor relating the fraction of the force exerted on the gas by the potential Φ to the force that is exerted by the gas on a migrator embedded in the disk. This fraction is Figure 1. SMBH accretion disk model used in our simulations (Sirko & Goodman 2003). From top to bottom are plotted the midplane temperature T, surface density Σ (in g cm −2 ), disk aspect ratio h (H/r), optical depth τ, and Toomre Q as a function of Schwarzschild radius R s . The top axis represents the translation from Schwarzschild radius to parsecs for a 10 8 M e SMBH.
given as The turbulent potential, Φ, is taken to be the sum of n=200 independent, scaled oscillation modes where ψ is a dimensionless measure of the strength of the turbulent force in comparison to the migration forces (see Section 2.2). It is related to the Shakura & Sunyaev (1973) viscosity parameter α by Baruteau & Lin (2010) as where h is the aspect ratio of the disk and comes from the mode lifetime being set by the speed of sound. In our model h is not constant, but to fix the scaling in Equation (13) we set h=0.05. Magnetohydrodynamic simulations of accretion disks suggest typical values for α of 10 −3 -0.1 (Davis et al. 2010). A value of α=0.01 gives us ψ=4.25×10 −4 . In Equation (12), Λ c,m is one oscillation mode defined as Each oscillation mode is defined by m, an azimuthal wavenumber chosen from a log-normal distribution between 1 and 64, and c denotes the initial center of the perturbation. The position c is given in cylindrical coordinates r c and f c selected from uniform distributions from the inner boundary to the outer boundary of the disk and from 0 to 2π, respectively. The z coordinate is assumed to be small enough to have a negligible effect. Ω c is the Keplerian angular velocity at r c . The mode evolves as a function of = + t t t 0 , where t 0 is the time when the mode comes into existence. The lifetime of the perturbation is p D = ( ) t r mc 2 , 1 5 c s which represents the sound-crossing time for each mode. The radial scale of the perturbation is chosen from a Gaussian distribution and scales as σ=πr c /4m. At the beginning of the simulation there are n=200 modes. When one mode expires another mode is created so that there are always 200 modes. Ogihara et al. (2007) showed that all modes m>6 can be left out of the summation to determine the total potential Φ. We use this simplification in our model and only include Φ perturbations where m<7. Equation (10) is used in our model to calculate the turbulent force on a given migrator at position (r, θ) as a function of the local speed of sound, Keplerian angular velocity, surface density of the gas, and time.
We note that when the net vertical magnetic flux of the disk is not sufficiently large, spiral acoustic waves or even radiation stresses dominate angular momentum transport and accretion power instead of MRI turbulence (Jiang et al. 2017). While the perturbations generated through these mechanisms will not be identical to those produced by MRI turbulence, as modeled above, we anticipate they will have qualitatively the same effect on our simulations (see Section 4).

Eccentricity and Inclination Dampening
Tanaka & Ward (2004) have shown that the gas disk exerts a force on migrators that acts to dampen their orbital eccentricity, e, and inclination, i, leading to the co-planar circularization of orbiters. They give the timescale where m is the mass of the migrator and a is the semimajor axis of the migrator. We follow the timescales given in Cresswell & Nelson (2008) for eccentricity and inclination, respectively: where ò=e/h and l=i/h. The resulting forces acting on these timescales as a function of position and velocity of an orbiting body are wherer andẑ are unit vectors in the r and z directions, respectively.

N-body Code
We use the Bulirsch-Stoer N-body code described by Sándor et al. (2011) that was modified by Horn et al. (2012) to include the additional forces outlined above in Sections 2.2-2.4. The total force acting on each sBH in our simulation is The forces acting from the gas disk, F mig , F damp and F turb , are calculated at the beginning of each Bulirsch-Stoer timestep and not recalculated during the modified midpoint method used to calculate F nbody . However, the Bulirsch-Stoer timestep is a small fraction of the dynamical timescales of the sBHs and is reduced during close encounters. Therefore holding these forces from the gas disk constant throughout each Bulirsch-Stoer timestep does not have a significant effect on the simulations.
Our simulations consider two sBHs to have formed a new sBHB once two conditions have been met. First, they must approach each other within a mutual Hill radius, where m i and m j represent the masses of the two sBHs and r i and r j represent their distances from the SMBH. Second, the relative kinetic energy of the binary, where μ is the reduced mass of the binary, and v rel is the relative velocity between the two sBHs, must be less than the binding energy, Due to the complex and poorly understood interactions between sBHBs and the gas disk within the Hill sphere, for simplicity once a gravitationally bound sBHB forms, our model assumes that it is merged. Indeed it is likely given the conditions of our simulations that all sBHBs will merge within approximately 10-500 yr (Baruteau & Lin 2010), which is a short timescale compared to any dynamical timescales. However, escapes from within a mutual Hill sphere are of course possible. We discuss the merging of sBHBs in our simulations further in Section 5.

Initial sBH Populations
In this section we describe the two models for the initial sBH populations used in our simulations, which are outlined in Table 1. We choose the number of sBHs in each model based on the lower limit of about 10 3 sBHs within 0.1 pc of a SMBH estimated by Antonini (2014) based on the distribution of S-Star orbits around Sgr A å . This estimate is consistent with the population of O(10 4 ) sBHs within 1pc of Sgr A * inferred by Hailey et al. (2018). Assuming sBHs are uniformly distributed throughout the disk, we estimate that around 1% of sBHs in an AGN disk will be within the inner 1000 au (≈0.005 pc). Both of our models therefore include 10 sBHs within roughly 1000 au.
The gravitational wave decay lifetime of an sBH a few hundred astronomical units from an SMBH in a gas-free nucleus is (Peters 1964 Using m 1 =10 8 M e , m 2 =30 M e , e 0 =0.05, and a 0 =650 au as fiducial values used in our runs (see below), gives a decay time of approximately 3.72×10 11 yr. Since this value is several orders of magnitude longer than the run time of our simulations, our models do not include the gravitational wave decay of the orbits of the sBHs around the SMBH.
Our three fiducial models (labeled F1-F3 in Table 1) contain 10 sBHs of uniform masses. This uniform mass distribution is different for each fiducial model, and ranges from 10 M e in F1 to 30 M e in F3. The innermost sBH has an initial semimajor axis of 500 au. The semimajor axis of each successive sBH is separated by 30 R mH from the one before it (see Equation (22)). These initial positions are chosen to create a distribution of sBHs around the migration trap found at roughly 667 au by Bellovary et al. (2016). We note that this initial distribution is somewhat arbitrary; however, these fiducial runs are mainly used as a baseline example to show how sBHs of different masses and initial semimajor axes that are initially not under each other's gravitational influence can migrate to form sBHBs in an AGN disk.
In our second set of models the masses of the sBHs vary in a more physically realistic manner. We draw them from the initial mass function for massive stars given by Kroupa (2002), by drawing from a Pareto power-law probability distribution of sBHs with a probability density where a=1.35, m 0 is a scale factor of 5 M e , and x is a mass that is drawn from the distribution. In our two lower-mass runs, denoted LMA and LMB, a randomly generated mass is rejected if it is greater than 15 M e , so the masses of the sBHs range from 5 to 15 M e . In our two higher-mass runs, denoted HMA and HMB, the mass is allowed to range from 5 to 30 M e . Despite being denoted higher-and lower-mass runs, the total mass of the higher-mass runs does not always exceed that of the lower-mass runs because of random variation. This is the case for LMB, for example, which has the highest total mass of 100 M e . The initial semimajor axes for the sBHs in these models are chosen randomly from a uniform distribution ranging from 300 to 1000 au. We do not use an initial-final mass relation for the sBHs (i.e., Fryer et al. 2012) which would require us to make assumptions of the metallicity and supernova explosion model of our simulations. However, our distribution of initial masses for our sBHs remains similar to what they would be if such a relation had been used.
For all models, the initial eccentricities and inclinations of the sBHs are selected randomly from a Gaussian distribution. The mean value for the initial eccentricity is 0.05, with a standard deviation of 0.02. Selections are made until the value is positive. The mean value for the inclination is 0 with a standard deviation of 0°. 05, and the absolute value of the randomly selected value is used. The initial mean anomaly and pericenter values are chosen randomly from a uniform distribution ranging from 0 to 2π.
For the variable mass models the distance between sBHs is calculated based on the randomly generated positional coordinates. If any two sBHs are within 10 au of each other, a new distribution is generated until no two sBHs are within 10 au of each other.
The masses of the sBHs remain constant over the course of the simulations, i.e., the sBHs are not accreting gas. This is a realistic simplifying assumption based on the Eddingtonlimited accretion rates, which would give a mass doubling time of about 40Myr. Since our simulations are only run for 10Myr and most of the mergers take place within the first few megayears, the additional mass due to accretion is insignificant, to both the mass of the object and the migration rate. Gas accretion onto the sBHs could have a significant effect on the gas disks around the sBHs (i.e., feedback). However, these back-reactions have not been well quantified and so we defer the study of the effects of gas accretion to future work. These models were run for 10 Myr which is within the range of estimated lifetimes for an AGN disk (Haehnelt & Rees 1993;King & Nixon 2015;Schawinski et al. 2015). However, the final orbits of all sBHs in all seven models are established in less than 3 Myr, and these orbits remain stable for the remainder of the run. Over longer periods of time we would expect more sBHs to migrate inwards toward the SMBH from the outer disk. These sBHs may perturb the stable resonant orbiters or sBHs in the migration trap. We defer investigation of this evolution to future work. Figure 2 shows the migration history for runs F1, F2, and F3 in the top, middle, and bottom panels, respectively (see Table 1). In the top two panels of Figure 2 the figures on the left show the migration history from the start of the run to shortly after the final merger. The orbits of the remaining sBHs stay the same until the end of the 10 Myr run. The figures on the right in the top two panels are zoomed-in views of mergers for runs F1 and F2.

Fiducial Model
The bottom left panel of Figure 2 shows the main period of mergers for the F3 run. The orbits of the four remaining sBHs remain the same for over 5 Myr. However, a turbulent mode (see Section 2.3) opens up near the remaining orbiters at around 5.3 Myr, causing the 60 M e sBH to form an sBHB with the 180 M ☉ sBH. This merger is shown in the bottom right panel of Figure 2. The turbulent mode continues to cause the semimajor axes of the orbits of the sBHs in the migration trap to oscillate. The oscillations are more distinctive in the semimajor axis of the 30 M e sBH, because it is significantly less massive than the 240 M e sBH.
In these fiducial runs it is clear that more massive bodies migrate faster toward the migration trap, as expected since the migration torque is proportional to the square of the mass of the orbiter and so the acceleration is linearly proportional to mass. Thus, more massive sBHs reach the migration trap more rapidly. For example, the sBHs in model F3 all reach the migration trap or nearby resonant orbits in roughly 14 kyr, whereas it takes the sBHs in model F1 around 140 kyr. In all cases the last sBHs to reach the migration trap region are the innermost sBHs. These innermost sBHs have the slowest migration rates because within 1000 au of the SMBH the aspect ratio of the disk increases with proximity to the SMBH (see Figure 1). The higher aspect ratio of the inner disk also means that the innermost sBHs will remain on eccentric orbits longer than sBHs since the damping force, F damp,r is inversely proportional to h 4 (see Equations (16)-(19)). Figure 3 shows the growth of sBHs through mergers over time. Massive bodies approaching or reaching the migration trap encounter each other at high rates. Since binaries form at greater rates as sBHs migrate toward the migration trap, the faster migration rate of the more massive bodies leads to faster sBHB formation in the more massive fiducial models. For example, F2 and F3 both have four sBHBs form within the first 10 kyr, whereas it takes nearly 50 kyr for a sBHB to form in F1. Figure 4 shows the eccentricity of all 10 sBHs over the first 200 kyr for the F1 run. While the initial eccentricities of the sBHs' orbits are dampened by the gas within the first 10 kyr, these eccentricities can actually delay sBHB formation at earlier times in our simulations. sBHs that are on eccentric orbits may pass within a Hill radius of each other, but because their orbits have different pericenter phases their relative velocities are great enough that the relative kinetic energy of the two sBHs remains greater than their binding energy (see Section 2.5).
Oscillations in the eccentricity of the orbits of the sBHs that occur later in the run are due to interactions between sBHs. As the sBHs migrate into closer proximity with each other they will be pulled toward each other. This feature can be seen in Figure 2 as little spikes in the semimajor axes of the orbiters. These spikes can be periodic if they occur when two orbiters with similar semimajor axes are in phase with each other. The change in semimajor axis drives the eccentricity of the sBHs. The gas disk will dampen these eccentricities, leading to a decrease in eccentricity until another close pass occurs. These interactions are what cause the oscillations in Figure 4. The eccentricity of the sBH orbits rarely increases to more than 10 −2 . This feature is common in our simulations and is discussed further in Section 4.2.
When sBHs pass within just a couple of Hill radii of each other, whether on not their relative kinetic energy is low enough to form an sBHB, the effective semimajor axis of their orbits around the SMBH often spike dramatically as their orbits are strongly perturbed from Keplerian orbits, as can be seen in Figure 2. However, this should be interpreted as a dramatic change in velocity rather than position.
In our runs the most massive sBH consistently ends up closest to the migration trap. However, in some cases, such as in the F2 run, no sBH ends up precisely in the migration trap. Instead the most massive sBH ended up roughly 2.5au away from the migration trap. At these small distances, the migration torque is very minimal, and the dynamics due to the high density of sBHs in the region play a larger role in determining the orbits' positions. Less massive sBHs end up either on Trojan or resonant orbits that exchange angular momentum with the other sBHs. These final configurations tend to be stable on megayear timescales. However, MRI turbulence can lead to sBHB formation even after these stable orbits are established if it knocks a sBH out of resonance, as happened in the F3 run. Encounters with other objects either being ground down into the disk, or migrating inward from further out in the disk might also disturb the steady configurations over longer timescales.

Varying Masses
Figures 5-8 show the migration histories for runs LMA, LMB, HMA, and HMB, respectively (see Table 1). The top panels of all four figures show the migration histories of the simulation up until all sBHs have reached stable orbits. In each simulation the sBHs remain at these final radii for the remainder of the 10 Myr run. The bottom two panels of these figures show examples of sBH interactions.
As the more massive sBHs migrate through the disk they overtake less massive sBHs and frequently form sBHBs. The time that elapses before the first binary capture of the simulation varies among the four runs from a few hundred years to roughly 20 kyr due to the randomly generated initial positions and eccentricities. Even if two sBHs have similar initial positions at the beginning of a simulation, if the orbits of the sBHs are too eccentric, the relative kinetic energy of the two sBHs that approach each other within 1 R mH may be higher than their binding energy, preventing them from forming an sBHB (see Section 4.1). Figure 9 shows the build-up in mass of the sBHs due to mergers for runs LMA (top left), LMB (top right), HMA (bottom left), and HMB (bottom right). In our simulation two sBHs are considered merged as soon as they form an sBHB (i.e., approach each other within 1 R mH ; see Section 2.5). Six to eight mergers occur in each run. The most massive sBH at the end of each run ranges from 45 to 65 M e , which represents 60%-65% of the total mass of the run.
The time that elapses before all sBHs reach the migration trap also varies and depends on the random generation of positions and masses. The smaller the initial semimajor axis of an sBH the longer it will take to migrate toward the trap, especially if it has a smaller initial mass.
Dynamical effects can produce some exceptions. For example, in the LMB run (see Figure 6), there are three sBHs with very small initial semimajor axes ranging from 310 to 320 au. Being in such close initial proximity causes the sBHs to interact with each other from the start, but they do not immediately form an sBHB. The least massive sBH only has a mass of 5 M e and after the three-body interaction ends up on its own at approximately 300 au. This low-mass sBH left alone in a region with a very low migration rate takes nearly 3 Myr to finally make it to the migration trap. The two more massive sBHs (8 M e and 14 M e ) end up in a stable horseshoe co-orbit as modeled by Cresswell & Nelson (2006), who found that it was common for planets in a protoplanetary disk to become coorbital, occupying either horseshoe or tadpole orbits that survived for the duration of their runs. Figure 10 shows the relative phase, semimajor axes, and ratio of the orbital period around the SMBH for these two co-orbital sBHs. Over a period of thousands of orbits the phase difference between the two sBHs oscillates between 180°and 20°. When the phase difference is at a minimum the two sBHs swap radial positions. Occasionally the migration rate of the more massive, 14 M e , sBH is large enough compared to the migration rate of the less massive, 8 M e , sBH that it overtakes it while the two are out of phase. However, the two sBHs still swap radial positions when they are closest to being in phase. As a result the 8 M e sBH migrates at the rate of the 14 M e sBH, which means the 8 M e sBH reaches the migration trap at nearly double the rate it would alone.
In the HMA run, the cyan and purple sBHs in the center panel of Figure 7 are also on a horseshoe co-orbit until the orbit is destabilized by the presence of a 26 M e sBH, which the cyan sBH merges with. The co-orbital tadpole (i.e., Trojan) orbits that were observed by Cresswell & Nelson (2006) are seen in runs F2 and F3 (see Figure 2).
In all cases, after several hundred kiloyears one sBH becomes massive enough to dominate the region closest to the migration trap and lock all other less massive sBHs in highorder resonant orbits. sBHs migrating toward the trap at later times will either merge with the sBHs already populating resonant orbits (Figure 8), or end up on their own resonant orbit (Figure 7). Figure 11 shows the semimajor axes (upper panels) and eccentricities (lower panels) of the sBHs in or near the migration trap for the LMA (top left), LMB (top right), HMA (bottom left), and HMB (bottom right) runs. As in the F2 run (see Section 4.1), in the LMA, LMB, and HMB runs no sBH ends up exactly in the migration trap. Instead the most massive sBH ends up 1-2.5 au from the migration trap, where it becomes locked in a resonant orbit with the other sBHs. The semimajor axes of the sBHs around the SMBH spike periodically as the sBHs on resonant orbits exchange angular momentum with each other and the sBHs in the migration trap get pushed back into resonance. The sudden change in the orbit's semimajor axis causes a spike in eccentricity which is then dampened by the gas. Figure 12 shows one example from the HMB run of these interactions of two sBHs on a 27:28 resonance. When the phase difference between the sBH in the migration trap and the sBH on a resonant orbit is zero, the two are pulled toward each other by their mutual gravitational attraction. This temporarily drives an increase in the eccentricity of their orbits, before it is gradually dampened once again by the gas disk.
These orbits remain stable for 9 Myr to the end of runs LMA, HMA, and HMB, suggesting that trapping sBHs in resonant orbits around a migration trap could prevent more massive sBHs from building up. However in the LMB run, as in the F3 run (see Section 4.1), a perturbative force caused by disk turbulence pushes the 23 M e sBH out of resonance so that it merges with the 42 M e sBH in the migration trap. Therefore disk turbulence could provide a mechanism to break resonances, and create more massive sBHs. Horn et al. (2012) showed that increasing levels of disk turbulence makes this mechanism even more efficient. Batygin & Adams (2017) worked out an analytic solution for the breaking of resonances by turbulence for protoplanetary disks and found that the disruption of resonances by turbulence depends most strongly on the migrator-central mass ratio. For the migrator-central mass ratios and other relevant parameters in our simulations, their analytic solution agrees with our conclusion that turbulence could play a role in disrupting resonances.
The initial inclinations of the sBHs were very small, and all sBHs were quickly ground down into flat orbits in less than 50 yr. The initial eccentricities played a role in our models in preventing early sBHB formation, but were also a transient effect and were dampened by the gas in roughly 10 kyr. Larger initial values for inclination and eccentricity would likely delay sBHB formation because they would increase the relative kinetic energy of two sBHs. However, these larger inclinations and eccentricities will eventually be dampened by the gas disk, and as sBHs are ground down into the disk and their orbits are circularized, they would start to form sBHBs with other sBHs at later times.

Summary and Discussion
We have simulated the migration of compact objects in a model AGN disk (Sirko & Goodman 2003), using an analytic Figure 4. Eccentricities over the first 200 kyr of all 10 sBHs in the F1 run. Initial eccentricities are quickly dampened by the gas in less than 10 kyr; however, interactions between the sBHs in close proximity drive the eccentricity of the sBHs' orbits as they are pulled toward each other. This eccentricity is then dampened until another close passage occurs. model developed from simulations of the migration of protoplanets in protoplanetary disks. We have found that migration due to gas torques in AGN disks can provide an efficient mechanism to create a population of hard compact object binaries remarkably quickly, replicating the results of Horn et al. (2012) for protoplanets in a protostellar disk, but for the case of sBHs in an AGN disk.  . Migration in the LMB run, with the same notation as Figure 5. The top figure shows the first 4.5 Myr, which is the period during which binary formation occurs, and all sBHs migrate toward the migration trap to stable orbits where they remain for the rest of the 10 Myr run. The middle figure is a zoomed-in view of the first period of binary formation and the bottom figure is a zoomed-in view of the interaction between the three innermost sBHs, two of which end up co-orbital. Figure 7. Migration in the HMA run, with the same notation as Figure 5. The top figure shows the first 600 kyr which is roughly the period during which binary formation occurs, and all sBHs migrate toward the migration trap to stable orbits where they remain for the rest of the 10 Myr run. The middle figure shows a zoomed-in view of binary formation that breaks apart two co-orbital sBHs and the bottom figure shows a zoomed-in view of a later period of binary formation. In the top panel the 6 M e sBH is the last to reach the region of the migration trap, because it has a small initial semimajor axis. When it reaches the trap it ends up on its own resonant orbit, instead of merging with other sBHs. where N BH is the number of sBHs in an AGN disk, N GN is the average number density of galactic nuclei in the universe, f AGN is the fraction of galactic nuclei with AGNs that last for time τ AGN , f d is the fraction of sBHs that end up in the AGN disk, f b is the fraction of sBHs that form binaries, and ò represents the fractional change in N BH over one full AGN duty cycle. Using our finding that, within the inner 1000au of an AGN disk 60%-80% of sBHs form sBHBs in the lifetime of our AGN disk, we can use 0.6-0.8 as an upper limit on f b , giving an upper limit on the merger rate of 72 Gpc −3 yr −1 . This value is an upper limit because, although our model assumes a uniform distribution of sBHs throughout the disk, sBHs in the outer disk, further from the migration trap, may merge less frequently. In addition, this upper limit assumes that sBHs orbiting in the retrograde direction would have similar sBHB formation rates, which is unlikely because migration torques on retrograde orbiters should be much weaker. We defer a more realistic prediction of retrograde orbiter merger rates and merger rates of sBHs in the outer disk to future work. Uncertainties in AGN disk structure result in a wide variety of plausible theoretical models to describe these disks. However, migration traps should occur in any disk where there is a rapid change in the surface density gradient (Bellovary et al. 2016). Such rapid changes are likely to occur in most actual disks, since radiation pressure is expected to inflate the inner disk. This paper is intended to highlight the qualitative behavior of objects at the migration trap. Regardless of the location of these migration traps, whether they are at 331R s as in the Sirko & Goodman (2003) model or about 225R s as in the Thompson et al. (2005) model (Bellovary et al. 2016), most of the binary formation will take place in the immediate vicinity of the migration traps. We expect the qualitative behavior around the migration trap to be similar regardless of AGN disk model, although having a migration trap at a different distance from the SMBH as in Thompson et al. (2005) will affect how long it will take sBHs to migrate to the trap. Additionally, different disk models have different surface densities, which will affect migration rates. If these surface densities are lower than in Sirko & Goodman (2003), as they are in Thompson et al. (2005), the migration rates will be lower. We defer simulations of migration in alternative AGN disk models to future work.
We highlight that, although we have taken the compact objects to be sBHs here, similar results apply to any objects embedded in the AGN disk, including neutron stars, white dwarfs, or main-sequence or evolved stars, although their typically lower masses will result in slower migration rates. Our demonstration of how quickly binaries can form in AGN disks may help us to understand the behavior of other objects embedded in AGN disks. For example, Davies et al. (1998) attributed the observed lack of red giant stars in the Galactic center to direct collisions during single-binary encounters. We might suggest a simple alternative, albeit analogous mechanism motivated by our results in this paper: main-sequence turn-off stars efficiently form (or are exchanged into) compact binaries, such that they form common envelope binaries (or some other variation of the myriad possible binary evolution pathways) when the turn-off star evolves up the giant branch, preventing it from evolving into a normal red giant star. In short, myriad binary and stellar exotica could form in AGN disks. These additional compact objects could also contribute non-negligibly to subsequent binary mergers and interactions (Leigh et al. 2016), and even produce exotic populations that might contribute to the total light distribution in galactic nuclei nonnegligibly, once the gas disk has dissipated and the SMBH is no longer actively accreting at high rates.
One assumption of our model is that sBHBs merge as soon as they form. These binaries actually harden due to gas torques on a timescale that depends on the distribution of gas in the Hill sphere of the binary, and which also involves the complicated effects of accretion onto the sBHB and the resulting feedback. We justify our assumption by comparing the migration timescale to the binary hardening timescale. Baruteau et al. (2011) modeled the hardening of binaries in a gas disk. Their models showed that it takes roughly 1000 orbits of binary stars around the binary's center of mass to reduce the semimajor axis of the binary by a factor of two if the binary is rotating in the prograde direction with respect to its orbit around the central mass, and only 200 orbits for retrograde rotation.
We assume that after the binary's semimajor axis has been halved 20 times, the sBHB separation is small enough that gravitational radiation will rapidly merge the sBHB to form a single sBH of mass m i + m j . The binary inspiral time due to gravitational wave emission alone (Peters 1964), neglecting any gas hardening effects, exceeds the binary hardening timescale of (4-200)×10 3 orbits as long as the binary eccentricity e<0.9995. Note that this estimate may be a significant underestimate of the actual time to merger, since gas hardening may become less efficient as the binary shrinks. However, we have also neglected the possibility of hardening encounters due to tertiary objects in the disk, which will accelerate the rate of binary hardening (Leigh et al. 2016(Leigh et al. , 2018. Both of these complications will require further study in future work. Given our assumptions, Figure 13 shows the approximate radial dependence of the timescales of mergers for sBHBs rotating in prograde and retrograde directions, for sBHBs orbiting in the prograde direction through the disk. In our simulations these timescales will be equivalent for all mass sBHBs because sBHs are considered to form an sBHB when they approach each other within a mutual Hill radius, which is µ + ( ) m m i j 1 3 . For comparison, in Figure 13 we plot the time for 10 M e , 20 M e and 50 M e sBHs to migrate from their current radial location to the SMBH due to migration torques. Recall that the migration torques vary as a function of radius and temperature, surface density, and disk aspect ratio at each radius. Since the migration timescales of these objects are at least an order of magnitude larger than the time it would take for an sBHB of the same mass to merge, we can see that the likelihood of a tertiary encounter from another sBH is low.
This low likelihood is important because, while a tertiary encounter could accelerate a binary merger (Leigh et al. 2018), the third sBH would be ejected in the process, and because the sBHB will already be merged in our simulation, it is not possible for a third body to gain energy from a three-body encounter. However, our simulations do permit binary formation to occur via three-body interactions in a limited set of realistic circumstances. That is, three initially isolated sBHs Figure 10. Two sBHs from the LMB run in a stable horseshoe co-orbital configuration. The top panel shows the relative phase between the sBHs, the middle panel shows the semimajor axes of the two sBHs, which are labeled by their current masses, and the bottom panel shows the ratio of orbital periods around the SMBH. could end up in a sufficiently small volume that their mutual gravitational attraction dominates locally, and a chaotic threebody interaction ensues. If one star is ejected, the other two remaining sBHs could form a binary.
The dissipative effects of the gas actually enhance the probability of such three-body-mediated binary formation occurring. The critical orbital separation of an sBHB for which the kinetic energy of a third isolated sBH is equal to the orbital energy of the sBHB is known as the hard-soft boundary. Third-body encounters with hard binaries promote hardening, while with soft binaries they can promote ionization. In an AGN disk the hardsoft boundary for a sBHB in a circular orbit is (Leigh et al. 2018 where R H is the Hill radius, μ b is the reduced mass M 1 M 2 /(M 1 + M 2 ) of the binary, and M 3 is the mass of the third sBH. Since we consider sBHBs to be merged once they are within a Hill radius, as long as 12μ b /M 3 >1, the kinetic energy from a prograde tertiary sBH should not be enough to ionize an sBHB in our simulations.
Looking at examples in our simulations of sBHBs that have a close encounter with a third sBH on timescales shorter than the merger timescales in Figure 13, we find only one instance where a third sBH is massive enough that a HS,disk <R H . However in this case the third sBH never approaches closer than 10 R H from the sBHB, making it too distant to ionize the sBHB. Therefore in our models the ionization of our binaries by a prograde third body interaction appears to be rare.
Finally it is possible for one of two sBHs within a mutual Hill radius of each other to be ejected, even if the binding energy of the two sBHs is less than their relative kinetic energy and there is no tertiary interaction. Preliminary results from A. Secunda et al. (2019, in preparation) suggest that this is rare. In A. Secunda et al.
(2019, in preparation) the merger boundary is reset to 0.65R mH for the same runs as in this paper. This boundary was chosen to allow us to study some of the properties of the sBHBs we were forming without requiring unreasonably large computational resources, in the form of integration time. In four out of 32 cases, sBHs that would have merged under the criteria presented in this paper did not merge when the boundary for merger was 0.65R mH . Instead Figure 12. Zoomed-in view of the interactions between two sBHs from the HMB run on resonant orbits. The top plot shows the semimajor axis of the two sBHs, the middle plot shows the phase difference between the two sBHs, and the bottom plot shows the eccentricity of the two sBHs. When the phase difference reaches zero (represented by the vertical dashed line) the two sBHs are pulled toward each other by gravity, making their orbits more eccentric. This eccentricity is then dampened until they are pulled toward each other again. Figure 13. Various timescales in years are plotted as a function of radial distance from the SMBH in astronomical units. The blue, yellow, and green lines represent the approximate time for 10 M e , 20 M e and 50 M e sBHs to migrate from their current location to the SMBH due to only migration torques. The dashed and dotted black lines represent the merger time for an sBHB that forms when two sBHs are within a mutual Hill radius (see Equation (22)) for a prograde orbiting sBHB and a retrograde orbiting sBHB, respectively. The merger timescales are significantly shorter than the migration timescales, suggesting that the probability of the sBHB failing to merge due to an encounter with a tertiary body is low. these sBHs swapped orbits. This orbit swapping in place of sBHB formation is already seen in runs with the sBHB formation criterion set to 1 R mH (see the yellow and brown lines in the center panel of Figure 6). Additionally, those sBHs that failed to merge initially later were able to merge with other sBHs. Therefore, the merger histories of runs with a more stringent merger criterion were qualitatively identical to those presented above.
Future work that includes the relevant gas physics should evolve sBHBs to much smaller merger boundaries of 0.1-0.01R mH to further probe the poorly understood evolution of sBHBs in a gas disk. We do not do so here for simplicity's sake, since whether an sBHB in a gaseous accretion disk will be able to merge is still an open question. For example, recent work by Moody et al. (2019) found that the circumbinary disk around BHBs can actually exert a net positive torque on the BHB, causing its semimajor axis to increase. The properties of the sBHBs formed in our simulations should serve as useful, physically motivated inputs for future hydrodynamic simulations of BHB evolution in gas disks.
Our model is efficient at building up massive sBHBs on timescales far shorter than the lifetime of the AGN disks that host them. We argue that these sBHBs are likely to merge, producing gravitational wave events such as those observed by LIGO. However, future work using hydrodynamic simulations is needed to better describe the interactions between the gas disk and the sBHs, in particular examining the binary hardening timescale due to gas torques. More work is also needed to model the evolution of the sBH population as additional compact objects either drift inward or have their orbital inclination ground down into the inner region of the disk where they may be able to break resonances and form additional sBHBs and more massive sBHs. We have also completely ignored the role of retrograde orbiters in this paper, and the population of objects on retrograde orbits that can ionize binaries embedded in the disk. General relativistic effects are also not included in our model. Ultimately a full, three-dimensional, time-evolving AGN disk model should be used to provide the most accurate predictions for the merger rates of sBHs and the build-up of over-massive sBHs. In the meantime, constraints from the next few LIGO runs on mergers from this model channel should help put limits on models of AGN disks (in particular, the presence or absence of density gradients likely to produce migration traps), such as those used here.