The following article is Open access

The Fate of EMRI-IMRI Pairs in Active Galactic Nucleus Accretion Disks: Hydrodynamical and Three-body Simulations

Peng Peng, Alessia Franchini, Matteo Bonetti, Alberto Sesana, and Xian Chen

Published 2025 August 8 • © 2025. The Author(s). Published by the American Astronomical Society.
The Astrophysical Journal, Volume 989, Number 1Citation Peng Peng et al 2025 ApJ 989 122DOI 10.3847/1538-4357/adef42

Download Article PDF

DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

Article metrics

304 Total downloads
0 Video abstract views

Dates

Received 2024 December 2
Revised 2025 June 16
Accepted 2025 July 10
Published 2025 August 8

Unified Astronomy Thesaurus concepts

Active galactic nuclei; Intermediate-mass black holes; Stellar mass black holes; Gravitational waves

Journal RSS

Create or edit your corridor alerts

What are corridors?

Abstract

Extreme-mass-ratio inspirals (EMRIs) and intermediate-mass-ratio inspirals (IMRIs) are important gravitational wave (GW) sources for the Laser Interferometer Space Antenna. It has recently been suggested that EMRIs and IMRIs can both form in the accretion disk of an active galactic nucleus (AGN). Considering the likely encounter between a stellar-mass black hole (sBH) and an intermediate-mass black hole (IMBH) during migration in the AGN disk, P. Peng & X. Chen (Paper I) showed that a gap-opening IMBH can drive a surrounding sBH to migrate synchronously. In this work, we extend the study in Paper I with a more sophisticated model. We first use 3D hydrodynamical simulations to study the coevolution of the disk and the migration of an sBH in the vicinity of an IMBH. We find that the gaseous torque, together with the tidal torque exerted by the IMBH, can drive synchronized migration until ∼10 Schwarzschild radii from the central supermassive black hole. We further use a relativistic three-body code to study the final fate of the sBH in the GW-dominated regime. We find that the sBH can be either captured or kicked out by the IMBH, which will result in either two subsequent IMRIs or an EMRI followed by an IMRI. These events will bring rich information about the formation and evolution of sBHs and IMBHs in AGNs.

Export citation and abstractBibTeX RIS

Previous article in issue

Next article in issue

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. Introduction

Extreme-mass-ratio inspirals (EMRIs) and intermediate-mass-ratio inspirals (IMRIs) are anticipated sources of gravitational waves (GWs) for milli-Hertz (mHz) GW detectors such as the Laser Interferometer Space Antenna (LISA; P. Amaro-Seoane 2018). An EMRI typically involves an SMBH with 10⁶–10⁷M_⊙ and a stellar-mass black hole (sBH) with 10–100 M_⊙, resulting in a mass ratio q < 10⁻⁴. In contrast, IMRIs can form from an intermediate-mass black hole (IMBH) with 10²–10⁴M_⊙ and a stellar-mass compact object, or from an SMBH and an IMBH, with a mass ratio q in the range 10⁻⁴ ≲ q ≲ 10⁻² (P. Amaro-Seoane et al. 2007).⁶ Since the GW-radiation timescale is proportional to q⁻¹, EMRIs and IMRIs could be detectable in the LISA band for years, accumulating 10⁴–10⁵ GW cycles and providing detailed information about the spacetime near a massive black hole (P. Amaro-Seoane et al. 2007).

One important place for EMRI formation is the accretion disk of an active galactic nucleus (AGN), which is known to be a possible breeding ground of sBHs (e.g., D. Syer et al. 1991; P. Artymowicz et al. 1993; L. Šubr & V. Karas 1999; V. Karas & L. Šubr 2001; Y. Levin 2003; J. Goodman & J. C. Tan 2004). Here, sBHs can either form from massive stars in the disk’s outskirts or be captured from the nuclear star cluster (D. Syer et al. 1991; J. Goodman & J. C. Tan 2004; Y.-F. Jiang & J. Goodman 2011; G. Fabj et al. 2020; Y.-X. Chen & D. N. C. Lin 2024; Y. Wang et al. 2024). Once embedded in the disk, sBHs might migrate toward the central SMBH due to hydrodynamical effects such as density wave interactions (P. Artymowicz et al. 1993; Y. Levin 2003), thermal torque (E. Grishin et al. 2024), and head or tail winds in the geometrically thick disk (S. K. Chakrabarti 1993; B. Kocsis et al. 2011; F. J. Sánchez-Salcedo 2020). sBHs reaching close enough to the SMBH will produce EMRIs, with some recent calculations suggesting that these “wet” EMRIs can significantly contribute to the overall cosmic EMRI rates (Z. Pan & H. Yang 2021; Z. Pan et al. 2021; A. Derdzinski & L. Mayer 2023).

IMRIs can also form in the AGN disks, from disk-embedded IMBHs inspiraling into the central SMBH or IMBHs capturing sBHs. The IMBHs, on the one hand, can be brought into the galactic nucleus through galaxy mergers or inspiraling globular clusters (M. Volonteri et al. 2003; A. Mastrobuono-Battisti et al. 2014). Then, these IMBHs could be captured by the accretion disk, similarly to the capture of sBHs (P. B. Ivanov et al. 1999). On the other hand, the IMBHs could also grow from the aforementioned sBHs accumulated in the AGN disk. The growth is due to either the accretion of the surrounding gas (B. McKernan et al. 2011, 2012; H. Tagawa et al. 2021; Y.-X. Chen & D. N. C. Lin 2023), or frequent binary formation (e.g., H. Tagawa et al. 2020; J. Li et al. 2022; T. C. N. Boekholt et al. 2023; C. Rowan et al. 2023) and merger (e.g., B. McKernan et al. 2012; I. Bartos et al. 2017; N. C. Stone et al. 2017; Y.-P. Li et al. 2022b; R. Li & D. Lai 2022) with the other sBHs in the disk. Formation and growth of IMBHs is expected to be particularly efficient at specific radii within the AGN disk, known as the “migration traps” (J. M. Bellovary et al. 2016; A. Secunda et al. 2019, 2020; P. Peng & X. Chen 2021; E. Grishin et al. 2024), which are typically located at tens to hundreds of Schwarzschild radii (R_S) from the central SMBH. Due to the vanishing gaseous torques around those traps, sBHs will accumulate there, form binaries, and finally merge with the assistance of gas and/or a third body. Moreover, the merger remnants will remain in the accretion disk and participate in the next generation of mergers. Through such “hierarchical mergers” (Y. Yang et al. 2019b; D. Gerosa & M. Fishbach 2021), sBHs could gradually grow into the intermediate-mass range. Then, the IMBH can form an IMRI by either migrating inward and fall into the central SMBH or experiencing another merger with the sBHs in the crowded disk.

Within an AGN disk, sBHs and the IMBHs migrate differentially due to the different driving mechanisms, which could induce an encounter between the sBH and the IMBH. For example, an sBH can excite a perturbative gas density wave in the disk and migrate due to the back reaction, a mechanism called Type I migration (P. Goldreich & S. Tremaine 1980; W. R. Ward 1997). In contrast, for the IMBH, the tidal torque on the surrounding gas can be strong enough to open a gap in the disk (P. B. Ivanov et al. 1999; A. Gould & H.-W. Rix 2000). This gap-opening process is similar to what occurs around gas giants in protoplanetary disks (D. N. C. Lin & J. Papaloizou 1979a, 1979b, 1986; P. Artymowicz & S. H. Lubow 1994). Such ga ap-opening IMBH continues to exchange energy and angular momentum with the accretion disk by its tidal force, thus migrating toward the central SMBH (A. Gould & H.-W. Rix 2000; P. J. Armitage & P. Natarajan 2002; K. D. Kanagawa et al. 2018), a process called Type II migration. Type I and Type II migration rates can differ significantly due to the very different mass of the migrating object and different disk structure. So, if one has an sBH and an IMBH, the one with the slower migration rate can be possibly chased up by the faster one, during their migration.

When they encounter, the sBH and the IMBH interact with each other, which may be analogous to the interaction between a gap-opening gas giant and a sub-gas-opening terrestrial planet or planetesimal in a protoplanetary disk (e.g., J. M. Hahn & W. R. Ward 1996; A. Pierens & R. P. Nelson 2008; E. Podlewska & E. Szuszkiewicz 2008). In particular, Paper I (P. Peng & X. Chen 2023) shows that after the encounter, the migration of the sBH will be largely affected by the IMBH. The IMBH can drive the sBH to migrate synchronously with it, i.e., the sBH is forced to migrate inward with the same migration rate as the IMBH. Such a synchronized migration is driven by the strengthened gaseous torque near the gap (P. Peng & X. Chen 2023), as well as the tidal torque of the IMBH (e.g., H. Yang et al. 2019a). The sBH can migrate synchronously with the outside IMBH until they reach a distance of about 10 R_S from the central SMBH, where the sBH and the IMBH can form an EMRI and an IMRI with the central SMBH, respectively (P. Peng & X. Chen 2023). Since the EMRI and the IMRI in this scenario are close to each other, this kind of GW source could induce very different phenomena in the LISA detection compared to a single EMRI or IMRI, because of resonances, ejections, chaotic behavior, and other exotic phenomena.

The fate of the sBH is interesting because it affects the type of GW phenomena that might appear in the LISA detection. However, in Paper I, we could not assess this fate; i.e., whether the sBH is kicked out, it captured by the IMBH, or eventually falls into the SMBH. This was caused by the limited capacity of the N-body code used in Paper I, as well as the semianalytical 1D model used to simulate the disk evolution, which likely misses important evolutionary features related to the gaseous torque on the sBH. So, in this work, we make use of more sophisticated 3D hydrodynamical simulations which can follow the evolution of the disk as well as the migration of the sBH and the IMBH self-consistently. Furthermore, we use an accurate three-body dynamical code, in which general relativistic phenomena (including GW emission) are modeled in the post-Newtonian (PN) formalism (M. Bonetti et al. 2016, 2021), to which we add a custom gaseous drag implementation to explore the final evolution of the sBH.

The paper is organized as follows. In Section 2, we describe the initial condition and setups of the disk, sBH, and the IMBH in the hydrodynamical simulation. We try two different initial conditions of the disk, with and without the part of the disk inside the orbits of the IMBH (which are referred to as InnerDisk and NoInnerDisk), to study the viscous spread and convergence of the inner part of the disk. In Section 3, we present the results of InnerDisk, to show how the migration of the sBH is accelerated by the presence of the IMBH. We also show how the evolution of the inner part of the disk affects the migration of the sBH. In Section 4 we present the results of NoInnerDisk, in which initially there is no inner disk, to show that the sBH can still migrate synchronously with the IMBH even if the inner disk has been partially accreted. Based on the above results, in Section 5 we use the custom three-body dynamical simulation code to explore the final fate of the sBH. We summarize and discuss our results in Section 6.

2. Method

According to the analysis in Paper I (e.g., see Section 2 of Paper I), an IMBH can encounter an sBH inside its orbit due to differential migration, which is then followed by the synchronized migration. We use the term synchronized migration to identify the situation where the ratio between the migration rate and the semimajor axis $\dot{a}/a$ is the same for both the IMBH and the sBH, which implies that the relative change of semimajor axis a/a(t = 0) − 1 remains the same for the two black holes. In this process, the orbital radius of the sBH is kept close to but slightly smaller than the orbital radius of the IMBH. We therefore start our simulations assuming the sBH orbit to be close to the outer IMBH with the aim of studying how the subsequent migration of the sBH is affected by the IMBH and whether synchronized migration is also seen in a more realistic 3D treatment of the problem.

Similarly to Paper I, we will divide the whole evolution into two stages, depending on whether the IMBH is close enough to the central SMBH so that the GW radiation is strong enough to dominate the orbital decay. In the first stage, where the GW radiation is still negligible, we simulate the coevolution of the disk together with the migration of the sBH and the IMBH, by means of 3D smooth particle hydrodynamical (SPH) code PHANTOM (D. J. Price et al. 2018). In the simulations, the disk is modeled with millions of gas particles while the SMBH, IMBH, and the sBH are modeled as three sink particles (M. R. Bate et al. 1995). The SMBH is put at the center, while the sBH and the IMBH are placed on two circular orbits with similar semimajor axis. We will introduce the detailed setups and initial conditions of the hydrodynamical simulations in the next subsections.

In the second stage, where GW radiation starts to dominate the orbital decay, we need accurate three-body dynamical simulations to model the PN evolution of the system and the possibly chaotic motion of the sBH, while introducing the effect of the disk as an effective external dissipative force. The detailed setup for these simulations will be given in the Section 5.

2.1. Hydro Simulations’ Initial Conditions

In the hydrodynamical simulations, we choose two kinds of initial conditions for the gas disk, referred to as InnerDisk and NoInnerDisk hereinafter. The two initial conditions mainly differ at the “inner disk” part, i.e., the part of the disk inside the orbit of the IMBH. The InnerDisk simulations start with a smooth disk, which extends from close to the SMBH to beyond the orbital radius of the IMBH. Conversely, in the NoInnerDisk simulations the gas initially only exists outside the orbit of the IMBH. We refer to this disk with the term “outer disk” hereinafter. Some of the gas is able to flow through the gap carved by the IMBH and form an inner disk. We choose to also simulate this kind of initial condition because we found that in InnerDisk the surface density of the inner disk gradually decreases because the gas accreted by the SMBH is only partially replenished by gas flowing across the IMBH orbit from the outer disk. Since the sBH is located within the inner disk, its migration may change due to a decrease of the inner disk’s surface density. In order to assess the importance of the inner disk in the migration of the sBH, we therefore perform two simulations with only the outer disk in order to set a lower limit on the inner disk surface density, which is determined by the amount of material that is able to flow through the IMBH orbit reaching the inner regions.

In both InnerDisk and NoInnerDisk cases, the initial surface density is described with the power law ${\rm{\Sigma }}\,={{\rm{\Sigma }}}_{0}{(r/{r}_{0})}^{-p}$ , where r is the orbital radius and r₀ = 1 is the distance unit in the simulation and we change p according to the different setup we use, assuming p = 1 for the InnerDisk and p = 1.5 for the NoInnerDisk case. The physical distance corresponding to r₀ is 100 R_S. The unit of time t in the simulations is P₀, i.e., the Keplerian orbital period at r₀. We run two different simulations with an initial inner disk InnerDisk, one with Σ₀ = 10⁻⁴ and the other one using Σ₀ = 10⁻⁵ in units of ${M}_{{\rm{SMBH}}}/{r}_{0}^{2}$ , where M_SMBH is mass of the SMBH. We refer to these two cases as InnerDisk_Σ4 and InnerDisk_Σ5, respectively. In NoInnerDisk we only use Σ₀ = 10⁻⁵, since we want to simulate the condition when the inner disk is largely accreted as a “lower limit” for its surface density. We choose p = 1.5 in the NoInnerDisk case, since the outer disk mass is more concentrated near the orbit of the IMBH, effectively making the inner disk formation faster. The exact value of p does not affect the final result qualitatively since the mass of the inner disk is mainly determined by the viscous evolution. The disk extends from 0.2 to 5.0 times the IMBH orbital radius in InnerDisk and from 1.0 to 2.5 in NoInnerDisk. Again, we choose a smaller extension of the outer disk in NoInnerDisk to focus on the evolution of the inner disk.

The equation of state for the disk is locally isothermal, i.e., the temperature (as well as the sound speed) of the gas is constant with time at a specific location. The sound speed c_s is also described by the power law ${c}_{{\rm{s}}}={c}_{{\rm{s}},{\rm{0}}}{(r/{r}_{0})}^{-0.5}$ , which means the scale height h/r is constant over the disk radial extent. We choose c_s,0 so that h/r = 0.02, to simulate a relatively thin AGN disk. The viscosity of the gas is described by the Navier–Stokes viscosity model (D. J. Price et al. 2018), where we take the shear kinematic viscosity to be proportional to the sound speed and the scale height through the viscous parameter α (N. I. Shakura & R. A. Sunyaev 1973). We choose the viscous parameter of the disk to be α = 0.02 in all simulations, in order to ensure that the IMBH can successfully open a gap in the disk. We increased the number of particles until we were able to match the accretion rate onto a single central object computed from the simulations with the theoretical value for a steady state disk, i.e., $\dot{M}=3\pi \nu {\rm{\Sigma }}$ .

We found the minimum number of particles needed to resolve the chosen disk viscosity and scale height to be N = 3 × 10⁶. We therefore use N = 3 × 10⁶ in the InnerDisk simulations to simulate the disk, while in the simulations without the inner disk NoInnerDisk we increase the number of particles to N = 10⁷ to better resolve the formation of the inner disk from the material that is able to cross the gap. The parameters and initial conditions for different simulation models are summarized in the Table 1.

Table 1. Initial Conditions and Parameters for Different Simulation Models

Name	Σ₀	p	r	N	a_sBH
`InnerDisk_Σ4`	10⁻⁴	1	0.2–5	3 × 10⁶	0.75

`InnerDisk_Σ5`	10⁻⁵	1	0.2–5	3 × 10⁶	0.75

`NoInnerDisk_075`	10⁻⁵	1.5	1–2.5	1 × 10⁷	0.75

`NoInnerDisk_078`	10⁻⁵	1.5	1–2.5	1 × 10⁷	0.78

Note. The first column is the name of different models. InnerDisk_Σ4 and InnerDisk_Σ5 belong to the same category of simulations InnerDisk but have different surface densities. NoInnerDisk_075 and NoInnerDisk_078 belong to the NoInnerDisk category but have different initial r_sBH. The second column Σ₀ is the surface density at r = 1. The third column p is the power index of the initial surface density profile. The fourth column is the initial radius extension of the disk. The fifth column N is the gas particle numbers in the simulations. The sixth column represents the initial orbital radius of the sBH.

Download table as: ASCII Typeset image

2.2. Parameters of the Sink Particles

As mentioned above, the SMBH, the IMBH, and the sBH are modeled as three sink particles with mass 1.0, 10⁻³, and 2 × 10⁻⁵ respectively, in units of the SMBH mass 10⁶M_⊙. Note that the physical mass of the IMBH is 10 times larger than the mass of IMBH assumed in Paper I. This choice is dictated by the fact that it is harder to resolve accretion for a smaller IMBH mass in the simulations.

The initial orbital semimajor axis of the IMBH is a_IMBH = r₀ = 1, corresponding to a physical distance of 100 R_S. Since we do not take the relativistic effects into account and our disk hydrodynamics are scale free, the choice of the distance unit does not affect the results. This means that we can extend the results in Section 3 and 4 to any a_IMBH as long as GW radiation is not important. The sBH is added into the simulation after the shape of the gap opened by the IMBH has reached a steady state in order to avoid the effect of the transient gap-opening process on the sBH migration. The initial orbital semimajor axis of the sBH a_sBH is chosen to be 0.75 for the InnerDisk simulations. For the NoInnerDisk simulations, we choose two initial semimajor axes a_sBH = 0.75, 0.78 to study the effect of the initial location on the possibility of synchronized migration between the IMBH and sBH. We refer to these two different cases with NoInnerDisk_075 and NoInnerDisk_078. The initial location of the sBH is already close to the edge of the gap (the width of gap is around 0.2 in simulation units). The initial orbital eccentricities of the IMBH e_IMBH and sBH e_sBH are zero.

The accretion radius of the SMBH is set to be 0.03, corresponding to 3 R_S. The accretion radius of the IMBH is 0.25 times its Hills radius, which is around 0.025 in simulation units. The IMBH mass does not increase significantly during the simulation with this accretion radius. We performed one test with a smaller softening length, i.e., 0.03 times the IMBH Hill radius, and confirmed that the interaction between the IMBH and the sBH does not change significantly, which means that our choice of softening would not affect the resonance interaction between the IMBH and the sBH. For the sBH, we assume that its mass does not change during the simulations, since we do not resolve the accretion onto the sBH as it would require an extremely small sink radius and therefore a prohibitively high number of particles.

The gravitational force between sinks and gas particles is softened when gas particles get too close to the sinks. For both the SMBH and the IMBH, the softening radius is set to the accretion radius. The sink particles’ orbit will change with time as a result of their interaction with the gas particles. Therefore, the migration of the sBH and the IMBH is simulated self-consistently. We have performed a test with only one sBH and a central SMBH, and found that these simulations can recover the theoretical Type I migration rate (P. Goldreich & S. Tremaine 1980; W. R. Ward 1997). We have also performed tests with only one IMBH and a central SMBH, and the migration rate of IMBH agrees with the scaling relation found in previous simulations (K. D. Kanagawa et al. 2015, 2018).

3. Accelerated Migration and Disk Evolution

In this section, we show the results of the simulations belonging to the InnerDisk set: InnerDisk_Σ4 and InnerDisk_Σ5. The upper panel of Figure 1 shows the column density plots of the InnerDisk_Σ5 simulations at 0 P₀, 10 P₀ and 150 P₀, respectively. We can see the IMBH gradually carves a gap around the orbit. As mentioned in Section 2, we first let the IMBH open a gap in the smooth disk before inserting the sBH in the simulation. It takes around 100 P₀ for the gap structure to reach a quasi-steady shape, at which point the gap has a depth of ∼1.5 dex, as shown in Figure 2. We then place the sBH on a circular orbit very close to the IMBH and near the inner edge of the gap, as shown by the red dotted line in Figure 2. We chose the initial location for the sBH to be a_sBH = 0.75, since an sBH at a smaller separation would just evolve under Type I migration while an sBH that lies even closer to the IMBH would experience a too strong tidal interaction.

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Surface density plots of the disk with the IMBH embedded in it. The x-axis and y-axis range from −2.5 to 2.5. Brighter colors represent higher surface density and the color bar covers 5 orders of magnitude. The upper and lower panels represent the `InnerDisk_Σ5` and `NoInnerDisk` simulation at t = 0, 10, and 150 P₀, respectively.
Download figure:
Standard image High-resolution image

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Disk surface density profile when the gap structure converges (after 100 P₀, normalized by Σ₀). The x-axis and y-axis show the radius in unit of r₀ and the surface density of the disk in logarithmic coordinates The locations of the sBH and the IMBH are marked with the red dotted line and the black dashed line, respectively.
Download figure:
Standard image High-resolution image

We first show the migration rate of the sBH and analyze its driving mechanisms. For comparison, we also simulate the migration of a single sBH in the same initial location, but in an unperturbed disk, i.e., in the absence of the IMBH. Then, we extend the simulation to a slightly longer timescale, to see whether the disk evolves significantly and whether this affects the sBH migration.

The sBH migration rate in InnerDisk_Σ4 is shown in Figure 3, where the blue (orange) solid (dashed) line represents the migration in presence (absence) of the IMBH. By comparing the overall decrease of the blue solid line with the orange dashed line, we can see that the migration of the sBH is ∼4 times faster when the IMBH is also present. So, the existence of the IMBH accelerates the inward migration of the sBH, qualitatively consistent with what we have found in Paper I. We can also see that in presence of an IMBH the semimajor axis of the sBH does not decrease monotonically, but instead oscillates with time as a consequence of the gravitational tidal pull of the IMBH.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** The relative change of the orbital semimajor axis with time for the sBHs, with (blue solid line) and without (orange dashed line) the IMBH in the simulation, for `InnerDisk_Σ4`. The x-axis represents time in units of P₀. The y-axis shows the change of a_sBH relative to the initial value.
Download figure:
Standard image High-resolution image

In the following, we will briefly describe the behavior of the mechanisms involved in the sBH migration: tidal torque, strengthened Type I torque, and interfering density wave. Then, we will analyze the importance of these mechanisms according to the results. Additional details and equations about how these driving mechanisms affect the migration can be found in Paper I and H. Yang & Y.-P. Li (2024).

The tidal torque from the IMBH oscillates with time, as the orbital phase angles of the sBH and the IMBH change. The time averaged effect is a net sBH migration only when the two orbital periods P_sBH and P_IMBH are in one of the mean-motion resonances, i.e., are a simple integer ratio of each other (C. D. Murray & S. F. Dermott 1999). When the orbit of the sBH lies at the mean-motion resonance with the IMBH orbit, the average tidal effect tends to keep the integer ratio of the orbital timescales (i.e., the condition of mean-motion resonance) and to increase the orbital eccentricity of the sBH. Therefore, an IMBH with a faster migration rate tends to speed up the sBH migration in order to keep the orbital timescale ratio. The strength of this effect is proportional to the mass of the IMBH, since the tidal force is proportional to the IMBH mass.

The second effect to consider is the Type I torque from the back reaction of the gas density wave excited by the sBH (P. Goldreich & S. Tremaine 1980; W. R. Ward 1997). Its strength essentially depends on the detailed profile of the disk, in particular on whether there are some asymmetries in the flow around the sBH (S. J. Paardekooper et al. 2011). So, the torque can be stronger when the sBH is located around the edge of the gap, where the surface density is a steep function of the radius. In general, this kind of gas torque drives the sBH to migrate inward, damps the orbital eccentricity, and does not oscillate with time. The strength of this effect is proportional to the mass of the sBH and the surface density of the disk, since a heavier sBH or a more massive disk will induce a stronger density wave.

Besides the strengthened Type I migration torque, the sBH can also be affected by the density wave excited by the IMBH. Indeed, a recent study about protoplanetary disks found that when there are two planets inside the disk, their migration is shaped not only by their own density waves but also by the density waves excited by the companion (H. Yang & Y.-P. Li 2024). This mechanism is known as interfering density wave, and its detailed behavior is rather complex. It can drive the sBH either to migrate inward or outward, depending on the detailed dynamics and disk parameters. When driving the sBH inward, it also damps its orbital eccentricity. Since the interfering density wave is excited by the IMBH and acts on the sBH, its behavior is somewhat similar to the tidal torque, i.e., the strength depends on the IMBH mass and oscillates with time. On the other hand, since it is a gaseous torque, it also depends on the surface density, similar to the Type I torque. For the above reasons, these particular sets of simulations are not suitable for the verification of the scaling relations in the classical migration scenarios.

To understand whether the tidal torque directly from the IMBH or the two gaseous effects that accelerate the sBH migration, we separate the force acting on the sBH in the simulation into two contributions: the force from the gas particles and the force from the sink particles. We integrate in time the torques that the gas particles and the IMBH exert on the sBH separately, and compare them with the orbital angular momentum change of the sBH during the simulation. We note that, since the change of the specific orbital angular momentum ${L}_{{\rm{sBH}}}\propto {a}_{{\rm{sBH}}}^{1/2}{(1-{e}_{{\rm{sBH}}}^{2})}^{1/2}$ also includes the effect of the orbital eccentricity change, we also calculate the evolution of ${a}_{{\rm{sBH}}}^{1/2}$ separately and compare it with the above terms. We perform the same calculations for both our chosen values of disk surface density, i.e., InnerDisk_Σ4 and InnerDisk_Σ5. The results are shown in the Figures 4 and 5, respectively.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Integrated torques and specific orbital angular momentum change of the sBH for `InnerDisk_Σ4`. The dashed red and the solid blue line represent the time integrated torque exerted by gas particles and by the IMBH, respectively. The orange line represents the total evolution of the orbital angular momentum of the sBH, while the black line represents the contribution of the semimajor axis evolution on the angular momentum change. The x-axis represents the time in the unit of P₀. The y-axis shows how large the specific angular momentum has changed, or how large the torque has contributed.
Download figure:
Standard image High-resolution image

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** Same as Figure 4, but for `InnerDisk_Σ5`.
Download figure:
Standard image High-resolution image

In Figure 4, we can see that the black and the orange line, representing the contribution of the semimajor axis change to the angular momentum and the evolution of the angular momentum itself respectively, vary almost identically over time. This means that the change in angular momentum is mostly due to the migration of the sBH (evolution of a_sBH), rather than the eccentricity change. The small discrepancy represents the contribution of eccentricity damping on the angular momentum. This discrepancy decreases with time as the orbital eccentricity is damped by the gas torque. We can then compare the orange line with the blue and dashed red lines, to see the contribution of tidal torque and the gaseous torque on the migration. The angular momentum is wiggling while decreasing, similarly to the results in Figure 3. The decreasing trend matches quite well the red line, which represents the amount of the angular momentum dissipated by the gas. This naturally implies that most of the angular momentum is extracted by gaseous torques. The wiggling part of the orange line matches in amplitude and phase the blue line, meaning that the tidal torque dominates the evolution of a_sBH on short timescales, i.e., within the first ∼10 P₀. Compared to the red line, the overall decrease of the blue line is much smaller, which means the accumulated angular momentum change driven by the tidal torque is much smaller than that driven by the gaseous effect. We can therefore conclude that in InnerDisk_Σ4, the gaseous torque drives the accelerated migration of the sBH.

In Figure 5, we show the results of the same analysis performed on the InnerDisk_Σ5 simulation, in which there is less gas and thus the gaseous torque is weaker. We can see that, compared to InnerDisk_Σ4, the overall angular momentum decreasing rate (orange line), as well as the contribution of gaseous effect the migration (dashed red line), are smaller (notice the different y-axis scale compared to Figure 4). The reason is that the surface density is 10 times smaller in this case, so the migration rate of the IMBH and the sBH, driven by the gaseous effect, are also roughly 10 times smaller. Furthermore, the smaller Σ₀ also results in weaker eccentricity damping by the gas. So, the discrepancy between the black and orange line, which represents the contribution of the eccentricity damping on the angular momentum evolution, hardly decreases in ∼100 P₀. Relatively, the overall decrease of the blue line is more significant, which means the contribution of tidal torque on the migration is more important in this case but still comparable with the gas induced migration.

The results of the two simulations demonstrate that both the gaseous torque and the tidal torque can drive the accelerated migration. The gas torque dominates in case of high disk surface density, while tidal torque plays a major role when the disk surface density is low. We note that here we do not further reduce the surface density of the gas, since it is harder to keep the whole system stable when the surface density is too low. The reason is that the stability relies on an equilibrium between orbital-eccentricity damping by gas and eccentricity excitation by tidal torque. When the strength of the tidal torque is much larger than gaseous torque, the orbital eccentricity of the sBH increases rapidly and the synchronized migration breaks down at an earlier stage. We further analyze the contribution of the two types of the gaseous torque, i.e., the strengthened Type I torque and the interfering density wave, on the accelerated migration. We found that the density wave excited by the IMBH dominates the evolution of the sBH orbit, at least initially. But as the orbital eccentricity is damped by the gas, the strength of the interfering density wave will decrease so that the strengthened Type I torque will dominate.

In Figure 6, we show the long-term migration rate of the sBH compared to that of the IMBH in InnerDisk_Σ4. To show the sBH migration more clearly, we fit the average decreasing trend of a_sBH with an exponential function as ${\dot{a}}_{{\rm{sBH}}}/{a}_{{\rm{sBH}}}=0.0006\,{\rm{\exp }}(-0.013\,t/{P}_{0})$ , where ${\dot{a}}_{{\rm{sBH}}}$ is the change rate of sBH semimajor axis. The comparison of fitting and raw evolution is shown in the inset of Figure 6. In the main figure, the red dashed line and the orange solid line represent the migration rate of the sBH and the IMBH respectively. By comparing the red dashed and the orange solid line, we can see that initially the migration rate of the sBH is higher. However, the migration rate of the sBH decreases significantly faster compared to that of the IMBH. Around t = 230 P₀, the migration rate of the sBH becomes smaller than the IMBH rate. The exact time at which this transition happens may be affected by the fitting of the very complex decay behavior. As time progresses, the migration rate of the sBH is expected to become even smaller than in the case without an IMBH, as suggested by the trend of the red dashed line compared to the black dotted line. A similar but less evident phenomenon also happens in the case of a lower surface density disk, i.e., in our simulation InnerDisk_Σ5.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** Migration rate evolution of the sBH and IMBH. For clarity, we fit the semimajor axis evolution of the sBH with a simple exponential function, which can capture the general decreasing trend, shown by the inset. In the main figure, the red dashed line and the orange solid line represent the migration rate of the sBH and the IMBH, respectively. The black dotted line represents the migration rate of sBH in the case without IMBH, for comparison. In the inset, the blue solid line represents the raw evolution of a_sBH, while the red dashed line represents the fitted decreasing trend of a_sBH. The x-axis shows the time in unit of P₀. The y-axis of the main figure shows the migration rate normalized with semimajor axis, in logarithmic coordinates. The y-axis of inset shows the relative change of a_sBH.
Download figure:
Standard image High-resolution image

The above decreasing of the sBH migration rate results from the fact that the disk evolves significantly, especially the surface density of the inner disk compared to the overall disk profile. In Figure 7, we show the disk surface density profile evolution over 200 P₀ in InnerDisk_Σ4. Different lines represent the profiles at different times. We can see that the surface density of the inner disk decreases by about 0.5 dex, while the surface density around the gap and outer disk decreases more slowly and is essentially constant at r > 1.8. This is a natural result of the accretion of the gas in the inner disk on the central SMBH. In general, the surface density of the inner disk will continue to decrease, although this part of the disk will be partially replenished by the gas that is able to flow through the gap carved by the IMBH.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** Disk surface density profile evolution after the sBH is put into the simulation for `InnerDisk_Σ4`. The blue dotted, orange dashed, and green solid lines represent the surface density profile right after, 100 P₀ after, and 200 P₀ after the sBH is inserted in the simulation, respectively. The locations of the sBH at different times are marked with dots with different colors accordingly. The x-axis shows the radius. The y-axis shows the surface density in logarithmic coordinate.
Download figure:
Standard image High-resolution image

As we have shown in Figures 4 and 5, the inner disk surface density determines the strength of the gaseous torque, which contributes significantly to the migration of the sBH. It is, however, difficult to make predictions on the synchronized migration on longer timescales because the sBH migration rate is continuously decreasing due to the disk’s evolution. In order to place a lower limit on this migration rate, we initialize the same simulation but without the inner disk. This would allow us to measure the migration rate of the sBH due to the sole effect of the gas that is able to flow through the IMBH gap, effectively producing an inner disk. We present our results in the next section.

4. Lower Limit of Inner Disk and Synchronized Migration

Since the inner disk will be gradually accreted, which will then affect the migration of the sBH, it is necessary to investigate whether there will still be enough gas in the inner disk after a long enough time (i.e., of the order of the migration timescale) to sustain the sBH migration and possibly induce a synchronized migration with the IMBH. We therefore ran two simulations without the part of the disk that lies inside the IMBH orbit in order to place a lower limit on the sBH migration rate driven by the material that is able to overcome the IMBH gap. The parameters of the simulations are shown in Table 1. We initialize an outer disk with initial surface density Σ₀ = 10⁻⁵, since we aim at placing a lower limit on the sBH migration rate.

Figure 8 shows the evolution of the disk surface density and in particular how much gas is able to flow through the gap and feed an inner disk. We show the disk column density for the case NoInnerDisk where we only include the IMBH in the lower panels of the Figure 1. The left, middle and right plots represent the snapshot at 0 P₀, 10 P₀, and 150 P₀. We can see the gas in the outer disk gradually flows in and the IMBH gradually carves a gap around its orbit. The mass of the inner disk gradually increases and is clearly separated from the outer disk by the gap. After its formation, the inner disk evolves very little and the profile finally converges to a steady state after roughly 150 P₀.

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** Disk surface density profile evolution for the result of `NoInnerDisk`. The blue, orange, green, and red dashed lines represent the disk profile 10, 20, 50, and 100 P₀ after the start of the simulation. The purple, brown, pink, and gray solid lines represent the disk profile 150, 200, 250, and 300 P₀ after the start of the simulation. The x-axis shows radius. The y-axis shows the surface density in logarithmic coordinate.
Download figure:
Standard image High-resolution image

We note here that we do not replenish the outer disk with a steady inflow of material, i.e., our disk has a finite mass, so the whole disk surface density drops very slowly due to the accretion. Nevertheless, we find the inner disk shape to reach a quasi-steady state equilibrium between the material coming from the outer disk and the material that is accreted onto the SMBH.

We insert the sBH into the simulation after the formation of the IMBH-carved gap, at t = 150 P₀. We explored two different initial positions of the sBH, i.e., r = 0.75 and r = 0.78, and we refer to them as NoInnerDisk075 and NoInnerDisk078. The disk profile evolution is almost the same for these two simulations. We increase the initial sBH separation to r = 0.78 in order to potentially aid the synchronized migration of the two objects in the disk. We here note that the results in NoInnerDisk075 show that the sBH still migrates more slowly than the IMBH, somewhat in contrast with the results of Paper I, since the migration rate is sensitive to the initial choice of the sBH separation. Nevertheless, when the initial separation decreases due to the slower migration of the sBH, the synchronized migration condition can be reached, as we will show in the following.

The results of NoInnerDisk_078 are shown in Figure 9. We computed the average of the sBH migration rate over ∼15 orbits to reduce the wiggling induced by the very close IMBH tidal torque. The dashed blue and the solid orange lines represent the migration of the sBH after the average procedure and the migration of the IMBH, respectively, while the dotted green line represents the migration rate of the sBH induced by the gaseous torque. Comparing the sBH migration rate here ${\dot{a}}_{{\rm{sBH}}}\,\sim -3\times 1{0}^{-6}\,{a}_{{\rm{sBH}}}$ to that in Figure 5 $\,{\dot{a}}_{{\rm{sBH}}}\sim -1\times 1{0}^{-5}\,{a}_{{\rm{sBH}}}$ , we can see that the migration rate here is several times smaller. The reason is that in this case the inner disk is less massive than the case in Figure 5, which results in an ineffective gaseous torque. By comparing the decrease of the dashed blue line and the orange line, we can see the migration rate of the sBH in this case is almost the same as the IMBH. So, the sBH migrates synchronously with the IMBH. Furthermore, in this case, the sBH migration does not slow down a lot, like found in Figure 6, since in NoInnerDisk the whole disk profile converges. So, we can say that in this case the synchronized migration between the sBH and the IMBH can be kept all along the migration.

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** The migration of the sBH and IMBH in `NoInnerDisk_078`, as well as the migration rate of the sBH contributed by the gaseous effect. The dashed blue, dotted green, and the solid orange lines represent the evolution of a_sBH, the contribution due to the gas and the evolution of a_IMBH. The x-axis shows the time in unit of P₀. The y-axis shows the relative change of semimajor axis compared to the initial value.
Download figure:
Standard image High-resolution image

We note that during the synchronized migration, the sBH could be captured into the mean-motion resonance with the IMBH. For example, the ratio a_sBH/a_IMBH remains constant ∼0.78 during the migration of both objects, while a_sBH and a_IMBH evolve significantly. The value of this ratio indicates that the sBH is trapped in the 3:2 resonance with the IMBH. The success of capture is consistent with the mean-motion resonance theory, which requires the migration timescale to be longer than the typical resonance trapping timescale (e.g., C. D. Murray & S. F. Dermott 1999; A. C. Quillen 2006). Given our simulation parameters and considering the first-order resonance, the migration timescale should be longer than about $0.1\,{q}_{{\rm{IMBH}}}^{-4/3}\,{P}_{0}\sim 1000\,{P}_{0}$ to allow resonant capture (e.g., A. C. Quillen 2006). Indeed, the migration timescale of our system is ∼10⁵ P₀, and therefore satisfies this criterion.

In conclusion, with both the assistance of the gaseous torque and the tidal torque, the sBH can migrate synchronously with the IMBH, even when the inner disk has been partially accreted. For the parameters of our simulation, the equilibrium point at which synchronized migration occurs is 0.78 a_IMBH. The exact value for synchronized migration might depend on the exact profile and width of the gap, as well as on the IMBH mass. We caution that it might be still possible to induce synchronized migration, even if the initial separation between the sBH and the IMBH is larger, because the IMBH would migrate faster than the sBH, therefore reducing their relative separation over time. One possible reason why we do not find synchronized migration for an initial separation of r = 0.75 is that we would need to evolve the system for a significantly longer time, which is computationally expensive.

Note that, because of the nature of the torques driving the migration, this conclusion does not depend on the specific physical scale and can be extended to any physical radius along the migration process as long as the GW radiation does not dominate the dynamics, provided that the disk surface density profile has reached a steady state. Indeed, all the related timescales (e.g., migration, viscous, and accretion timescales) will scale with the same factor, which will not change the migration behavior of the sBH. Therefore, once the IMBH has encountered one sBH inside its orbit, it will migrate inward synchronously until the GW radiation starts to be important. We study this later stage of the evolution by means of custom three-body simulations in the next section.

5. Later Evolution in GW-dominated Scheme

The synchronized migration shown in the last section will continue so long as the gaseous and tidal torques are the dominant physical mechanisms governing the dynamics of the system. As the IMBH and sBH orbit shrink, however, GW emission will become more and more efficient and will eventually dominate their orbital evolution around the SMBH. Since the shrinking rate ${\dot{a}}_{{\rm{GW}}}$ is a strong function of the mass, the IMBH will start to inspiral much faster than the sBH, catching up with it. At this point, the orbits of the sBH will be more and more chaotic, according to the tests done in Paper I. Whether the sBH will be kicked out, fall into the SMBH, or be captured by the IMBH, is unclear. So, here we use a three-body numerical integration code, implemented with PN dynamics and a fictitious force mimicking gaseous drag, to explore the final fate of the sBH in the GW-dominated phase and the related detectable phenomena.

5.1. Dynamical Code and Gaseous Force

We use the code developed in M. Bonetti et al. (2016, 2021). The code integrates the equations of motion of a three-body system evolving the PN dynamics up to order 2.5PN, by directly solving the three-body trajectories from the PN Hamiltonian. This allows us to self-consistently calculate the GW radiation between the sBH and the IMBH when they get close to each other, as well as the orbital decay of the sBH and the IMBH around the SMBH. The initial orbital radius of the IMBH is set to a_IMBH = 30 R_S, where the orbital decay timescale due to the GW radiation becomes shorter than the migration timescale driven by the gas. The initial orbital radius of the sBH is a_sBH = 20 R_S, such that a_sBH/a_IMBH = 0.67. This is slightly smaller than the value we get in the last section and this allows a relaxation stage before the synchronized migration. The initial orbital eccentricity of the sBH is e_sBH = 0.05, and for IMBH it is e_IMBH = 0.01, although the exact value does not affect the result qualitatively. We will keep the initial phase angle of the IMBH fixed to 0, while trying different values between 0 and 2π for the sBH, to test the effect of the randomness on the results. In the following, we use R_S as the length unit, and 1 yr as the time unit.

Since there are no gas particles in the dynamical code, we need to implement the gaseous torque on the sBH and the IMBH through the introduction of a fictitious force, i.e., an extra ad hoc force added to the simulation besides the gravity among the three-body system. This fictitious force will produce the same effect of driving migration and damping eccentricity as the “true” gaseous drag we get in the hydrodynamical simulation. So, we use the gaseous torque results presented in the last two sections to calibrate the force strength and how it depends on the orbital elements. We then test the effect of this fictitious force in the dynamical code for the two-body case, i.e., simulation with only the IMBH and SMBH or the sBH and SMBH, to confirm that it reproduces the effect of the gaseous torque.

We keep the above gaseous force model constant for most of the time in the simulation, since we assume that during the evolution covered by the dynamical simulation, the surface density of the inner disk does not change a lot. This assumption is reasonable according to the simulation result in A. Cerioli et al. (2016), which simulates how the inner disk is squeezed by the IMBH when the orbital decay of the IMBH is dominated by GW radiation. They found that the squeezing effect on the disk surface density is negligible until the IMBH reaches r ∼ 5 R_S. In our simulation, most of the interesting dynamical interplay between the IMBH and the sBH occurs at larger distances. So, it is reasonable to assume a constant gaseous force model in our study. Conversely, when the orbit of the sBH is so chaotic that the location is deep into the region of the gap of the IMBH, we decrease the gas force by one dex to model the low gas density in the gap, since the gas surface density inside the gap is about one dex lower than the surface density outside the gap (see Figure 8). Furthermore, when the distance between the sBH and the IMBH is too small (smaller than 0.03 a_IMBH), we will turn off the gaseous force term to avoid any unpredictable spurious interactions, since we do not know how the gaseous force will affect the evolution when the motion of the sBH is dominated by the IMBH gravity.

5.2. Outcome 1: Swap and Ejection

One example of the system’s evolution is shown in Figure 10. The initial sBH phase angle for this case is 0.6π. We can see that initially the sBH is not migrating with the IMBH, since the IMBH is still far away, while the sBH orbital eccentricity is growing due to the perturbation. Around t = 2 yr, since the sBH orbit enters one of the mean-motion resonances (C. D. Murray & S. F. Dermott 1999) with the IMBH (see Section 3 for details), its orbital eccentricity e_sBH suddenly becomes very large, while a_sBH decreases significantly, due to the strong tidal torque effect when the mean-motion resonance occurs. After that, e_sBH starts to slowly decrease, due to the damping of the orbital eccentricity by both the gas and the GW radiation. Around t = 6.5 yr, when the orbit of the IMBH is close to the sBH orbit again, similar phenomena take place. So, we can see that although the sBH orbital evolution is more irregular in the GW-dominated regime, the tidal torque from the IMBH, together with the orbital-eccentricity damping effect, can drive the evolution of a_sBH roughly synchronously with a_IMBH. However, around t = 7.7 yr, the IMBH’s decay by GW radiation is so fast that the tidal torque is no longer strong enough to synchronize the sBH’s evolution. Eventually, the IMBH catches up with the sBH, inducing strong perturbations in its orbit. The orbit of the sBH becomes so chaotic that it intersects with the orbit of the IMBH, swapping the hierarchy of the three-body system. The IMBH then eventually merges with the central SMBH, while the sBH is left behind, moving onto a highly eccentric orbit around the SMBH.

Figure 10. Refer to the following caption and surrounding text. — **Figure 10.** Migration of the sBH and the IMBH in the GW-dominated scheme, for the case with initial sBH phase angle 0.0π. The orange and the blue lines represent the orbital radius of the IMBH and the sBH. The x-axis represents the time in unit of yr. The y-axis shows the orbital radius in unit of R_S.
Download figure:
Standard image High-resolution image

To better understand the behavior of the sBH when it is “kicked away” by the IMBH, in Figure 11(a) we zoom in the orbital radius evolution around t = 7.7 yr. In Figure 11(b), we instead show the relative distance between the sBH and the IMBH dr, verifying the large impact of the IMBH’s gravity on the motion of the sBH. From Figure 11(a), we can see that around one orbital apocenter, the sBH suddenly “departs” from its original orbit, toward the orbit of the IMBH. The reason is that at this point the distance between the sBH and the IMBH becomes less than the Hills radius of the IMBH, which is ${q}_{{\rm{IMBH}}}^{1/3}\,{a}_{{\rm{IMBH}}}=0.1\,{a}_{{\rm{IMBH}}}\sim 1$ (J. G. Hills 1988), so the sBH feels the gravitational pull of the IMBH. This close encounter event corresponds to the sudden drop in the distance down to dr/R_S < 1 in Figure 11(b).

Figure 11. Refer to the following caption and surrounding text. — **Figure 11.** (a): Orbital radius of the sBH and the IMBH in the final stage, for the case with initial sBH phase angle 0.0π. The orange and the blue lines represent the orbital radii of the IMBH and the sBH, respectively. The x-axis represents the time in unit of yr. The y-axis shows the orbital radius in unit of R_S. (b): Distance between the sBH and the IMBH in the binary formation stage, for the case with initial sBH phase angle 0.0π. The x-axis represents the time in unit of yr. The y-axis shows the distance between the sBH and the IMBH in units of R_S, on logarithmic scale.
Download figure:
Standard image High-resolution image

Within the Hills radius of the IMBH, the sBH forms a binary with the IMBH, corresponding to the stage around t = 7.70763–7.70768 yr in Figure 11(a) and Figure 11(b). Due to the tidal force of SMBH, the binary orbit is irregular, with the closest distance around dr ≲ 0.2. This binary quickly breaks down, and the sBH is kicked onto a highly eccentric orbit beyond the IMBH. In fact, this kind of transient stage has been broadly studied in planetary dynamics in the so-called “Jacobi capture” problem, in which two planets encounter with each other, with typical impact parameter ∼Hills radius, and then form a binary with orbit inside the Hills sphere (e.g., T. C. N. Boekholt et al. 2023; M. Dodici & S. Tremaine 2024). The binary orbit in the Hills sphere can be chaotic due to the nature of such a three-body system, without a stable orbital apocenter and pericenter distance. Due to the chaotic nature and the instability of such a system, if there is no dissipation on the binary’s orbital motion, the system eventually breaks down, kicking away one of the planets. The exact lifetime of the binary is also subject to randomness. In our case, the orbital energy of the binary is partially dissipated by the GWs emitted by the sBH-IMBH pair. But the distance of the sBH-IMBH pair is still large enough that the effect on the final outcome of the encounter is negligible, and the sBH escapes from the Hills sphere due to the instability. During the chaotic interaction, the kinetic energy of the IMBH relative to the SMBH is partially transported to the sBH, so after the transient stage the sBH is kicked to an orbit with a much larger orbital radius. In this case, the sBH has experienced a failed Jacobi capture process.

Although the sBH is eventually kicked away, it has already formed a “temporary” EMRI with the central SMBH which is inside the mHz GW band (e.g., LISA, P. Amaro-Seoane 2018), since the orbital radius of the sBH has already migrated to as low as ∼5 R_S. In fact, from Figure 10, we can see that the sBH stays in 5 –10 R_S for about one year. After the sBH is kicked away, the IMBH forms an IMRI, and finally merges with the central SMBH. So, in this a case, although the sBH is finally kicked out, we can still possibly detect a “temporary, failed EMRI” and an IMRI in sequence.

5.3. Outcome 2: Binary Formation

We examine here a case in which the system forms a stable sBH-IMBH binary induced by the Jacobi capture, followed by the coalescence of the IMBH-sBH binary, thus forming a “light IMRI’. The initial sBH phase angle for this case is 0.0π.

Similarly to Figure 11, we show the evolution of the orbital radii of the sBH and the IMBH, before and during the chaotic three-body dynamical interaction, in Figure 12(a). We also show the relative distance between the sBH and the IMBH in Figure 12(b). Similarly to the previous case, once the distance between the sBH to the IMBH is smaller than the Hills radius (around t = 7.5992 yr), the sBH is “diverted” toward the IMBH and forms a binary with it. In this stage, the motion of the sBH is very chaotic, without a stable apocenter and pericenter distance, similar to what we saw in the last section.

Figure 12. Refer to the following caption and surrounding text. — **Figure 12.** (a): Orbital radius of the sBH and the IMBH in the final stage, for the case with initial sBH phase angle 0.6π. The orange and the blue lines represent the orbital radii of the IMBH and the sBH respectively. The x-axis represents the time in unit of yr. The y-axis shows the orbital radius in unit of R_S.(b): Distance between the sBH and the IMBH in the binary formation stage, for the case with initial sBH phase angle 0.6π. The x-axis represents the time in unit of yr. The y-axis shows the distance between the sBH and the IMBH in unit of R_S, on logarithmic scale. The inset shows the zoom in of this figure between t = 7.599464–7.599471 yr.
Download figure:
Standard image High-resolution image

However, in this case, the minimum distance between the sBH and the IMBH is as low as dr ∼ 0.01 (versus ∼0.2 in the previous case). This corresponds to ∼10 times Schwarzschild radius of the IMBH, and the associated GW radiation is strong enough to dissipate the orbital energy within just few encounters. So, after such encounters, the orbit of sBH-IMBH binary becomes stable and we have a successful Jacobi capture process producing a genuine sBH-IMBH binary, i.e., a “light IMRI.”

The distance between the sBH and the IMBH after the stable binary formation is shown in the Figure 12(b) inset. We can see that compared to the motion in the transient stage, the orbit in this stage is more regular, showing stable apocenter and pericenter passages. Due to the very strong GW radiation, the sBH quickly inspirals into the IMBH in just few orbits. Due to restriction of the simulation setup, the code cannot continue to calculate the motion of the IMBH after the merger of the sBH-IMBH binary. But since the IMBH is already very close to the SMBH (∼10 R_S), it will then inspiral into the SMBH within about one year. So, in this case, two successive IMRI events, i.e., the sBH-IMBH binary and IMBH-SMBH binary, with different chirp masses, could be detected by the GW detector within one year.

5.4. Statistical Result

In the above two cases, all the initial conditions are the same, except the initial phase angle of the sBH. Still, the long-term evolution, especially the dynamical behavior of the sBH during and after it is captured in the Hills sphere of the IMBH, is very different. This highlights that the dynamical evolution and the final fate of the sBH is very sensitive to the initial condition due to the chaotic nature of the three-body problem. Specifically, different initial phase angles of the sBH will result in different impact parameters when the sBH first enters into the Hills sphere, which can further induce a very different set of close approaches and lifetime in the transient stage.

To explore the probability of the different outcomes (swap and sBH ejection versus sBH-IMBH binary formation), we perform a set of simulations with different initial phase angles. More specifically, the initial phase angle of the IMBH is kept as 0, while the initial phase angle of the sBH is chosen in 0–2π with interval 0.02π. The other initial conditions, such as the semimajor axis and the eccentricity, are kept the same. We ran a total of 100 simulations, with the final outcomes are summarized in the Figure 13. Strips with different colors and styles represent different kinds of outcome, as labeled. To highlight the different kinds of evolution, in Figure 14 we show the radial evolution for 10 of them, showing the migration of the IMBH and the sBH, as well as the final fate of the sBH. The orange line and the blue line represent the orbital radius of the IMBH and the sBH, respectively. The numbers of different outcomes in Figure 14 are roughly consistent with the statistics for all 100 simulations.

Figure 13. Refer to the following caption and surrounding text. — **Figure 13.** Distribution of different kinds of outcomes for different simulations. The blue, green (with dash), red (with backslash), and orange (with forward slash) strips represent the outcome of ejection, binary formation, sBH still within ∼10 R_S after ejection, and complete EMRI before IMRI.
Download figure:
Standard image High-resolution image

Figure 14. Refer to the following caption and surrounding text. — **Figure 14.** Ten examples of the sBH and IMBH orbital evolution in the GW-dominated scheme, with sBH initial phase angle as 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, and 1.8 π. The orange and blue lines represent the orbital radius of IMBH and sBH, respectively. The x-axis represents the time in unit of yr. The y-axis shows the orbital radius in unit of R_S. The text in every subfigure represents the initial phase angles of the sBH.
Download figure:
Standard image High-resolution image

Considering the whole set of 100 simulations, about 44% of the simulations end with binary formation and merger of the light IMRI, shown by the dashed green strips in Figure 13. They are represented by the case with initial phase angle of 0.0π, 0.8π, 1.4π, and 1.6π in Figure 14. In all of these cases, the light IMRI merges within a few orbits after the binary formation. The light IMRI formation and merger happen when the IMBH migrates to ≲15 R_S, so after the merger the IMBH will fall into the SMBH within ∼3 yr (merger timescale, M. Maggiore 2008). In this kind of situation, the GW detectors can detect two successive IMRI events with different chirp mass (∼10² and 10⁴ M_⊙, respectively). We also point out that the GW waveform of the light IMRI could be very different from an isolated IMRI, since the merger event happens very close to the central SMBH (e.g., A. Torres-Orjuela et al. 2019; D. J. D’Orazio & A. Loeb 2020; X. Zhang & X. Chen 2023; Z. Zhang & X. Chen 2024).

Conversely, in 46% of the all cases, the sBH and the IMBH has not merged finally, shown by the blue strips in Figure 13. The sBH is ejected after a transient binary stage. These events are represented by the cases with initial phase angle as 0.2π, 0.4π, 0.6π, 1.0π, and 1.8π in Figure 14. In most of the ejection cases, one temporary EMRI event will happen before the IMRI, since the orbit of the sBH has already enters the mHz GW band (∼10 R_S) before being kicked out. Then, the IMBH will inspiral into the SMBH within a few years for these cases. Thus, a mHz GW detector could also detect a “temporary failed EMRI” and an IMRI in sequence. In only ∼4% of all cases, the sBH is kicked out too early, so that only a single IMRI will happen. An example is shown by the case of 0.4π in Figure 14, in which the sBH is kicked out around ∼20 R_S.

Finally, there are some exotic cases not included in the poll discussed above. For example, in ∼7% of all simulations, the sBH just falls into the SMBH before being kicked or captured, as shown by the red strips (with backslash) in Figure 13. One of these cases is shown in Figure 14, corresponding to a phase angle of 1.2π. In this kind of case, there will be a complete EMRI event followed by an IMRI event. In another ∼3% of cases, although the sBH is kicked out, the orbital radius of the sBH is still not very large after ejection (e.g., ∼10 R_S). These cases are shown by the orange strips (with forward slash) in Figure 13. So, in these cases there will be an IMRI event followed by an EMRI event.

We have seen that in most cases two GW events can be detected successively, either two IMRIs, EMRI-IMRI, or IMRI-EMRI. In less than half of these cases, the sBH-IMBH binary can form. So, the probability of the sBH-IMBH binary formation is of the order of 0.1, which is roughly consistent with the previous works on the binary formation in the AGN disk. As mentioned earlier, the binary formation process shown in this work can be described by the Jacobi capture scenario, in which an sBH-IMBH binary with irregular motion can form when the sBH enters into the Hills sphere. Whether the binary could become stable and merge depends on whether the energy dissipation due to the GW radiation is strong enough. Since the strength of the GW radiation is severely dependent on the distance between the sBH and the IMBH (M. Maggiore 2008), the final result will be roughly determined by the closest distance between the sBH and the IMBH during the transient stage. To reduce the orbital energy in several binary orbital timescales, the closest distance between the sBH and the IMBH should be as low as ∼0.01 (also ∼0.01 times the Hills sphere) (e.g., J. Li et al. 2022). According to the previous works, the probability that the closest distance is smaller than 0.01 can be ≲0.1 (e.g., T. C. N. Boekholt et al. 2023; M. Dodici & S. Tremaine 2024), which should also be the rough binary formation probability. In our case, the binary formation probability is slightly higher, which may result from the minor contribution on the energy dissipation by the gaseous torque. Additionally, we note that the gaseous force prescription adopted in our analysis may introduce some uncertainty in the calculated binary formation probability. As mentioned earlier, we have calibrated the gaseous force using the results from hydrodynamical simulations, for different orbital radius and eccentricity of sBH. Such mock force should behave well, at least in the early stages of N-body simulation. However, when the motion of the sBH is totally dominated by the gravity of the IMBH (i.e., when the sBH deviates from the original orbit), the force prescription may no longer be accurate, so we turn off the force when the sBH is too close to the IMBH. There could be some caveats by simply turning off the force, since the gas force could assist the merger of sBH and IMBH, which means that the probability of the final IMBH-sBH binary formation and merger could be affected.

6. Conclusion

In this work, we extended the work in Paper I by investigating the interaction between a gap-opening IMBH and a surrounding sBH in the AGN disk. Using 3D hydrodynamical simulations, we found that the migration of sBHs around the inner edge of the gap opened by the IMBH can be significantly accelerated by both gaseous and tidal torques. We further found that, although the inner disk is partially accreted, the gaseous torque together with the tidal torque of IMBH can induce synchronized migration until approximately 10 Schwarzschild radii from the central SMBH.

We further used a relativistic (PN motion) three-body simulation setup to study the final fate of the sBH in the GW-dominated regime and found that in most cases the sBH can be either captured or ejected by the IMBH. In both cases, two successive GW events can happen. For the ejection case, a temporary EMRI will enter the LISA band, followed by an IMRI event. For the capture case, two IMRIs will happen successively: the first forms from the sBH-IMBH binary, the second corresponds to the IMBH inspiralling into the SMBH, with a time interval ≲3 yr.

Although we only run simulations with the IMBH mass m_IMBH = 1000 M_⊙ due to the resolution limitation, our conclusions should also apply to less massive IMBHs, provided m_IMBH is sufficient to open even a shallow gap in the disk. In the initial stage, where GW radiation is not significant, both the gaseous and tidal torques can still drive synchronized migration, even for cases with a smaller gap-opening IMBH. For example, the Type I torque will be effectively strengthened as long as the disk profile is perturbed by the IMBH. In addition, the interfering density wave and tidal torque will always be effective when the orbital timescale ratio between the sBH and the IMBH is approximately an integer. In the second stage, where GW radiation dominates, a smaller IMBH will result in slower migration, generally enhancing the stability of synchronized migration. The ejection or capture process should still occur, but in a region closer to the SMBH, which makes it even easier to produce two successive EMRI/IMRI events. Nevertheless, if the mass of the outer black hole is too small, the synchronized migration could break down before the GW-dominated regime. For example, previous works have investigated a system in which the outer black hole is an sBH, instead of an IMBH. Although synchronized migration can occur initially, it does break down at ∼100 R_S (H. Yang et al. 2019a). Notably, even a relatively low-mass sBH (with mass ratio compared to SMBH ≳5 × 10⁻⁴) can partially open a gap. Therefore, although gap opening is not a necessary condition for synchronized migration, it is a natural outcome when the synchronized migration could extend into the GW-radiation regime.

In this work, we only study the case in which there is only one sBH around the IMBH. However, there could be more than one sBH. For example, in Figure 4 of A. Secunda et al. (2020), we can see that when the IMBH forms in the migration trap, there is more than one sBH around the IMBH. The system remains dynamically stable, unless another sBH migrates from the outside to the migration trap to disturb the whole system. In addition, the distance between the sBHs is larger than their Hill radii, which avoids any chaotic dynamical interaction among these small sBHs. When the IMBH opens a gap and starts to migrate inward, the system could be disturbed. However, the migration timescale is much longer than the typical dynamical timescale, so the process should be quasi-static, and the system should remain dynamically stable, which means the synchronized migration shown in this work could, in principle, still occur. However, since the dynamics of few-body systems could be very complex, we would still need to explore this problem with numerical simulations in order to draw meaningful conclusions on synchronized migration in this regime. This will be the subject of future work.

Estimating the event rate for such successive EMRI/IMRI events is challenging due to many theoretical uncertainties, including the number of sBHs within an AGN accretion disk and the formation channels of IMBHs. We can roughly estimate the event rate as follows, if the IMBHs are produced by successive mergers of the sBHs inside migration traps. Recent theoretical works, without considering migration traps, predict an EMRI event rate in AGNs of 10–10⁴ yr⁻¹ (Z. Pan et al. 2021). But if there is a migration trap in the disk, the sBHs will accumulate in the trap rather than form EMRIs with central SMBH. These sBHs could then form IMBH by repeated merger. Since the formation of a gap-opening IMBH that can leave the migration trap will consume 10 − 100 sBHs roughly, the formation rate of the IMRI, so as the successive EMRI/IMRI events mentioned in this work, would be about 1–10² yr⁻¹ in this scenario.

Acknowledgments

This work is supported by the National Key Research and Development Program of China (grant No. 2021YFC2203002) and the National Natural Science Foundation of China (grant No. 12473037). A.F. acknowledges support provided by the “GW-learn” grant agreement CRSII5 213497 and the Tomalla Foundation.

Footnotes

6
Although the q = 10⁻⁴ threshold is somewhat arbitrary, there are at least two reasons to define EMRIs and IMRIs as separate categories. The first has to do with their astrophysical nature. While EMRIs involve sBHs and SMBHs, which are observationally well-established objects, IMRIs involve the observationally elusive IMBHs. The second has to do with their GW modeling. While EMRI can be (at least to some extent) modeled in the point-particle limit, this is not true for the IMRI, which makes it much more complicated to construct accurate waveforms.

Please wait… references are loading.