Orbital Migration of Interacting Stellar Mass Black Holes in Disks around Supermassive Black Holes II. Spins and Incoming Objects

The mass, rate, and spins of merging stellar-mass binary black holes (BBHs) detected by aLIGO and Advanced Virgo provide challenges to traditional BBH formation and merger scenarios. An active galactic nucleus (AGN) disk provides a promising additional merger channel, because of the powerful influence of the gas that drives orbital evolution, makes encounters dissipative, and leads to migration. Previous work showed that stellar mass black holes (sBHs) in an AGN disk migrate to regions of the disk, known as migration traps, where positive and negative gas torques cancel out, leading to frequent BBH formation. Here we build on that work by simulating the evolution of additional sBHs that enter the inner disk by either migration or inclination reduction. We also examine whether the BBHs formed in our models have retrograde or prograde orbits around their centers of mass with respect to the disk, determining the orientation, relative to the disk, of the spin of the merged BBHs. We find that orbiters entering the inner disk form BBHs, that are commonly asymmetric in mass, with sBHs on resonant orbits near the migration trap. When these black holes reach 80 solar masses, they form BBHs with sBHs in the migration trap, which over 10 Myr reach about 1000 solar masses. We also find that 68% of the BBHs in our simulation orbit in the retrograde direction, which implies our merger channel could easily produce the low values of the dimensionless aligned spin, chi_eff, observed in a majority of BBH mergers detected thus far.


INTRODUCTION
The high rate of stellar mass binary black hole (BBH) mergers inferred from Advanced Laser Interferometer Gravitational Wave Observatory (aLIGO) (Aasi et al. Corresponding author: Amy Secunda asecunda@princeton.edu 2015) and Advanced Virgo (Acernese et al. 2014) detections, as well as the high masses and low spins of many of these mergers (Abbott et al. 2019), has inspired much debate over which merger channels could best produce these mergers. Isolated binary evolution (Belczynski et al. 2016(Belczynski et al. , 2010Postnov & Yungelson 2014;De Mink & Mandel 2016) and dynamical formation in globular clusters (Benacquista & Downing 2013;Wang et al. 2016), are two such potential merger channels. Conversely, Hopman & Alexander (2006), O'Leary et al. (2009), Antonini & Rasio (2016), and Rodriguez et al. (2016) suggest that over-massive stellar mass black holes (sBHs) are most likely to form in galactic nuclear star clusters. The apparent sBH cusp at the center of the Milky Way observed by Hailey et al. (2018) lends further weight to this possibility.
Here we focus on the scenario proposed by McKernan et al. (2014McKernan et al. ( , 2018, who have suggested that the gas disks in active galactic nuclei (AGN) are especially favorable locations for the BBH mergers detectable by aLIGO and Virgo. McKernan et al. (2014) point out that gas disks act to decrease the inclination of intersecting orbiters and harden existing binaries. Additionally, the orbiters in a gas disk exchange angular momentum with the surrounding gas leading to a change in the semi-major axes of their orbits, known as migration (Goldreich & Tremaine 1979). As sBH orbiters in the disk migrate with speed dependent on their mass, they encounter each other and form BBHs, especially since sBHs orbiting in the same direction will encounter each other at relative velocities far smaller than in a gas-free star cluster (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 & Ward 2004). However, Paardekooper & Mellema (2006) showed that in the more realistic case of a radiatively inefficient disk with an adiabatic midplane, for some values of the radial density and temperature gradient the torque from the disk can instead cause outward migration.  used analytic arguments and numerical simulations to model the sign and strength of migration. They found that at the boundaries between regions of inward and outward migration the torques cancel out, leading to an orbit with zero net torque where the migration halts (Lyra et al. 2010), referred to as a migration trap. Bellovary et al. (2016) applied the  migration torque model to two steady-state AGN disk models derived by Sirko & Goodman (2003) and Thompson et al. (2005) and showed that migration traps were present in both.
Building on this work, Secunda et al. (2019), hereafter Paper I, used an augmented N-body code (Sándor et al. 2011;Horn et al. 2012) that incorporates the  migration torques to simulate the migration of ten sBHs in a Sirko & Goodman (2003) AGN disk. Paper I found that migrating objects encountered each other at a rapid rate, causing 60-80% of sBHs in their simulations to form BBHs over 10 Myr. From these simulations Paper I estimated an upper limit of the merger rate parameterized in McKernan et al. (2018) of 72 Gpc −3 yr −1 . The most massive models in Paper I produced over-massive sBHs of around 100 M . Tagawa et al. (2019) ran N-body simulations combined with a semi-analytical model of a Thompson et al. (2005) AGN disk and found a sBH merger rate of ∼0.02 -60 Gpc −3 yr −1 , with binary formation due mainly to dynamical friction and dynamical interactions. The upper limits of these two papers are similar because both papers used a broadly similar disk model, and made similar assumptions about the lifetime of the disk. Tagawa et al. (2019) is different from Paper I and this paper, though, in that they only simulated the outer regions of their AGN disk, far beyond the location of the migration trap for a Thompson et al. (2005) disk. Yang et al. (2019) carried out Monte Carlo simulations of sBHs ground down into alignment with an AGN disk for a variety of supermassive black hole (SMBH) masses and accretion rates and found a BBH merger rate of ∼4 Gpc −3 yr −1 . However, it is difficult to compare their rate with Paper I, because they did not include the sBHs initially co-planar with the disk when it formed.
In this paper, we use the simulation developed in Paper I of the inner disk region surrounding the migration trap of a Sirko & Goodman (2003) disk to further investigate the potential of an AGN disk to provide a significant merger channel for BBH mergers detectable by aLIGO and Virgo. First we examine the impact of sBHs whose inclinations are ground down by the disk until they are co-planar with the disk (McKernan et al. 2014, as well as sBHs that are migrating inward from beyond 1000 AU (see §2.2 and §3.1). Sándor et al. (2011) found that in a proto-planetary disk adding migrating Mars-mass bodies helped disrupt resonant convoys and led to mergers and the formation of Super-Earths. Since the mass ratio of planets to stars in a proto-planetary disk is similar to the mass ratio of sBHs to the SMBH in an AGN disk, we wish to investigate if adding our own additional orbiters will produce similar effects. Second, we use our models from Paper I to predict what fraction of BBHs will orbit around their center of mass in the prograde versus retrograde direction relative to the direction of rotation of the gas disk (see §2.3 and §3.2). This prediction has implications for the spins of the sBHs resulting from the mergers of these BBHs, and therefore the value of the dimensionless aligned spin, χ eff (McKernan et al. , 2019a, which will help constrain the contribution of mergers of sBHs in AGN disks to aLIGO and Virgo detections. Paper I used the unsaturated torque model from , which only considered fully unsaturated, static torques. Including saturation, as in , has been shown to lead to more rapid inward migration for orbiters outside a narrow mass range, and also introduce a mass dependency to both the location and existence of migration traps (Hellary & Nelson 2011;Coleman & Nelson 2014;Dittkrist et al. 2014). Additionally, Paardekooper (2014) and Pierens (2015) found that co-rotation torques on orbiters in disks can often be dynamic, which can act to stall inward migration and boost outward migration.
However, unsaturated, static torques are a reasonable approximation for the torques present in the inner radiative region of an AGN disk, which should be turbulent and have a high viscosity. Nelson & Papaloizou (2004), Baruteau & Lin (2010), Uribe et al. (2011), and Baruteau et al. (2011) all found that co-rotational torques in turbulent disks are subject to stochastic turbulent fluctuations that keep the co-rotational torque unsaturated even in locally isothermal simulations. More recent simulations from Guilet et al. (2013), Uribe et al. (2015), and Comins et al. (2016) corroborate these results, although because the width of the co-rotational region for smaller objects has yet to be resolved, it is possible for saturation to occur if turbulent fluctuations are strong enough to wipe out the horseshoe turns of these smaller orbiters.
Dynamical co-rotation torques are also likely not to play a large role in the migration of orbiters in an AGN disk, with a Shakura & Sunyaev (1973) viscosity parameter α 10 −2 , as long as their masses remain below ∼30 M (see Paardekooper 2014, Paper I). We therefore limit the initial mass of an orbiter in our simulations to 30 M . However, the masses of orbiters in our incoming black holes run, described in Sections 2.2 and 3.1, can rapidly build up to over 100 M . These massive orbiters, which are typically formed close to the migration trap, may end up with boosted outward migration. We discuss these massive orbiters further in §4.
Finally, Benítez-Llambay et al. (2015) found that a heating torque results from a proto-planet's accretion luminosity. Originally, heating torque was found to counteract inward migration, but even more importantly heating leads to significant inclination and eccentricity pumping (Masset 2017;Eklund & Masset 2017). The heating torque resulting from sBHs should be quite different than the torque resulting from planets, as sBHs do not have a surface to heat via accretion shocks. Within the Hill sphere, sBHs will be surrounded by accretion disks if local gas is optically thin, or radiation supported envelopes if local gas is optically thick. In either case, the luminosity can heat up the surrounding AGN disk gas, in analogy with the mechanism found by Benítez-Llambay et al. (2015). We defer further investigation of the important question of sBH feedback and its effect on the dynamics of a sBH in a gas disk to future work.
The outline of this paper is as follows. In §2 we describe the methods used in this paper, including a brief outline of the N-body code used ( §2.1), and how the Nbody code has been altered in order to simulate sBHs that migrate into the inner disk from beyond 1000 AU ( §2.2), and to examine the direction of the orbits of BBHs formed in our simulations ( §2.3). In §3 we discuss our results, including the the impact of additional sBHs migrating inward ( §3.1), and the orientations of the orbits of the BBHs in our simulations ( §3.2). We discuss the implications and caveats of our findings in §4.

METHODS
In this Section we begin by briefly outlining the set up of the N-body code used in Paper I ( §2.1). Next we detail how we update our N-body code in order to first, simulate incoming sBHs that are either ground down into the disk or migrate inward from beyond 1000 AU ( §2.2) and second, determine whether the BBHs formed in our simulations orbit in the retrograde or prograde direction relative to the orbit of the gas disk ( §2.3). These updates allow us to further examine how efficient an AGN disk would be at producing the BBH mergers that have been detected by aLIGO and Virgo.

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 a migration force modeled using the analytic prescription of , a dampening force from dynamical friction derived from the timescales for eccentricity and inclination dampening given by Cresswell & Nelson (2008), and a perturbation spectrum that models a turbulent force driven by magnetorotational instability following Laughlin & Bodenheimer (1994) and Ogihara et al. (2007). These forces depend on the local AGN gas disk properties, including temperature, surface density, opacity, and scale height, as well as the gradients of temperature and surface density. A detailed derivation of these forces can be found in Paper I.
The AGN gas disk properties in our models are taken from the Sirko & Goodman (2003) AGN disk model, which is a modified Keplerian viscous disk model (Shakura & Sunyaev 1973), with a high accretion rate fixed at Eddington ratio 0.5. We use a SMBH mass of  (2003) SMBH accretion disk model used in our simulations. 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 parameter Q as a function of Schwarzschild radius Rs. The top axis represents the translation from Schwarzschild radii to parsecs for a 10 8 M SMBH. M = 10 8 M . The total mass of the disk integrated out to 2 × 10 5 AU is 3.7 × 10 7 M . 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.
Our simulations neglect forces exerted by the sBHs on the gas disk aside from those implicitly modeled by the migration torques, the effects of accretion onto the SMBH and the sBHs, and general relativistic effects. We consider a BBH to have formed in our simulations when two conditions are met. The first is that two sBHs must be 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. The second is that 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 at the Hill radius, We refer to these criteria as the standard merger criteria. In our runs examining the ratio of prograde to retrograde binaries, the first standard criterion is tightened to 0.65 R mH or 0.85 R mH to help us examine the orbits of the BBHs (see §2.3).
Due to the complex and poorly understood interactions between BBHs and the gas disk within the Hill sphere, for simplicity we consider a BBH merged as soon as it forms. As a result, in our simulations BBH formation and mergers between two sBHs are equivalent. Given the conditions of our simulations it is likely these BBHs would merge within approximately 10-500 yr of forming (Baruteau & Lin 2010), which we showed in Paper I is a short timescale compared to any dynamical timescales. In Paper I we also showed that a third migrator is unlikely to ionize a BBH in our simulation, which would prevent the merger from taking place. Nonetheless it is possible for orbiters to escape once they meet our standard merger criterion. We discuss this possibility in Section 3.2, where we actually evolve BBHs down to 0.65 R mH in order to examine the direction of the BBH's orbit. Finally, we ignore the role of kicks due to gravitational recoil at BBH merger. Such kicks will typically cause small perturbations of O(10 2 km s −1 ) to orbital velocities of O(10 4 km s −1 )(a/10 3 R g ) −1/2 , so we anticipate that kicks should not significantly change our results, although they may cause prompt electromagnetic counterparts to the mergers discussed here (McKernan et al. 2019b).
Our runs each begin with ten sBHs, as suggested by Antonini (2014), who used the distribution of S-star orbits around Sgr A to estimate that at least 10 3 sBHs reside within 0.1 pc of the SMBH. This estimate is consistent with the population of O(10 4 ) sBHs within 1 pc of Sgr A inferred by Hailey et al. (2018). We assume that sBHs are evenly distributed throughout the disk and estimate that there should be roughly ten sBHs within 1000 AU (≈ 0.005 pc).
We perform three fiducial runs, F10, F20, and F30, each with ten sBHs having uniform masses of 10, 20, and 30 M sBHs, respectively. In these runs the sBH closest to the SMBH has an initial semi-major axis of 500 AU, and each successive sBH is separated by 30 R mH from the one before it. This separation ensures that initially the impact of the sBHs on each other is minimal.
We then run four more realistic mass distributions, LMA, LMB, HMA, and HMB, with Kroupa (2002) initial mass functions and maximum masses of 15 M and 30 M for the LMA/LMB runs and HMA/HMB runs, respectively. We select masses from the Kroupa (2002) initial mass function, by drawing from a Pareto power law probability distribution of sBHs with a probability density where a = 1.35, the scale factor m 0 = 5M , and x is a mass that is drawn from the distribution. In runs LMA and LMB if a mass generated from the Pareto distribution is greater than 15 M then a new one is generated to ensure all sBHs have masses of less than 15 M . The same is done in the HMA and HMB runs for masses greater than 30 M . The initial semi-major axes of the sBHs in these runs are selected randomly from a uniform distribution between 200 and 1000 AU. For all models, the eccentricities and inclinations are chosen randomly from Gaussian distributions with means 0.05 and 0, respectively, and standard deviations 0.02 and 0.05 radians, respectively. The absolute value of the inclination generated is used, and if the eccentricity selected is negative a new one is selected until it is positive. The initial phases and arguments of pericenter are randomized.
For runs LMA, LMB, HMA, and HMB, the distance between all sBHs is calculated using the values generated. If any two sBHs are within 10 AU new initial positions will be generated to prevent two sBHs from being within less than a couple R mH initially. All sBHs in all seven runs are initialized orbiting in the prograde direction around the SMBH. However, an isotropic distribution of orbits would suggest that roughly half of sBHs should initially be orbiting in the retrograde direction when the disk forms. We defer a more thorough investigation of retrograde orbiters to future work, because we expect the differences in orbital evolution between prograde and retrograde orbiters to be non-trivial. For example, preliminary work by Secunda et al. (in prep) suggests that retrograde orbiters in the disk would ex-perience a rapid increase in eccentricity and decrease in semi-major axis that would lead to a build up of orbiters very close to the SMBH, or even coalescence with the SMBH. Additionally Fabj et al. (2020), show that retrograde sBHs orbiting on inclined orbits with respect to the disk, are unlikely to get captured by the disk, reducing the number of co-planar retrograde orbiters relative to prograde orbiters.
The simulations 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). The names and initial masses (or ranges of initial masses) of the sBHs for each of the seven runs are shown in the first two columns of Table  1.

Incoming Black Holes
Our simulated disk extends only to 10 3 AU. However, assuming sBHs are distributed somewhat uniformly in the inner ∼ 0.1 pc of an AGN disk, additional sBHs should migrate inward from beyond 1000 AU over time. Additionally, sBHs on inclined orbits that intersect the gas disk will be ground down into the plane of the disk (Bartos et al. 2017;McKernan et al. 2018).
Preliminary results by Fabj et al. (2020) show that for a Sirko & Goodman (2003) AGN disk, sBHs are preferentially ground down into the disk at small inclination angles (< 15 • ). sBHs with large initial semi-major axes (∼ 10 5 R g ) lose orbital energy and typically end up ∼ 10 3 R g from the SMBH in < 1 Myr for small inclination angles (≤ 10 • ). In several of their simulations, McKernan et al. (2019a) assume a relatively efficient grind-down rate of 100 sBH/Myr at radii < 5000 R g . Coupling a relatively high grind-down rate with a short migration time in the densest and thinnest part of the AGN disk, leads to an efficient re-population of the inner disk at < 1000 R g (McKernan et al. 2019a).
In our simulated disk with a SMBH of 10 8 M 1 R g ∼ 1 AU. Therefore here for simplicity we assume that sBH repopulation of the inner disk (< 10 3 AU) happens at a rate 1/0.1Myr. More massive sBH or BBHs will preferentially migrate in towards the trap. However here we simply assume that a 10 M sBH migrates inwards from 10 3 AU every 100 kyr starting at t = 100 kyr and apply this to all runs from Paper I (see Table 1 in §3.1). The results of these runs are described in §3.1.

Prograde Versus Retrograde Binaries
Whether a BBH is on a prograde or retrograde orbit around its center of mass with respect to its orbit around the SMBH will have implications for the spin of the sBH produced by the BBH merging. Thus far a majority of aLIGO and Virgo detections have been of BBH mergers with low values for their dimensionless aligned spins, χ eff (Abbott et al. 2019). Therefore, if this trend is confirmed by future aLIGO and Virgo detections, knowing the fraction of BBHs formed in our simulations that orbit in the prograde versus the retrograde direction will help determine if an AGN disk is a prominent channel for producing the mergers detected by aLIGO and Virgo.
In order to ascertain whether the BBHs formed in our simulations are on prograde or retrograde orbits about their centers of mass, we change our first standard merger criterion (see §2) to require that two sBHs approach each other within 0.65 R mH , to further track the evolution of the BBHs' orbits. This change allows us to determine the direction of the orbit of BBHs for four out of our seven runs (F10, F20, F30, and HMA), by both examining by eye the direction of orbit for the BBHs, and using a conservation of angular momentum argument (see §3.2).
However, for three runs, LMA, LMB, and HMB, changing the merger distance to 0.65 R mH led to two sBHs forming a close, very rapidly rotating BBH that greatly decreased the time step of that run. This rapidly rotating BBH is likely an aphysical result, because of our neglect of the interactions between the gas and the BBH within 1 R mH . As a result, for these three runs we use a merger distance criterion of 0.85 R mH to prevent this problem from occurring, and are forced to use only conservation of angular momentum to determine the orientation of the BBH's orbit (see §3.2).
In all runs, our second merger criterion, that the relative kinetic energy of two sBHs must be less than their binding energy in order to form a BBH, remains the same. The results of these runs are described in §3.2.

Incoming Black Holes
Incoming sBHs migrating into the inner disk every 100 kyr can either represent sBHs migrating inward from the outer regions of the disk, or sBHs whose orbits have been ground down by dynamical friction from inclined orbits into alignment with the disk. The initial set-up of the runs described here (see §2.1) are identical to the runs performed in Paper I, and therefore the merger histories of each run is identical to the merger histories of the runs from Paper I up until the first incoming sBH migrates in from beyond 1000 AU at 100 kyr. The early migration history of the F30 and HMA runs are shown in the left and right panels, respectively, of Figure 2. In the more massive runs, most if not all of the initial ten sBHs will have either formed BBHs or are on stable resonant Table 1. Models with sBHs migrating in from beyond 1000 AU. Column 1: Name of run; Column 2: initial masses (or range of masses) of bodies in M before mergers and additional migrators begin drifting in; Column 3: the total combined mass of all bodies in the run in M including the 1000 M from the additional migrators; Column 4: the mass of the most massive sBH at the end of the run in M ; Column 5: the time each run takes to form a black hole over 100 M in kyr.

Run
MsBH  Figure 2). In the less massive runs, migrators drifting in will have an affect on the migration history of the initial ten sBHs, as shown in the right panel of Figure 2, which shows the first 500 kyr of the HMA run, where most sBHs are initially less than 10 M . In all runs, the most massive sBH ends up in the migration trap within the first few 100 kyr and from then on typically forms a BBH with all sBHs that are more than ∼ 80 M , making it often the only black hole 80 M . Figure 3 shows the evolution of the semi-major axes (top panel) and eccentricities (bottom panel) of the orbiters for the first 2.6 Myr of the F10 run. After the initial few incoming sBHs reach the inner disk, the merger histories of these runs become somewhat periodic. The sBH on the outermost resonant orbit will continuously form BBHs and merge with sBHs migrating inward from the outer disk. As this outermost resonant orbiter becomes more massive, its semi-major axis will start to oscillate as the collective mass between itself and the sBH in the migration trap becomes great enough for the resonant sBH to feel a strong gravitational perturbation. The semi-major axis of the sBH in the migration trap will typically oscillate less than the semi-major axis of the resonant orbiter, because the sBH in the trap is usually significantly more massive than the resonant orbiter, and the oscillations are inversely proportional to mass. For example, in the top panel of Figure 3 the purple, 60-80 M , resonant orbiter oscillates a lot more than the green, 270 M sBH in the migration trap, before they form a BBH at around 2.5 Myr. In the HMA run, because the initial sBHs are less massive, and therefore take longer to migrate towards the migration trap, the first incoming sBH interacts with other orbiters before they have a chance to reach stable orbits. At 500 kyr the sBH in the migration trap, shown in green, is 80 M , the inner resonant orbiter, shown in purple, is 21 M , the innermost resonant orbiter, shown in pink, is 6 M , and the outer resonant orbiter, shown in pink, which was the first sBH to migrate in from beyond 1000 AU, is 20 M .
When the outermost resonant orbiter builds up enough mass, in this example around 80 M , the gravitational perturbation will be strong enough to perturb the orbiter out of resonance towards the migration trap where it will form a BBH with the sBH in the trap. The next sBH that migrates inward will end up on a resonant orbit until it too becomes massive enough to form a BBH with the sBH in the migration trap.
The eccentricities shown in the bottom panel of Figure 3 similarly follow a periodic evolution. The initial eccentricities of the additional orbiters migrating in from beyond 1000 AU, which can be seen as the single lines extending up beyond 10 −2 , are quickly dampened by the gas by around three orders of magnitude. However, as these orbiters migrate inward they begin to feel gravitational perturbations from the outer resonant orbiter which drives their eccentricities back up to a few times 10 −3 before the two orbiters form a BBH. These constant interactions with incoming orbiters keep the eccentricity of the outer resonant orbiter relatively high, also a few times 10 −3 . As the resonant orbiter becomes more massive the green sBH in the migration trap begins to be perturbed by the resonant orbiter, which also drives its eccentricity up to a few times 10 −3 . Finally, right before the resonant orbiter and the sBH in the trap form a BBH, the eccentricity of the resonant orbiter can reach as high as O(10 −2 ).
The general pattern shown in Figure 3 and described above repeats itself over the course of the entire run somewhat consistently. Occasionally, an incoming sBH may merge with another incoming sBH before it has time to merge with the resonant orbiter, for example the brown and yellow incoming sBHs at around 0.5 Myr and the green and grey sBHs at around 0.75 Myr. These mergers typically occur when either the resonant orbiter or the sBH in the migration trap become massive enough to produce strong gravitational perturbations that accelerate incoming sBHs towards the migration trap faster than they would typically migrate. When two incoming sBHs merge before either merges with the resonant orbiter, this merger can lead to multiple mergers in rapid succession, relative to the slower somewhat periodic mergers that typically occur, or dramatic changes in orbit. Both multiple mergers and changes in orbit can be seen in Figure 4 which shows the last megayear of the migration history of the HMA run, when there are two black holes in the trap, the more massive of which is 940 M . Because the black holes in the migration traps have greater masses at later times, multiple mergers are more common towards the end of our simulation runs. Figure 3. The migration history of the first 2.6 Myr of the F10 run with sBHs migrating into the inner 1000 AU of the disk every 100 kyr and the same notation as in Figure 2. The top panel shows the evolution of the semi-major axes of the sBHs as a function of time. Most of the initial population of the inner 1000 AU of the disk form BBHs and merge within the first 100 kyr, forming a 90 M sBH that orbits in the migration trap around 665 AU (shown in dark green). Then as sBHs drift into the inner disk they form BBHs with the sBH that is on an outer resonant orbit. When the resonant orbiter becomes massive enough (about 90 M ) it merges with the sBH in the migration trap. This pattern repeats itself three times over the first 2.5 Myr. The bottom panel shows the evolution of the eccentricities of the sBHs as a function of time. The thin lines towards the top of the plot show the initial eccentricities of the sBHs migrating inward from the outer disk which are randomly assigned. These initial eccentricities are quickly dampened to under 10 −3 until they are driven up to a few times 10 −3 right before they form BBHs. The eccentricities of the sBHs on outer resonant orbits (in orange or purple) remain over 10 −3 as their eccentricity is constantly being pumped by the additional orbiters migrating inward, with which they form additional BBHs. The dark green orbiter in the migration trap has its eccentricity dampened after each merger, and then gradually driven higher by the gravitational perturbations of the resonant orbiter as it becomes more and more massive until it reaches an eccentricity of over 10 −3 and merges with the resonant orbiter. Figure 4. The last megayear of the HMA run, with the same notation as in Figure 2. Because the more massive sBH in the migration trap, shown in dark green, is 940 M , it produces strong gravitational perturbations that make the migration history less routine than in the F10 run. For example, multiple instances of BBH formation take place in rapid succession at just past 9.5 Myr when the dark green, grey, and light green orbiters all merge, forming a 60 M sBH. sBHs also can end up on different orbits instead of merging. For example, just after 9.8 Myr, the purple sBH ends up in a trojan orbit with the light green sBH after interacting with the incoming purple and orange orbiters. Figure 5. The migration history for the HMA run from 2.2 to 2.6 Myr with sBHs migrating into the inner 1000 AU of the disk every 100 kyr and the same notation as Figure 2. The large number of sBHs on resonant and co-orbital orbits in this run leads to a rapid series of BBH formation set off by the sBHs migrating inward. Figure 6. An example of how incoming sBHs can destabilize orbits, leading to BBH formation, drawn from the F30 run with sBHs migrating into the inner 1000 AU of the disk every 100 kyr and the same notation as Figure 2. The left panel shows how the merger of an incoming sBH with the 70 M sBH on an outer resonant orbit, represented in pink, breaks the inner resonance of the 110 M sBH, represented in green. The resonance is broken when the relative true longitudes of the green sBH, pink sBH, and purple 250 M sBH in the migration trap approach zero. The right panel shows the true longitude of the pink and purple sBHs relative to the inner resonant orbiter shown in green. After the three have merged they form a 360 M sBH, which continues on in the migration trap on a trojan orbit with the orange 30 M sBH.
The evolution of our various runs also becomes less routine than the evolution shown in Figure 3 if many sBHs end up on resonant orbits. The migration history of the HMA run, which at around 2.4 Myr has six resonant orbiters, from 2.2 to 2.6 Myr is show in Figure 5. When there are many resonant orbiters at once, they generate a strong gravitational perturbation when they approach the same true longitudes, that is line up with each other with respect to the SMBH. Much like a more massive resonant orbiter would, these aligned resonant sBHs accelerate incoming sBHs inward at a faster rate than they would typically migrate. The large number of resonant sBHs can also lead to a rapid sequence of mergers as multiple resonant orbits can be destabilized by a single merger event. For example, the orange sBH in Figure 5 accelerates inward rapidly at around 2.42 Myr and forms a BBH with the purple outer resonant orbiter a relatively short time after the green orbiter that had migrated in before it. The purple outer resonant orbiter is now 70 M . When the true longitudes of the outer resonant orbiter aligns with the other orbiters it is able to decouple the 36 and 10 M co-orbital sBHs orbiting at 640 AU. Destabilizing these co-orbital sBHs leads to two additional mergers of sBHs on inner resonant orbits, and finally the merger of the outer purple resonant orbiter and the green orbiter in the migration trap. These five mergers occur within roughly 100 kyr of each other. Figure 6 shows another example of how an incoming sBH can break an inner resonance that had been stable for nearly 1 Myr, from the F30 run. During the initial evolution of the F30 run, before additional sBHs start migrating in (see left panel of Figure 2), two sBHs end up on a trojan orbit in the migration trap at around 11 kyr, with one sBH on an inner resonant orbit. The 250 M sBH in the trap is shown in purple, the 30 M sBH in the trap is shown in orange, and the 80 M sBH on the inner resonant orbit is shown in green.
At about 1 Myr the inner resonance is broken. The left panel of Figure 6 shows that at this time there are still two trojan orbiters in the migration trap and the second sBH to have migrated into the disk, shown in pink, is on an outer resonant orbit. This outer resonant sBH has merged with all of the subsequent sBHs migrating in from outside the disk. These mergers give the pink sBH a mass of 70 M when it merges with a sBH migrating in, shown in purple. The pink sBH is now 80 M , which makes it massive enough that when the differences between the true longitudes of the pur- Figure 7. An example of turbulence destabilizing orbits from the migration history of the HMA run with sBHs migrating into the inner 1000 AU of the disk every 100 kyr. This figure uses the same notation as Figure 2. The turbulent mode leads to first the merger of the orange 80 M sBH with the navy 30 M sBH which had previously been on a trojan orbit with the green sBH in the migration trap. Next, the navy black hole, now 110 M , merges with the green, 820 M black hole in the trap, forming a 930 M black hole. The turbulent mode can be seen in the oscillations of the semi-major axes of the black holes. Figure 8. The distribution of initial semi-major axes in AU of BBHs formed in all of our runs combined. There is a peak around 3 AU from BBHs with total masses around 10-20 M . This peak extends to around 4.5 AU from BBHs with total masses of around 30-90 M . There is a smaller peak centered around 9 AU from around 6.5 to 10 AU from BBHs in the migration trap with total masses of around 350-1000 M .
ple 250 M sBH in the trap, the true longitude of the green inner resonant orbiter, and its own true longitude approach zero, it gravitationally perturbs the green sBH off its resonant orbit towards the trap, where it merges with the purple sBH.
The difference in true longitude between the purple and green sBHs and the pink and green sBHs is shown in the right panel of Figure 6. The difference in true longitude between the green and pink sBHs are just barely past zero when the purple and green sBHs align, causing the resonance between them to be broken and the green sBH to merge with the purple sBH in the migration trap. The resulting 280 M sBH, represented in green, also ends up merging with the pink sBH shortly after. These merged sBHs continue along in a trojan orbit with the less massive, orange, sBH in the migration trap.
While the incoming sBHs in these runs were the most efficient mechanism for breaking resonances in our simulations, an additional mechanism for breaking resonances, which was also seen in the simulations in Paper I, is turbulence. An example of turbulence leading to the disruption of both a resonant and a trojan orbit is shown in in Figure 7. In this figure when a turbulent mode opens up in the vicinity of the migration trap the orange 80 M resonant orbiter merges with the navy, 30 M sBH, which had previously been on a trojan orbit with the green black hole in the migration trap. Next the black hole shown in navy, which is now 110 M , merges with the green black hole in the migration trap, forming a black hole of mass 930 M . The oscillations in semi-major axis from the turbulent mode are clearly visible, for example in the navy orbiter before it leaves the migration trap, and the pink orbiter once it reaches the migration trap. Figure 8 shows the distribution of the initial semimajor axes, in AU, of BBHs formed in all seven of our simulations combined. Because in our simulations two sBHs must be within a mutual Hill radius of each other to form a BBH (see §2.1), the initial semi-major axes of BBHs are related to the total masses of the BBHs. The distribution in Figure 8 shows the three most common classes of BBH formation in our simulations: 1. BBHs from two 10-20 M sBHs, which typically have initial semi-major axes of around 3 AU when they form near the migration trap. These BBHs either form early in the run from sBHs that are initially in the inner 1000 AU of the disk, or when two incoming sBHs form a BBH before they have the chance to merge with the resonant orbiter. The peak in the distribution in Figure 8 is from this first class of BBHs.
2. BBHs from incoming 10 M sBHs forming BBHs with resonant orbiters, which often grow to around 80 M before forming BBHs with the black holes in the migration trap. The BBHs in resonant or-bits form around 680-700 AU from the SMBH, which means they have initial semi-major axes up to around 4.5 AU. This produces the shoulder in the first peak in Figure 8 that extends to 4.5 AU.
3. Resonant orbiters merging with the black hole in the migration trap at 665 AU, which tend to grow in mass from 100-200 M at earlier times, to as massive as 800-900 M at later times. As a result, at earlier times BBHs form with semi-major axes of around 6.5 AU and at later times with semimajor axes as large as 10 AU. This third class leads to the second maximum in the distribution in Figure 8 between 6.5-10 AU.
The last two columns of Table 1 give the maximum masses formed and the time it takes to form a 100 M object for all seven runs. In our runs sBHs typically merge to form sBHs of over 100 M in no more than 1 Myr. These models, which simulate the realistic influx of sBHs being ground down into the disk and migrating inward from beyond 1000 AU, further suggest that over the course of an AGN disk lifetime, black holes of nearly 1000 M can form. These intermediate-mass black holes (IMBHs) can form because nearly all inwardly migrating sBHs form BBHs over these timescales.

Prograde Versus Retrograde Binaries
In our simulations described in §2.3 most sBHs rapidly (relative to our adopted timestep) close the distance between 1 R mH , where we consider them bound, and 0.65 R mH , where we consider them merged. As a result, occasionally it can be difficult to determine by visual inspection whether the orbit around their centers of mass of sBHs in binaries are prograde or retrograde. The left and center panels of Figure 9 show two examples of BBHs unambiguously orbiting in the prograde and retrograde direction, respectively. The right panel of Figure 9 shows an example of an orbit that is more ambiguous. In all three panels the sBHs' orbits are plotted in orange or black until they are within 1 R mH of each other and their binding energy is greater than their relative kinetic energy, i.e. our standard merger criteria, after which they are plotted in purple.
Since the orbits themselves are not always apparent, we also use an angular momentum argument to determine the direction of the BBHs' orbits. The sum of the angular momentum of the two sBHs, where m i , v i , and r i are the mass, velocity, and distance from the SMBH of the two sBHs, respectively, is equal to L CM + L bin . L bin is the angular momentum of the orbit of the two sBHs in a binary around each other, and, where m tot is the sum of the masses of the two bodies and v CM and r CM are the center of mass velocity and distance of the center of mass from the SMBH, respectively. Therefore if we subtract L CM from the sum of the angular momentum of the two sBHs, the difference will be the angular momentum of the orbit of the sBHs around each other. If the difference is positive (negative), the binary must be orbiting in the prograde (retrograde) direction. We used this method to verify the direction of orbit for three ambiguous BBH orbits in the four runs where the merger criterion was set to 0.65 R mH . This method was used exclusively for the three runs in which the merger criterion needed to be set to 0.85 R mH (LMA, LMB, and HMB, shown in italics in Table 2), because the time it took sBHs to move from within 1 R mH to within 0.85 R mH of each other was too short to inspect those orbits by eye.
Our conservation of angular momentum argument is limited by the fact that our N-body calculation includes a source of angular momentum from the migration torques that act on orbiters, causing them to migrate inwards and outwards. However, the differences in angular momentum are calculated for the first few timesteps after two sBHs meet the standard merger criteria, and the timescales of these timesteps are of order hours, whereas migration timescales are on the order of tens of thousands of years, for even the most massive sBHs in these models. Therefore angular momentum will be nearly conserved at least over several N-body timesteps after a BBH is first formed.
There were, however, several cases in which the difference L sum − L CM changed sign between the time when two sBHs satisfied our standard merger criteria, and when they met the updated merger criteria (i.e. 0.65 or 0.85 R mH ). The change in L sum − L CM is a result of a third body approaching close enough to either remove or contribute angular momentum to the binary. The amount of angular momentum exchanged is several orders of magnitude less than the angular momentum of any of the three orbiters, making it unlikely to dramatically harden or soften the BBH (see Paper I). Still, this exchange of angular momentum may be capable of affecting the evolution of the BBH. We defer further study of the interactions between BBHs and other orbiters to future work that better resolves the evolution of a BBH interacting with gas within its Hill sphere. In our study we use the initial value of L sum − L CM to determine the direction of the orbit of BBH. Figure 9. The positions of two sBHs in a BBH with respect to their center of mass. The left panel shows an example of a BBH from the F20 run that is unambiguously orbiting in the prograde direction, the center panel shows an example of a BBH from the F30 run that is unambiguously orbiting in the retrograde direction, and the right panel shows an example of a BBH from the HMA run whose orbit is ambiguous. One sBH is plotted in orange, the other in black, until they are within 1 RmH of each other, at which point they are both plotted in purple. Orbits are considered unambiguous when it is easy to tell from the points plotted in purple the direction that the two sBHs are orbiting around their centers of mass.
For runs where the merger criterion was 0.65 R mH and we could compare our angular momentum calculation with plots of the orbits there was agreement between our calculation and the plots in all but two cases. We have labeled these orbits as "Unknown." The number of BBHs orbiting around their center of mass in the prograde and retrograde directions, as well as the number of BBHs whose orbits were Unknown are shown in Table  2. The LMA, LMB, and HMB runs are shown in italics because these are the runs where we only trace the evolution of the sBHs down to a distance of 0.85 R mH , and so only angular momentum conservation is used to determine the direction of the BBH orbits.
Overall our models found 15 BBHs orbiting in the prograde and 32 BBHs orbiting in the retrograde direction. The total number of BBHs we examined in our simulations is small, only 49, making it difficult to make any statistically significant statement. However, if it were equally likely for prograde and retrograde BBHs to form in an AGN disk the expectation value would be 24.5, which 32 exceeds by roughly 1.5σ. Additionally, while individually the difference between the number of prograde and retrograde BBHs in each run is not statistically significant, in most cases, it is still worth noting that in all runs, with the exception of the HMB run, there were more retrograde BBHs formed than prograde. This difference is most significant for the F30 run, which had 8 retrograde BBHs and 0 prograde. Therefore we can likely rule out that BBHs orbit primarily in the prograde direction around their centers of mass, and at least conservatively say that it is an even split between the number of BBHs orbiting in the prograde and retrograde direction.
Finally, with the exception of the few sBHs which formed a close, very rapidly rotating BBH, changing the relative distance required for two sBHs to be merged from 1 R mH to 0.85 R mH or 0.65 R mH had little qualita-tive affect on the migration history of our various runs. Therefore, our assumption that two sBHs will merge once they are within a mutual Hill radius of each other remains valid down to 0.65 R mH , in most cases. Future work that resolves the relevant gas physics, and follows the evolution of BBHs down to merger boundaries of 0.1-0.01 R mH is needed to further probe the poorly understood evolution of BBHs in gas disks.

DISCUSSION
Our models that incorporate sBHs drifting in from beyond 1000 AU every 100 kyr build up over-massive sBHs of roughly 100 M in less than a megayear. Incoming sBHs frequently form BBHs with sBHs on resonant orbits near the migration trap. Since resonant orbiters near the trap typically have masses in the range 20 − 80 M , BBHs formed by in-migrating (∼ 10 M ) sBHs with this population will tend to be asymmetric in mass ratio (q ∼ 0.1 − 0.5). The ability of our model to produce asymmetric mass BBH mergers makes the envi-ronment of a migration trap in an AGN disk a promising location for the recent aLIGO merger GW190412 (The LIGO Scientific Collaboration & the Virgo Collaboration 2020).
Over 10 Myr our incoming sBH models build up IMBHs of nearly 1000 M . Our simulations produce such massive black holes because incoming sBHs are effective at perturbing other sBHs out of resonant or trojan orbits, leading nearly all sBHs in each simulation to form BBHs with the most massive black hole, located in the migration trap, over time. Therefore the inward migration of sBHs on orbits initially beyond 1000 AU and the grinding down of sBHs on inclined orbits by the gas disk make our model an even more efficient merger channel for aLIGO and Virgo detections than the scenario described in Paper I, where we only considered preexisting objects in the inner disk. We caution that both of these results come from using the static co-rotational torques derived by , but the hundreds of solar mass black holes that form in these runs would likely be subject to a dynamic co-rotation torque that could lead to different orbital evolution for these sBHs (Paardekooper 2014, Paper I). Additionally IMBHs may also become massive enough to open a gap in the disk and therefore be subject to Type II migration (Lin & Papaloizou 1986).
With these caveats in mind, our models provide a mechanism for building 1000 M IMBHs that may be detectable by the Laser Interferometer Space Antenna (LISA). Typically, the most massive black hole is in the migration trap and will form a BBH with any black hole over 80 M . The merger of an 80 M black hole and a black hole of several 100 M should be detectable first by LISA, if relatively nearby (redshift z 1), and then at merger by aLIGO. At binary separations detectable with LISA, mechanisms that harden binaries (e.g. gas hardening or tertiary encounters) could be testable. Occasionally in our simulations BBHs with small mass ratios and total masses of around 500 M form. Such BBHs are at the edge of the detection capabilities of LISA at z ∼ 1. In addition, if the roughly 1000 M IMBHs formed in our simulations undergo Type II migration inward, they could form SMBH-IMBH binaries. The extreme mass ratio inspirals (EMRIs) of these IMBHs should be detectable by LISA, depending on a variety of properties of the system.
A majority of aLIGO and Virgo detections have low values for the dimensionless aligned spin of merged sBHs, where S i and M i are the spin and mass, respectively, of each sBH in the merged BBH, andL is the unit vector in the direction of the binary angular momentum. Therefore if BBH mergers occurring in AGN disks are producing a significant fraction of aLIGO and Virgo detections, the sBHs produced by our model must either have low spins or mis-aligned spins. Whether a BBH is orbiting around its center of mass in the prograde or retrograde direction relative to the direction of its orbit around the SMBH will change the orientation of the spin of the sBH formed from the BBH merger. The resulting spins of BBH that merge in an AGN disk are distributed around 0.7 if the orbit of the BBH is prograde and -0.7 if the orbit is retrograde (Tichy & Marronetti 2008;McKernan et al. 2018;Varma et al. 2019). Therefore, because BBHs produced in our model are conservatively 50% prograde and 50% retrograde, second stage BBHs formed from the products of these earlier BBH mergers will often have mis-aligned spins. Our model could also produce low spin sBHs from the mergers of BBHs orbiting around their centers of mass in the retrograde direction relative to their orbit around the SMBH, which appears to be the preferred direction of orbit in our simulations. The sBHs produced from these mergers will have spins distributed around a ∼ −0.7 which will evolve towards low spin due to gas accretion at an average rate of (τ AGN /40 Myr)(ṁ/Ṁ Edd ), where τ AGN is the AGN lifetime andṁ/Ṁ Edd is the average gas accretion rate as a fraction of the Eddington rate . At sub-Eddington accretion rates these sBHs would spin down to low spin values on timescales on the order of the lifetime of the disk, O(10 Myr). However, at super-Eddington accretion rates, it would be possible for a ∼ −0.7 sBHs to spin down in a few megayears, while sBHs with spins a ∼ +0.7 should be spun up to maximal spins on similar timescales. Thus, the slight excess of BBHs with retrograde orbits in our simulations could also produce sBHs with low spins, and super-Eddington accretion onto sBH embedded in AGN disks is an important subject for future investigation.
Additionally, even with our conservative estimate of 50% prograde and 50% retrograde, it will still be more common to have sBHs with spins distributed around -0.7, because the timescale for retrograde BBHs to merge is a factor of 5 to 10 shorter than the timescale for prograde BBHs to merge (Baruteau et al. 2011). Therefore, a roughly even ratio of prograde to retrograde BBHs would lead to both the formation of BBHs having mis-aligned spins, and an over-abundance of sBHs that would evolve towards low spins.
The models presented in this paper help deepen our understanding of the AGN disk merger channel for aLIGO and Virgo. Periodically adding sBHs from the outer disk into our simulations of the inner disk within 10 3 AU increases both the number of high mass BBHs and the number of mergers in the vicinity of the migration trap. These mergers in the vicinity of the migration trap will often be asymmetric in mass, resembling the recent aLIGO detection GW190412 (The LIGO Scientific Collaboration & the Virgo Collaboration 2020). Additionally, our models illustrate that BBHs that merge in AGN disks will form with both prograde and retrograde spins. The mergers of the secondary sBHs formed from these mergers should have the low χ eff observed in aLIGO detections. Future work is needed to better understand the influence of different torque models on our simulations. Following the evolution of these BBHs down to within 0.1-0.01 R mH and incorporating realis-tic gas physics are important to better understand how BBHs traverse the LISA detection band and enter the aLIGO detection band. Finally, as the orbital evolution of retrograde orbiters becomes better understood (e.g. Secunda et al. (in prep.); Fabj et al. 2020) simulations should incorporate sBHs orbiting in the retrograde direction as well.